A passive seismic reverse-time migration imaging method based on vector wave field decomposition

Abstract

The application discloses a passive seismic reverse-time migration imaging method based on vector wave field decomposition, and the method comprises the following steps: obtaining multi-component passive seismic records and migration velocity models of a target area; separating P waves and S waves from the multi-component passive seismic records while preserving vector attributes, to obtain P wave records and S wave records; calculating travel time information based on the migration velocity models, and extracting direct P waves, direct S waves and corresponding conversion waves from the P wave records and the S wave records according to the travel time information; performing elastic wave field reverse-time extrapolation on the four types of records respectively as boundary conditions to obtain corresponding wave fields; and performing P-S imaging and S-P imaging by using the extrapolated wave fields to obtain underground structure images. Through the record layer vector wave field separation and target wave field screening, the physical correspondence and anti-crosstalk capability of the wave field separation are improved, and the non-physical mode conversion artifacts and source position interference problems caused by the direct reverse transmission of the coupled records in the prior art are solved.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, and in particular to a passive seismic reverse time migration imaging method based on vector wave field decomposition. Background Technology

[0002] Traditional subsurface structure imaging primarily relies on active-source seismic exploration, reconstructing the subsurface structure by artificially generating seismic waves and combining them with migration imaging methods. While this approach offers high imaging accuracy, it typically suffers from high acquisition costs, complex construction processes, and significant environmental disturbances, limiting its application in scenarios such as deep exploration, urban environments, monitoring of unconventional energy development, and large-scale, dense array observations. In contrast, passive seismic imaging utilizes signals from natural earthquakes, microseismic events, or other passive sources, offering advantages such as lower cost, environmental friendliness, and suitability for long-term continuous observation. In recent years, it has become an important development direction in the field of seismic imaging.

[0003] Currently, passive seismic reverse-time migration (RTM) imaging methods are mainly divided into two categories. The first category treats passive sources as pseudo-active sources. Given known parameters such as source location, seismogenic time, and focal mechanism, it performs imaging using the same processing approach as active source RTM. Theoretically, this method can describe wave propagation well, but the imaging quality heavily depends on the accuracy of source parameters. Deviations in source location or focal mechanism can significantly affect the imaging results. The second category does not rely on prior source information and directly utilizes the coherence between direct and converted waves in multi-component seismic records for imaging. This type of method typically uses multi-component seismic records acquired at the surface as boundary conditions for RTM. During RTM, Helmholtz decomposition is used to separate P-waves and S-waves in the coupled elastic wavefield. Cross-correlation or dot-product imaging conditions are then applied to the separated wavefield components. PS imaging is formed by the interference between the direct P-wave and the P-to-S converted wave in the coda, or SP imaging is formed by the interference between the direct S-wave and the S-to-P converted wave in the coda. This type of method does not rely on prior source information and is suitable for dense arrays and complex passive observation conditions. It is currently a commonly used technique in passive seismic reverse time migration imaging.

[0004] However, the second type of method mentioned above still has several shortcomings in practical applications. First, this type of method usually directly performs reverse-time backpropagation on multi-component coupled seismic records. The direct waves, converted waves, and different wave types in the original records are coupled with each other, which can easily lead to non-physical mode conversions during reverse-time extrapolation. This results in crosstalk and artifacts between PS and SP images in the imaging results, seriously affecting the clarity and reliability of subsurface structure imaging. Second, the Helmholtz decomposition method commonly used in existing technologies is based on divergence and curl operators for wavefield separation. In passive source reverse-time migration, it is prone to polarity reversal and may cause amplitude and phase distortion, making the physical meaning of the imaging results unclear and not conducive to the accurate identification of subsurface interface properties. Furthermore, the cross-correlation or dot-product imaging conditions commonly used in existing technologies are easily affected by the polarization direction angle when directly used for interferometric imaging between P-waves and S-waves, leading to attenuation of imaging amplitude. This makes it difficult to truly reflect the interface conversion response intensity and also has certain limitations in physical interpretation. Furthermore, in passive seismic records, direct waves typically have high energy. If direct waves and converted waves cannot be separated and selected in a targeted manner, the imaging results may retain a strong source location response, or even obscure the true subsurface interface information.

[0005] Therefore, how to effectively separate and utilize the coupled P-waves, S-waves, and converted wave fields in multi-component passive seismic records without prior knowledge of the seismic source information, suppress wave field crosstalk and artifacts caused by non-physical mode conversion, and improve the polarity accuracy, amplitude fidelity, and structural imaging clarity of converted wave imaging is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0006] The purpose of this invention is to provide a passive seismic reverse time migration imaging method based on vector wavefield decomposition. By separating the vector wavefield of the recording layer and filtering the target wavefield, the physical correspondence and anti-crosstalk capability of wavefield separation in multi-component passive seismic reverse time migration imaging are improved. This solves the problems of non-physical mode conversion artifacts and source position interference caused by direct back-transmission of coupled records in the prior art.

[0007] To address the aforementioned technical problems, a first aspect of this invention provides a passive seismic reverse-time migration imaging method based on vector wavefield decomposition, comprising the following steps: Step S1: Obtain multi-component passive seismic records and migration velocity models for the target area. The multi-component passive seismic records contain at least two orthogonal particle velocity components. Step S2: Separate the P-wave and S-wave of the multi-component passive seismic record while preserving vector attributes to obtain P-wave and S-wave records; Step S3: Calculate travel time information based on the offset velocity model, and extract direct P-wave and SP converted wave from the P-wave record and direct S-wave and PS converted wave from the S-wave record according to the travel time information; Step S4: Using the direct P-wave record, SP converted wave record, direct S-wave record, and PS converted wave record as boundary conditions, perform reverse time extrapolation of the elastic wave field to obtain the corresponding direct P-wave field, SP converted wave field, direct S-wave field, and PS converted wave field. Step S5: Perform PS imaging using the direct P-wave field and the PS converted wave field, and perform SP imaging using the SP converted wave field and the direct S-wave field to obtain an image of the underground structure of the target area.

[0008] Furthermore, after performing reverse-time extrapolation of the elastic wave field using the direct P-wave record, SP-converted wave record, direct S-wave record, and PS-converted wave record as boundary conditions, the method further includes: Step S41: Separate the P-wave and S-wave in the direct P-wave field, and extract the pure direct P-wave field component from the separation result; Step S42: Separate the PS converted wave field from the P-wave to the S-wave, and extract the pure PS converted wave field component from the separation result; Step S43: Use the pure direct P-wave wave field component and the pure PS converted wave field component as inputs for PS imaging in step S5.

[0009] Furthermore, after separating the PS converted wave field into P-wave and S-wave components and extracting the pure PS converted wave field components from the separation result, the process further includes: Step S421: Perform P-wave and S-wave separation on the direct S-wave field, and extract the pure direct S-wave field component from the separation result; Step S422: Separate the SP converted wave field from the P-wave to the S-wave, and extract the pure SP converted wave field component from the separation result; Step S423: Use the pure SP converted wave field component and the pure direct S-wave field component as inputs for SP imaging in step S5.

[0010] Furthermore, in step S5, PS imaging and SP imaging are performed, specifically by using amplitude correction normalization imaging processing. The amplitude correction and normalization imaging processing includes: Step S51: Retain the sign information of the dot product result of the two wavefield components involved in imaging; Step S52: Reconstruct the amplitude information using the product of the magnitudes of the two wave field components; Step S53: Normalize the illumination terms of the corresponding direct wave wave field components in the two wave field components. Step S54: Based on the symbol information, the reconstructed amplitude information, and the normalization processing result, generate the imaging result.

[0011] Furthermore, the illumination term corresponding to the direct wave field component in step S53 is specifically the autocorrelation of the corresponding direct wave field component.

[0012] Further, the reconstruction of amplitude information in step S54 includes: Step S541: Based on the geometric relationship that the vibration directions of P-wave and S-wave are nearly orthogonal in an isotropic elastic medium, the vector inner product of the two wave field components is separated into a sign term and an amplitude term. Step S542: Extract the amplitude term and use the amplitude term as the amplitude information of the amplitude correction normalization imaging process.

[0013] Further, step S3, which involves extracting the direct P-wave, SP-converted wave, direct S-wave, and PS-converted wave based on the travel time information, includes: Step S311: Based on the offset velocity model, calculate the theoretical arrival times of the direct P-wave, SP-converted wave, direct S-wave, and PS-converted wave respectively; Step S312: Based on the theoretical time, set corresponding extraction time windows for the P-wave record and the S-wave record respectively; Step S313: Within the extraction window, extract the direct P wave and the SP converted wave from the P wave record, and extract the direct S wave and the PS converted wave from the S wave record.

[0014] Further, the extraction of the PS converted wave from the S-wave record in step S3 includes: Step S321: In the S-wave record, a cutoff window is set for the first preset time period after the direct S-wave, the first preset time period corresponding to surface waves or other non-target strong interference wave trains; Step S322: Extract the PS converted wave during the wake period outside the cutoff window.

[0015] Further, in step S4, the direct P-wave record, SP-converted wave record, direct S-wave record, and PS-converted wave record are used as boundary conditions for elastic wave field reverse-time extrapolation, including: Step S401: Using the direct P-wave record and the PS converted wave record as the first imaging branch, perform elastic wave field reverse time extrapolation to obtain the direct P-wave wave field and the PS converted wave field; Step S402: Using the direct S-wave record and the SP converted wave record as the second imaging branch, perform elastic wave field reverse time extrapolation to obtain the direct S-wave wave field and the SP converted wave field; The first imaging branch is used for subsequent PS imaging, and the second imaging branch is used for subsequent SP imaging. The reverse time extrapolation of the two branches is performed independently.

[0016] Further, in step S5, obtaining an image of the underground structure of the target area includes: Step S501: For each seismic source, each event, or each time period, the PS imaging results obtained from the first imaging branch are superimposed to obtain a PS superimposed image; Step S502: For each seismic source, each event, or each time period, the SP imaging results obtained from the second imaging branch are superimposed to obtain an SP superimposed image; Step S503: Based on the PS overlay image and the SP overlay image, obtain the underground structure image of the target area.

[0017] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described passive seismic reverse time migration imaging method based on vector wave field decomposition.

[0018] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described passive seismic reverse time migration imaging method based on vector wave field decomposition.

[0019] The above-described technical solutions of the embodiments of the present invention have the following beneficial technical effects: 1. By separating the P-wave and S-wave of the recording layer vector and pre-screening the direct and converted waves based on travel time information, the coupled seismic records are decomposed into four types of target wavefields with clear physical types and corresponding to the imaging branches before reverse time extrapolation. This cuts off the path of non-physical mode conversion of non-target wavefields during reverse time extrapolation, effectively eliminating wavefield crosstalk and source position interference introduced by direct backpropagation of coupled records in traditional methods, and significantly improving the clarity and reliability of passive seismic reverse time migration imaging results. 2. The solution involves separating the P-waves and S-waves again after the reverse-time extrapolation of the recorded separation and extracting the pure wave field components. This ensures that the wave fields participating in the imaging are strictly corresponding in terms of physical type, eliminating the secondary crosstalk artifacts introduced by the residual coupled wave field in the imaging stage. At the same time, it avoids the polarity reversal and amplitude distortion problems caused by divergence-curl operations in traditional Helmholtz decomposition, ensuring the polarity correctness and physical interpretability of the imaging results. 3. By adopting amplitude correction normalization imaging conditions that preserve symbol information, reconstruct amplitude by multiplying by modulus and perform illumination normalization, the imaging amplitude attenuation and geometric diffusion caused by the angle between the vibration directions of P-wave and S-wave are effectively compensated while maintaining the interface polarity information. This overcomes the problem of insufficient amplitude fidelity in passive source converted wave imaging under traditional dot-matrix imaging conditions, making the imaging results more accurately reflect the true conversion response intensity of the subsurface interface. Attached Figure Description

[0020] Figure 1 This is a flowchart of the passive seismic reverse time migration imaging method based on vector wave field decomposition provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of a passive source earthquake reverse-time migration method based on vector decomposition of P-wave and S-wave records provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the horizontal components of reconstructed wavefield snapshots of different seismic data provided in embodiments of the present invention. Figure 1 ; Figure 4 This is a schematic diagram of the horizontal components of reconstructed wavefield snapshots of different seismic data provided in embodiments of the present invention. Figure 2 ; Figure 5 This is a schematic diagram of the horizontal components of reconstructed wavefield snapshots of different seismic data provided in embodiments of the present invention. Figure 3 ; Figure 6 This is a schematic diagram of the horizontal components of reconstructed wavefield snapshots of different seismic data provided in embodiments of the present invention. Figure 4 ; Figure 7 This is a schematic diagram of the geometric relationship between the propagation direction of P-waves and S-waves and the vibration direction in an isotropic elastic medium provided by an embodiment of the present invention. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments and the accompanying drawings. It should be understood that these descriptions are merely exemplary and not intended to limit the scope of the invention. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0022] To facilitate understanding, the key terms involved in this invention will first be explained. Passive source reverse-time migration imaging without source information refers to a method that directly utilizes passive seismic records for reverse-time extrapolation and structural imaging without prior knowledge of the source location, time of origin, and source mechanism. P-wave / S-wave separation is the process of decomposing the coupled elastic wave field into P-wave and S-wave components. Direct wave is the wave field that propagates directly from the source to the receiver without target interface transformation. Converted wave is the wave field formed by wave mode transformation when seismic waves propagate at the subsurface interface, including PS converted waves and SP converted waves. PS imaging is an image of the subsurface structure formed using the correspondence between direct P-waves and PS converted waves. SP imaging is an image of the subsurface structure formed using the correspondence between direct S-waves and SP converted waves. Amplitude correction and normalization imaging conditions are imaging conditions that correct the amplitude while preserving imaging symbol information and combine this with normalization processing to generate the imaging result.

[0023] The method provided in this invention aims to solve the problems of wavefield crosstalk, non-physical mode conversion artifacts, and amplitude distortion in existing passive seismic reverse-time migration imaging. By performing finer separation of P-waves, S-waves, and their converted wave components in the seismic record, and combining this with amplitude correction and normalization imaging conditions suitable for passive source converted wave imaging, low-crosstalk, high-fidelity imaging of subsurface interfaces and structures is achieved, resulting in imaging results with correct polarity, more balanced amplitude, and clearer physical meaning. Another objective of this invention is to improve the applicability and stability of passive seismic imaging under complex subsurface structural conditions, ensuring that the obtained imaging results derive more from the real interface information carried by the seismic record itself. This provides a more reliable technical means for high-resolution imaging of targets such as faults, fractures, crustal and lithospheric discontinuities.

[0024] The passive seismic reverse-time migration imaging method based on vector wavefield decomposition in this invention takes multi-component passive seismic records and migration velocity models as input. It performs vector wavefield pre-separation and wavefield type screening on the seismic records, and then performs elastic wavefield reverse-time extrapolation and imaging on the screened wavefield components to finally obtain an image of the subsurface structure. The following details each step in the above method in a specific scenario.

[0025] Please refer to Figure 1 and Figure 2 The first aspect of this invention provides a passive seismic reverse-time migration imaging method based on vector wavefield decomposition, comprising the following steps: Step S1: Obtain the multi-component passive seismic record and migration velocity model of the target area. The multi-component passive seismic record contains at least two orthogonal particle velocity components.

[0026] First, the target area refers to the geological region where subsurface structural imaging is required, such as fault zone development areas, unconventional oil and gas reservoir fracturing areas, crustal and lithosphere discontinuity detection areas, or volcanic and geothermal activity monitoring areas. Multi-component passive seismic records can originate from natural seismic signals recorded by broadband seismic arrays, short-period dense arrays, or downhole geophone arrays deployed at or near the surface, or from microseismic events induced during mining, reservoir impoundment, or hydraulic fracturing. These records contain at least two orthogonal particle velocity components, typically an east-west component and a north-south component; in some implementation scenarios, a vertical component may also be included to form a three-component seismic record. The migration velocity model is a spatial distribution description of the propagation velocity of seismic waves in the subsurface medium of the target area, encompassing both P-wave and S-wave propagation velocities. It can be obtained through active-source seismic exploration, background noise tomography, receiver function inversion, or previous research results; this velocity model will be used for subsequent steps in travel time calculation and elastic wave field inverse extrapolation.

[0027] Step S2: Separate the P-wave and S-wave from the multi-component passive seismic record while preserving vector attributes, to obtain the P-wave record and the S-wave record.

[0028] In multi-component seismic records, P-waves and S-waves are often coupled, and the particle velocity component contains both P-wave and S-wave vibration components. To eliminate mutual interference between different wave types before imaging processing, a separation method that preserves wavefield vector properties is used to process the coupled records at the recording level. Specifically, this separation can be achieved by introducing auxiliary wavefield information to progressively reconstruct the particle velocity component associated only with the P-wave from the coupled wavefield, and then removing it from the original coupled wavefield to obtain the S-wave particle velocity component that retains amplitude and phase information. Alternatively, an equivalent vector wavefield decomposition scheme can be used to achieve the above separation. Compared with the traditional Helmholtz decomposition based on divergence and curl operators, this processing method avoids introducing polarity reversal and amplitude distortion during the separation process, providing high-quality input data for subsequent wavefield screening and imaging. After separation, two types of independent records with clearly defined physical properties are obtained: the P-wave record mainly corresponds to the P-wave propagation response in the medium, and the S-wave record mainly corresponds to the S-wave propagation response in the medium.

[0029] Step S3: Calculate travel time information based on the offset velocity model, and extract direct P-waves and SP-converted waves from P-wave records and direct S-waves and PS-converted waves from S-wave records based on the travel time information.

[0030] In passive earthquake observation, the hypocenter is typically located at a certain depth below or around the target area. As the seismic waves radiated from the hypocenter propagate to the surface receiver, some reach the receiver via a direct path, forming direct waves; others undergo wave mode conversion at the interface below, forming converted waves. This step utilizes the migration velocity model and known receiver locations to calculate the theoretical arrival times of direct P-waves, direct S-waves, and converted PS and SP waves (after conversion at the target interface) from the hypocenter to each receiver using ray tracing or solving equations. Based on these theoretical arrival times, a time window centered on the arrival time of the direct P-wave is set in the P-wave record to extract the direct P-wave, and a time window centered on the arrival time of the converted PS wave is set in the coda to extract the SP converted wave. Similarly, a time window centered on the arrival time of the direct S-wave is set in the S-wave record to extract the direct S-wave, and a time window centered on the arrival time of the converted SP wave is set in the coda to extract the PS converted wave. This decomposes the original seismic record into four independent recording components with clear physical meaning, enabling effective differentiation between the target wavefield and non-target wavefield involved in subsequent imaging.

[0031] like Figure 3 , Figure 4 , Figure 5 and Figure 6 As shown, Figure 3 The P-wave obtained from vector decomposition. Figure 4 The S-wave is obtained by using coupled P-wave and S-wave records as virtual sources. Figure 5 This refers to the vector-decomposed P-wave obtained by using the vector-decomposed P-wave record as a virtual source. Figure 6 The vector-decomposed S-wave is obtained by using the vector-decomposed S-wave record as a virtual source. Figures 3-6 This illustrates that when P-wave / S-wave separation is not performed beforehand in the seismic record, non-physical mode transitions can occur near the receiver plane during retrograde extrapolation. These non-physical transitions can be significantly reduced when extrapolating using either the separated P-wave or S-wave records separately.

[0032] Step S4: Using the direct P-wave record, SP-converted wave record, direct S-wave record, and PS-converted wave record as boundary conditions, perform reverse time extrapolation of the elastic wave field to obtain the corresponding direct P-wave field, SP-converted wave field, direct S-wave field, and PS-converted wave field.

[0033] Elastic wavefield reverse-time extrapolation refers to the process of reconstructing the propagation of seismic waves in the subsurface medium by extending the wavefield backward along the time axis using the elastic wave equation or the acoustic approximation equation as the governing equation and the surface receiving records as boundary conditions. In this step, the four types of records are independently used as boundary conditions for reverse-time extrapolation to avoid non-physical mode conversion caused by the coupling of different types of wavefields during the extrapolation process. Specifically, the wavefields extrapolated from the direct P-wave record and the PS converted wave record constitute the base wavefield of the PS imaging branch; the wavefields extrapolated from the direct S-wave record and the SP converted wave record constitute the base wavefield of the SP imaging branch. The extrapolation processes of the two imaging branches are independent of each other, corresponding to the P-wave to S-wave conversion response and the S-wave to P-wave conversion response at the subsurface interface, respectively, ensuring the correspondence between the imaging target and the wavefield used from a physical mechanism perspective.

[0034] Step S5: Perform PS imaging using the direct P-wave field and the PS converted wave field, and perform SP imaging using the SP converted wave field and the direct S-wave field to obtain the underground structure image of the target area.

[0035] The basic principle of imaging is that there is an inherent coherence relationship between the direct wave field and the converted wave field near the actual subsurface interface. By applying appropriate imaging conditions, this coherence relationship can be extracted from the extrapolated wave field and transformed into a spatial distribution image of the subsurface structure. Specifically, PS imaging utilizes the converted wave energy information generated when the incident P-wave is converted into an S-wave at the subsurface interface, while SP imaging utilizes the converted wave energy information generated when the incident S-wave is converted into a P-wave at the subsurface interface. By superimposing the imaging results obtained from the two imaging branches, an image of the subsurface structure of the target area can be obtained. This image can reflect the spatial location and morphology of structures such as faults, fractures, lithological interfaces, and lithospheric discontinuities.

[0036] In this invention, the separation of P-waves and S-waves in seismic records and extrapolated wavefields can be achieved not only through vector decomposition but also through other equivalent decomposition methods that maintain the correspondence between wavefield types. The extraction of direct waves and converted waves can be achieved not only through travel-time-based time windows but also through arrival-time picking, cross-correlation discrimination, or other wavefield identification methods. This invention can perform PS imaging alone, SP imaging alone, or both PS and SP imaging simultaneously followed by joint interpretation. This invention can be applied to the processing of two-dimensional or three-dimensional passive seismic data, and is also applicable to natural earthquakes, microseismic events, and other passive source records. Amplitude correction and normalization imaging conditions can also be achieved using sign preservation and amplitude correction methods with equivalent principles.

[0037] Existing methods directly perform reverse-time backpropagation on multi-component coupled seismic records, separating P-waves from S-waves and imaging during the backpropagation process. In the original records, P-waves, S-waves, and their converted waves are coupled to each other, easily leading to crosstalk between PS and SP waves and artifacts caused by non-physical mode conversions in the imaging results. Existing Helmholtz decomposition methods, based on divergence and curl operators for wavefield separation, are prone to polarity reversal in passive source reverse-time migration, potentially causing amplitude and phase distortion. When existing cross-correlation or dot-product imaging conditions are directly applied to interferometric imaging between P-waves and S-waves, they are easily affected by the polarization direction angle, leading to attenuation of imaging amplitude. In passive seismic records, direct waves typically have high energy; if direct waves and converted waves are not selectively separated, the imaging results may retain strong source location responses, even obscuring true subsurface structural information. This invention effectively solves the above problems by first separating P-waves and S-waves at the recording level and further distinguishing between direct waves and converted waves, then performing vector decomposition on the wave field after reverse time extrapolation to obtain pure wave field components for imaging, and employing amplitude correction to normalize imaging conditions.

[0038] Through the above steps, the present invention completes vector wavefield separation and wavefield type screening at the recording level, cuts off the path of non-physical mode conversion generated by non-target wavefields during the reverse time extrapolation process, effectively eliminates wavefield crosstalk and source position interference caused by direct back-transmission of coupled records in the prior art, and significantly improves the clarity and reliability of passive seismic reverse time migration imaging results.

[0039] Specifically, step S3 extracts the direct P-wave, SP-converted wave, direct S-wave, and PS-converted wave based on the travel time information, including: Step S311: Based on the offset velocity model, calculate the theoretical arrival times of the direct P-wave, SP-converted wave, direct S-wave, and PS-converted wave respectively.

[0040] Under passive seismic observation conditions, the hypocenter is usually located at a certain depth below the target area. For example, in lithospheric discontinuity detection, hypocenters are often distributed in the transition zone between the crust and mantle, while in microseismic monitoring, they are located near fracturing and modified sections. When the hypocenter wavefield propagates upward to the surface receiving array, the propagation paths and velocities of different wave types are different. Direct P-waves propagate along the direct path from the hypocenter to the receiving point at the P-wave velocity, and their theoretical arrival time can be obtained by dividing the straight-line distance from the hypocenter to the receiving point by the P-wave velocity. Direct S-waves propagate along the same direct path at the S-wave velocity, and their theoretical arrival time can be obtained by dividing the straight-line distance from the hypocenter to the receiving point by the S-wave velocity. SP-converted waves are converted from incident S-waves into P-waves at the subsurface target interface. Their propagation path is divided into two segments: an upward S-wave segment propagating at the S-wave velocity and an upward P-wave segment propagating at the P-wave velocity. The theoretical arrival time is the sum of the travel times of the two segments. Correspondingly, the PS-converted wave is transformed from the incident P-wave into an S-wave at the target interface, and the theoretical arrival time is also obtained by superimposing the travel times of the P-band and S-band. In actual calculations, ray tracing or numerical solution methods based on equations can be used to determine the propagation paths and corresponding theoretical arrival times of the four types of wavefields in the migration velocity model. This process combines the propagation characteristics of the wavefield with velocity model information, providing a quantitative basis for the accurate setting of subsequent time windows.

[0041] Step S312: Based on the theoretical arrival time, set the corresponding extraction time windows for P-wave records and S-wave records respectively.

[0042] In P-wave records, a time window with a certain width is set centered on the theoretical arrival time of the direct P-wave at each receiving point to extract the direct P-wave. Simultaneously, a time window is set centered on the theoretical arrival time of the SP converted wave during the corresponding coma period to extract the SP converted wave. In S-wave records, a time window is set accordingly centered on the theoretical arrival time of the direct S-wave to extract the direct S-wave, and a time window is set centered on the theoretical arrival time of the PS converted wave during the coma period to extract the PS converted wave. The time width of the time window is determined comprehensively based on the duration of the source wavelet and the dispersion broadening effect of the subsurface medium on the wavefield, typically taking a range of several principal periods before and after the theoretical arrival time. For example, for natural earthquake signals with a principal frequency of several hertz, the time window width can be set to the order of several seconds; for microseismic signals with a principal frequency of tens of hertz, the time window width is correspondingly shortened. In complex cases where the arrival times of different wavefields are close to each other or overlap, polarization analysis or waveform cross-correlation criteria can be further combined to help determine the time window boundaries.

[0043] Step S313: Within the extraction window, extract the direct P wave and SP converted wave from the P wave record, and extract the direct S wave and PS converted wave from the S wave record.

[0044] The extraction operation extracts the signal segments within the time window from the original recording, forming four independent recording components. After extraction, each recording component contains only wavefield information consistent with its physical type. The direct P-wave recording primarily retains the longitudinal wave propagation response without interface conversion; the direct S-wave recording primarily retains the transverse wave propagation response without interface conversion; the PS-converted wave recording primarily retains the converted wave energy from P-wave to S-wave at the interface; and the SP-converted wave recording primarily retains the converted wave energy from S-wave to P-wave at the interface. This wavefield filtering method based on travel time information ensures that the wavefield processed in the subsequent reverse-time extrapolation step is clearly defined in terms of physical type. It avoids non-physical mode conversions caused by the coupling of different wavetype components during reverse-time extrapolation, laying a solid foundation for the physical correspondence and anti-crosstalk capability of the final imaging results.

[0045] Through the above-mentioned wavefield extraction and time window selection process based on theoretical arrival time, the original coupled seismic record is decomposed into four types of wavefield components with clear physical types and strict correspondence with subsequent imaging branches. This effectively achieves the physical separation of target wavefield and non-target wavefield at the recording level, providing high-quality input data for subsequent inverse time extrapolation and imaging processing.

[0046] Specifically, step S3, extracting the PS converted wave from the S-wave record, includes: Step S321: In the S-wave record, set a cutoff window for the first preset time period after the direct S-wave. The first preset time period corresponds to the surface wave or other non-target strong interference wave train.

[0047] In passive earthquake observation, when seismic waves generated by natural earthquakes or microseismic events propagate to the Earth's surface, in addition to the target body wave phase, surface waves, such as Rayleigh waves and Love waves, are also generated that propagate along or near the surface. Surface waves typically have higher energy and propagate at lower speeds than shear waves, so their arrival time is often later than that of the direct S-wave, forming a large-amplitude, low-frequency, and dispersed wave train during the tailwave period of the S-wave record. Under dense array observation conditions, surface wave trains may last for tens of seconds or even longer, and their intensity is much higher than the weak PS converted wave signal generated at the target interface. Furthermore, in some complex observation environments, strong non-target interference wave trains may also exist in the S-wave record, caused by nearby human activities, industrial machinery vibrations, or fluid transport in wells. The presence of these interference wave trains will seriously affect the identification and extraction quality of PS converted waves. To eliminate such interference, this step sets a cutoff window within a specific time period after the arrival time of the direct S-wave in the S-wave record. The time range of this cutoff window is determined comprehensively based on the dominant frequency band and group velocity characteristics of surface waves in the observation area, as well as the activity periods of known interference sources. For example, for fundamental surface waves commonly observed in regional seismic observations, their arrival time can be estimated using regional velocity models and epicentral distance. The truncation window start point is set slightly before the theoretical arrival time of the surface wave, and the truncation window end point is set at the moment when the energy at its tail decays to the background noise level. For well fluid migration interference in microseismic monitoring scenarios, the truncation window time period can be estimated based on known fracturing operation sequences and fluid migration velocities. Signal segments within the truncation window will be removed in subsequent processing and will no longer participate in the extraction of converted waves.

[0048] Step S322: Extract the PS-converted waveform during the wake period outside the cutoff window.

[0049] After the truncation window is set, the extraction range of the PS converted wave is limited to the region before or after the truncation window in the wake period. The wake period before the truncation window, i.e., the time window after the arrival of the S-wave and before the arrival of the surface wave, is dominated by body wave scattering and interface converted waves, making it the preferred range for extracting PS converted waves. The wake period after the truncation window can also be included in the extraction range when the surface wave energy has attenuated to a level comparable to the converted wave signal. The start and end times of the extraction window are set based on the theoretical arrival time of the PS converted wave, taking into account the arrival error range caused by the uncertainty of the migration velocity model. Typically, the theoretical arrival time is used as the center, with an appropriate width of time window margin before and after. For example, in crustal structure detection applications, the theoretical arrival time of the PS converted wave can be calculated using the velocity model and the source depth. The width of the extraction window can be set to a range of several seconds before and after this arrival time to cover the travel time residual caused by the uncertainty of the velocity model. By using this tailwave period extraction method under truncated window constraints, the PS converted wave signal can be effectively separated from the recording background dominated by surface waves and strong interference, avoiding the introduction of false structural responses by strong interference wave trains in subsequent inverse time extrapolation and imaging.

[0050] By actively setting a cutoff window during the interference wave train period after the direct S-wave and extracting the PS converted wave during the tail wave period outside the cutoff window, the contamination of the converted wave signal by surface waves and other non-target strong interference wave trains is effectively eliminated, so that the PS converted wave record participating in subsequent inverse time extrapolation and imaging has a higher signal-to-noise ratio and a clearer physical correspondence.

[0051] After performing inverse-time extrapolation of the elastic wavefield using four types of seismic records as boundary conditions, in order to further eliminate wave type coupling that may be reintroduced during the extrapolation process due to wavefield propagation, and to ensure that the two wavefield components involved in imaging are strictly corresponding in terms of physical type, this invention continues to perform wavefield re-decomposition and pure component extraction processing on the extrapolated wavefield after inverse-time extrapolation.

[0052] Specifically, after performing reverse-time extrapolation of the elastic wave field using the direct P-wave record, SP converted wave record, direct S-wave record, and PS converted wave record as boundary conditions in step S4, the process also includes: Step S41: Separate the P-wave and S-wave of the direct P-wave field and extract the pure direct P-wave field component from the separation result.

[0053] Although step S2 has already performed vector separation of P-waves and S-waves in the seismic data at the recording level, during the retrograde extrapolation of the elastic wavefield, the wavefield is affected by factors such as medium inhomogeneity, interface reflection, and transmission as it propagates backward in the subsurface medium. During propagation, P-waves may partially convert to S-waves at velocity interfaces, and S-waves may similarly partially convert to P-waves, resulting in a certain degree of wave mode coupling in the extrapolated wavefield. If this coupled wavefield directly participates in imaging, it will introduce residual non-physical mode conversion artifacts into the imaging results. To address this issue, this step performs P-wave and S-wave separation again on the extrapolated P-wave wavefield, using the same vector wavefield separation method as in step S2. From the separation results, a pure direct P-wave wavefield component containing only P-wave vibration components is extracted. This component physically corresponds precisely to the P-wave particle motion without wave mode conversion.

[0054] Step S42: Separate the P-wave and S-wave of the PS converted wave field, and extract the pure PS converted wave field component from the separation result.

[0055] The extrapolated wavefield of PS-converted waves also faces the problem of re-coupling during propagation. In actual observation scenarios, PS-converted waves are converted from incident P-waves at subsurface interfaces. During reverse-time extrapolation, this wavefield propagates from the surface to the subsurface, and wave mode conversion may occur again at each velocity interface along the way, resulting in the extrapolated wavefield containing both P-wave and S-wave components. This step separates the P-waves and S-waves in the extrapolated PS-converted wavefield, extracting a pure PS-converted wavefield component containing only S-wave vibrations. This component accurately corresponds to the P-to-S-converted wave response generated at the target interface. For example, in microseismic hydraulic fracturing fracture monitoring scenarios, the PS-converted wave energy formed after the microseismic signal generated by fracturing is usually weak after conversion at the reservoir interface. If non-converted P-wave components are mixed into the extrapolated wavefield, it will seriously interfere with the imaging accuracy of the fracture response. Through further separation and purification, the S-wave information directly related to the fracture interface conversion can be effectively preserved.

[0056] Step S43: Use the pure direct P-wave wave field component and the pure PS converted wave field component as the input for PS imaging in step S5.

[0057] After the first two steps of separation and purification, the two wavefield components involved in PS imaging achieve a strict correspondence in physical properties—the pure direct P-wave wavefield component corresponds to the wavefield response of the incident P-wave that maintains its waveform during propagation, while the pure PS-converted wavefield component corresponds to the shear wavefield response generated after the same incident P-wave undergoes waveform conversion at the target interface. The coherence between the two is confined to different stages of the same physical conversion process, ensuring that wavefield interference under subsequent imaging conditions occurs only between wavefield components with a true physical conversion relationship.

[0058] By separating the P-waves and S-waves again and extracting the pure wave field components from the inverse-time extrapolated wave field, the wave fields that ultimately participate in imaging are strictly corresponding in terms of physical type. This eliminates the wave type coupling and residual crosstalk introduced again by the non-homogeneity of the medium during extrapolation propagation. At the same time, it avoids the polarity reversal and amplitude distortion problems that are easy to occur in the traditional divergence and curl decomposition during wave field separation, ensuring the polarity correctness and physical interpretability of the PS imaging results.

[0059] After completing the re-separation and pure component extraction of the extrapolated wave field of the PS imaging branch, this invention performs wave field re-decomposition on the extrapolated wave field of the SP imaging branch using the same processing logic, ensuring that the two input components of the SP imaging are also strictly corresponding in physical type.

[0060] After separating the P-wave and S-wave in the PS converted wave field in step S42 and extracting the pure PS converted wave field component from the separation result, the following steps are also included: Step S421: Separate the P-wave and S-wave of the direct S-wave field and extract the pure direct S-wave field component from the separation result.

[0061] Direct S-wave records are primarily characterized by shear wave particle vibrations. After retrograde extrapolation through the elastic wave field, due to the spatial variations in P-wave and S-wave velocities within the subsurface medium, S-waves also undergo partial mode conversion when crossing velocity interfaces. This results in the extrapolated wave field potentially containing a small amount of P-wave components converted from S-waves, in addition to the original S-wave components. This coupling is particularly pronounced near deep tectonic interfaces, such as the Mohorovičić discontinuity between the crust and upper mantle, where both P-wave and S-wave velocities exhibit significant abrupt changes, and the mode conversion effect during wave field propagation is even stronger. This step involves performing P-wave and S-wave separation operations on the direct S-wave extrapolated wave field, extracting the pure direct S-wave wave field component containing only shear wave vibrations from the separation results, and removing the P-wave components mixed in during propagation. This pure component physically corresponds precisely to the wave field response propagating in the form of shear wave particle motion without interface conversion, providing an accurate incident shear wave reference wave field for subsequent SP imaging.

[0062] Step S422: Separate the SP converted wave field from the P-wave to the S-wave, and extract the pure SP converted wave field component from the separation result.

[0063] SP-converted waves are derived from incident S-waves at the subsurface target interface. They propagate as longitudinal waves in the medium to the surface where they are received, with particle vibrations in the same direction as propagation. During retrograde extrapolation, this wavefield propagates in the reverse direction from the surface to the subsurface target interface, again facing the problem of wave mode conversion due to medium inhomogeneity. In seismic fault zone detection scenarios, the physical properties within the fault zone change drastically, and elastic waves undergo frequent wave mode conversions. The extrapolated SP-converted wavefield may contain significant shear wave components. If these are not separated, artifacts unrelated to the actual fault surface response will be generated in subsequent imaging. This step separates the P-waves from the S-waves in the extrapolated SP-converted wavefield, extracting a pure SP-converted wavefield component containing only longitudinal wave vibrations. This component accurately corresponds to the S-to-P-converted wave response generated at the target interface, ensuring the physical purity of the imaging input wavefield.

[0064] Step S423: Use the pure SP converted wave field component and the pure direct S wave field component as inputs for SP imaging in step S5.

[0065] After the above two steps of wavefield re-decomposition and pure component extraction, the two wavefield components involved in SP imaging achieve a strict correspondence in physical properties: the pure direct S-wave wavefield component corresponds to the wavefield response of the incident shear wave that maintains its waveform during propagation, and the pure SP converted wavefield component corresponds to the longitudinal wave field response generated after the same incident shear wave undergoes S-wave to P-wave conversion at the target interface. The coherence relationship between the two is limited to the incident end and the conversion end of the same physical conversion process, ensuring that wavefield interference under subsequent imaging conditions occurs only between wavefield components with a true physical conversion relationship.

[0066] By re-separating and extracting pure components from the extrapolated wave field of the SP imaging branch, the response information of the conversion of transverse waves to longitudinal waves at the target interface can be accurately separated from the wave field background that is recoupled during propagation. This eliminates the residual wave crosstalk introduced by the non-uniformity of the medium during extrapolation, ensuring that the SP imaging results and PS imaging results have consistent high reliability in terms of physical correspondence and polarity correctness.

[0067] In the process of extrapolating the elastic wave field in reverse time by using four types of seismic records as boundary conditions, this invention adopts a structured processing method of independent extrapolation by channel, which assigns the recording components of different physical types to two independent imaging branches to perform extrapolation operations separately.

[0068] In step S4, the direct P-wave record, SP converted wave record, direct S-wave record, and PS converted wave record are used as boundary conditions for the elastic wave field reverse-time extrapolation, including: Step S401: Using the direct P-wave record and the PS converted wave record as the first imaging branch, perform elastic wave field reverse time extrapolation to obtain the direct P-wave field and the PS converted wave field.

[0069] The two types of records processed by the first imaging branch are both physically related to the incidence of P-waves: the direct P-wave record corresponds to the wavefield response of the P-wave radiated from the seismic source that propagates directly to the surface receiver without interface conversion, while the PS-converted wave record corresponds to the wavefield response of the same P-wave that is converted into a S-wave at the target interface and continues to propagate to the surface. Both share the same incident P-wave wavefield in their physical origin, forming a pair of wavefields with an inherent causal relationship. In fault zone detection applications, when the incident P-wave passes through the fault surface, some of its energy is converted into a S-wave. This converted S-wave and the unconverted direct P-wave maintain a deterministic relationship in travel time and phase at the fault surface location. Assigning these two types of records to the same imaging branch for reverse-time extrapolation allows for the complete preservation of this co-source conversion wavefield correspondence in the extrapolated wavefield, providing physically matched wavefield inputs for subsequent PS imaging.

[0070] Step S402: Using the direct S-wave record and SP converted wave record as the second imaging branch, perform reverse time extrapolation of the elastic wave field to obtain the direct S-wave field and SP converted wave field.

[0071] The two types of records processed by the second imaging branch are both physically related to incident shear waves: the direct S-wave record corresponds to the wavefield response of a shear wave radiated from the seismic source that propagates directly to the surface receiver without interface transformation, while the SP converted wave record corresponds to the wavefield response of the same shear wave that is converted into a P-wave at the target interface and continues to propagate to the surface. Both share the same incident shear wave field in terms of physical origin, thus forming a pair of wavefields with an inherent causal relationship. In the scenario of lithospheric discontinuity detection, the teleseismic S-wave undergoes a wavetype conversion when it crosses the velocity interface within the lithosphere. The converted P-wave and the direct shear wave have a definite travel time relationship at the interface. Assigning these two types of records to the same imaging branch for reverse-time extrapolation ensures that the physical correspondence between the incident shear wave and its converted P-wave is fully preserved in the extrapolated wavefield.

[0072] In this system, the first imaging branch is used for subsequent PS imaging, and the second imaging branch is used for subsequent SP imaging. The inverse-time extrapolation of the two branches is performed independently. The inverse-time extrapolation processes of the first and second imaging branches are computationally independent and do not interfere with each other. Each branch uses records of different physical types as boundary conditions, avoiding the non-physical mode conversions caused by the coupling of different waveform components during extrapolation when all records are mixed and backpropagated in existing technologies. This structured processing method of independent extrapolation for each channel ensures, from a computational architecture perspective, that the response information of P-wave to S-wave conversion and S-wave to P-wave conversion at the target interface can be completely reconstructed in their respective independent wavefield channels. This provides a wavefield foundation with clear physical relationships and accurate waveform correspondence for the subsequent imaging processing of the two imaging branches.

[0073] Specifically, in step S5, PS imaging and SP imaging are performed, specifically using amplitude correction normalization imaging processing; amplitude correction normalization imaging processing includes: Step S51: Retain the sign information of the dot product result of the two wavefield components involved in imaging.

[0074] In traditional imaging, the dot product or cross-correlation operation is directly performed on the two wavefield components. The result contains both sign and amplitude information. The sign is determined by the angle between the vibration directions of the two wavefield components; the dot product is positive when the angle is less than 90 degrees and negative when it is greater than 90 degrees. In isotropic elastic media, the particle vibration direction of the P-wave is the same as the propagation direction, while the particle vibration direction of the S-wave is perpendicular to the propagation direction. The vibration directions of the P-wave and the S-wave generated at the interface are nearly orthogonal, which allows the sign of the traditional dot product result to accurately reflect the polarity of the impedance difference between the two media on both sides of the interface. This step retains this sign information separately from the dot product result, providing a basis for determining the interface polarity in the imaging result. This enables the subsequent imaging result to distinguish between positive and negative polarity impedance interfaces, such as distinguishing between discontinuities where velocity increases with depth and those where velocity decreases with depth.

[0075] Step S52: Reconstruct amplitude information by multiplying the magnitudes of the two wave field components.

[0076] In traditional dot product imaging, the imaging amplitude depends not only on the intrinsic energy of the two wavefield components but also on the modulation of the cosine of the angle between their vibration directions. Since the P-wave and S-wave vibration directions are nearly orthogonal, the value of this cosine term is relatively small. Directly using the absolute value of the dot product as amplitude information would result in a significant underestimation of the amplitude, artificially weakening the interface transition response intensity in the imaging results. This step multiplies the magnitudes of the two wavefield components to reconstruct the imaging amplitude information. The magnitude product depends only on the intrinsic energy of the two wavefield components and is unaffected by the angle between their vibration directions, thus more accurately reflecting the converted wave energy at the interface. In complex structural detection scenarios, such as imaging of steeply dipped faults or fracture zones, the local vibration direction relationship between the incident wave and the converted wave near the interface is complex. Spatial variations in the cosine term can introduce spurious amplitude modulation. Reconstructing the amplitude using the magnitude product can eliminate this geometric distortion.

[0077] Step S53: Normalize the illumination terms of the corresponding direct wave field components in the two wave field components.

[0078] In passive source seismic observations, the distribution of seismic sources is typically uneven, with significant differences in illumination intensity at the same imaging point below ground level due to sources at different azimuths and distances. Simultaneously, the reverse-time extrapolated wavefield is affected by geometric diffusion during propagation, resulting in stronger energy in shallower areas and weaker energy in deeper areas. Without normalization, the amplitude of the imaging results will primarily reflect illumination intensity rather than interface conversion efficiency. This step uses the autocorrelation of the corresponding direct wavefield component or other illumination metrics as a normalization factor to normalize the amplitude information reconstructed in step S52, achieving wavefield illumination compensation. This ensures that the response intensity at each imaging point in the final imaging result more evenly reflects the wave mode conversion capability of the interface itself, rather than being dominated by geometric factors such as source distribution and wavefield propagation path.

[0079] Step S54: Generate imaging results based on symbol information, reconstructed amplitude information, and normalization processing results.

[0080] By assigning the symbol information retained in step S51 to the normalized amplitude information, an imaging result that simultaneously contains the correct interface polarity and the true conversion response intensity can be obtained. In practice, the above calculation can be performed point-by-point for each imaging point, and the imaging results corresponding to each seismic source or event can be accumulated to finally form a complete underground structure image of the target area.

[0081] By separating and processing the symbolic information and amplitude information in the imaging operation, reconstructing the amplitude by multiplying the magnitudes and normalizing it with the direct wave illumination term, the amplitude correction normalization imaging processing effectively eliminates the amplitude attenuation caused by the angle between the vibration directions of P-waves and S-waves and the amplitude distortion caused by uneven illumination in traditional point-matrix imaging, while preserving the ability to distinguish interface polarity. This makes the imaging results more accurately reflect the true converted wave response intensity of the subsurface interface.

[0082] Furthermore, the illumination term corresponding to the direct wave field component in step S53 is specifically the autocorrelation of the corresponding direct wave field component. Autocorrelation calculation can provide an energy metric for the wave field component at zero delay, reflecting the illumination intensity of the direct wave at different imaging locations. During seismic wave propagation, the energy distribution of the same direct wave field at different spatial locations varies significantly due to the combined effects of geometric diffusion, medium absorption attenuation, and source radiation directionality. Directly using uncompensated imaging conditions will result in the imaging amplitude being dominated by illumination intensity, masking the interface's conversion efficiency information. Using the autocorrelation of the direct wave field component as a normalization factor allows for a quantitative measurement of the wave field energy at the imaging point. By dividing the reconstructed amplitude information by this illumination metric, amplitude distortion caused by differences in wave field propagation paths and uneven source distribution is effectively eliminated, resulting in a more balanced reflection of the interface's wave mode conversion capability in the final imaging result.

[0083] In step S54, the process of reconstructing amplitude information further utilizes the geometric relationship between the vibration directions of P-waves and S-waves in an isotropic elastic medium, which is achieved by explicitly separating the vector inner product into a sign term and an amplitude term.

[0084] For details, please refer to Figure 7 Step S54 reconstructs the amplitude information, including: Step S541: Based on the geometric relationship that the vibration directions of P-wave and S-wave are nearly orthogonal in an isotropic elastic medium, the vector inner product of the two wave field components is separated into a sign term and an amplitude term.

[0085] In isotropic elastic media, the particle vibration direction of a P-wave is parallel to the wave propagation direction, while the particle vibration direction of an S-wave is perpendicular to the wave propagation direction. Therefore, when co-originating P-waves and S-waves meet at the same imaging point, their particle vibration directions are approximately perpendicular to each other. Traditional dot product imaging directly calculates the vector dot product of the two wavefield components. The result includes both the cosine of the vibration direction angle reflecting the sign of the interface impedance difference and the amplitude information reflecting the magnitude of the wavefield energy. These two aspects are numerically inseparable, causing the imaging amplitude to deviate from the true interface conversion response intensity due to the modulation of the cosine term. This step explicitly decomposes the result of the vector dot product calculation into two independent parts: a sign term and an amplitude term. The sign term retains only the sign of the cosine of the angle, while the amplitude term extracts the energy measure of the two wavefield components themselves. In applications such as lithospheric discontinuity detection or fault zone structure imaging, the wave impedance difference of the target interface may be small, resulting in weak converted wave energy. If the imaging amplitude is further suppressed by the vibration direction angle, the weak interface signal is easily drowned out by background noise. By separating the sign from the amplitude, it can be ensured that the imaging amplitude is no longer affected by the angle cosine modulation.

[0086] Step S542: Extract the amplitude term and use the amplitude term as the amplitude information for amplitude correction normalization imaging processing.

[0087] The amplitude term is constructed by multiplying the magnitudes of the two wavefield components. The magnitude calculation depends only on the particle velocities of each wavefield component and is independent of the angle between their vibration directions. Therefore, using the amplitude term as amplitude information avoids the amplitude attenuation effect caused by the near orthogonality of the P-wave and S-wave vibration directions. This ensures that the response intensity at the interface in the imaging results depends only on the energy coherence of the direct wavefield and the converted wavefield at the imaging point, and is not affected by the numerical value of the geometric angle between their vibration directions. In scenarios such as microseismic crack monitoring or shallow structure detection in urban environments, the imaging target is usually a small-scale structure with a weak reflection coefficient. Using the amplitude term as amplitude information can effectively improve the detectability of weak converted responses and avoid amplitude underestimation caused by geometric angles.

[0088] By explicitly separating the vector inner product into a sign term and an amplitude term, and using the amplitude term as amplitude information, amplitude correction normalization imaging processing maintains the ability to distinguish interface polarity while eliminating the inherent amplitude attenuation effect due to the near orthogonality of the vibration directions of P-waves and S-waves, making the imaging amplitude more realistically reflect the energy intensity of the interface conversion wave.

[0089] After obtaining the PS and SP imaging results of each seismic source, event, or time period, this invention generates an image of the underground structure of the target area by superimposing and jointly applying the multi-source imaging results.

[0090] In addition, step S5 obtains an image of the underground structure of the target area, including: Step S501: The PS imaging results obtained from the first imaging branch are superimposed on each seismic source, event, or time period to obtain a PS superimposed image.

[0091] In passive seismic observation, available sources typically include multiple natural earthquake events or microseismic events at different azimuths and distances. The energy of the PS-converted waves generated by a single event is relatively limited, and the illumination effect on the same target interface varies depending on the source radiation pattern. In long-term observation scenarios with dense arrays, the array may continuously record seismic data for months or even years, accumulating hundreds to thousands of usable events. By spatially superimposing the PS imaging results corresponding to these events, the interface responses under different source illumination conditions can complement and enhance each other, effectively improving the imaging signal-to-noise ratio. Interface structures that are not obvious in a single event image gradually become clear after multiple superpositions. During the superposition process, the interface responses of each event's imaging results are enhanced in phase at the actual interface location, while random background noise and incoherent imaging artifacts are relatively weakened due to phase inconsistency, thus effectively improving the imaging signal-to-noise ratio.

[0092] Step S502: Superimpose the SP imaging results obtained from the second imaging branch for each seismic source, event, or time period to obtain a superimposed SP image.

[0093] The SP imaging branch records the response of incident shear waves as they are converted into P-waves at the interface, and is physically independent of the PS imaging branch. In natural earthquake observation, the incident wavefield of teleseismic events typically propagates upwards from below at a large incident angle when it reaches the receiving array, resulting in a rich shear wave component. This gives the SP imaging branch a unique advantage in certain tectonic environments. For example, in the detection of the Moho discontinuity or discontinuities within the lithosphere, teleseismic S-waves and their converted P-waves at large epicentral distances can illuminate the discontinuity at steep angles, and SP superimposed images can reveal interface details that are not easily visible in PS images. Similar to step S501, superimposing SP imaging results from multiple events can effectively enhance the signal-to-noise ratio of the SP conversion response and the structural complementarity under different event illumination conditions. Furthermore, by setting appropriate event selection criteria, such as signal-to-noise ratio thresholds, epicentral distance ranges, or focal depth ranges, quality control of the events participating in the superposition can be performed, further improving the reliability of the final superimposed image.

[0094] Step S503: Obtain the underground structure image of the target area based on the PS overlay image and the SP overlay image.

[0095] PS and SP superimposed images reflect the response characteristics of the same subsurface structure under different wavefield conversion mechanisms from two different perspectives: P-wave incident conversion and S-wave incident conversion. In certain tectonic environments, the response intensity of the two images to the same interface may differ. For example, in fluid-bearing fault zones, the converted waves generated by P-wave incident and S-wave incident have different sensitivities to fluid saturation. Comparing the two can reveal interface properties related to fluid distribution. In practical applications, PS and SP superimposed images can be interpreted as independent imaging profiles, or they can be juxtaposed and superimposed. By comprehensively identifying the interface locations jointly indicated by both images, the spatial distribution of subsurface faults, fissures, lithological interfaces, or lithospheric discontinuities in the target area can be inferred, thus obtaining a more comprehensive understanding of the subsurface structure than a single imaging type.

[0096] By superimposing and jointly interpreting the results of multi-source, multi-event PS and SP imaging respectively, the signal-to-noise ratio and structural continuity of passive source converted wave imaging are effectively improved, enabling the final underground structure image to more comprehensively and reliably reflect the spatial location, distribution morphology and physical property differences of each interface within the target area.

[0097] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor and a memory connected to the at least one processor. The memory stores instructions executable by the at least one processor, which, when executed by the at least one processor, cause the at least one processor to perform the aforementioned passive seismic reverse-time migration imaging method based on vector wavefield decomposition.

[0098] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described passive seismic reverse time migration imaging method based on vector wave field decomposition.

[0099] The embodiments of this invention aim to protect a passive seismic reverse-time migration imaging method based on vector wavefield decomposition, which has the following effects: 1. By separating the P-wave and S-wave of the recording layer vector and pre-screening the direct and converted waves based on travel time information, the coupled seismic records are decomposed into four types of target wavefields with clear physical types and corresponding to the imaging branches before reverse time extrapolation. This cuts off the path of non-physical mode conversion of non-target wavefields during reverse time extrapolation, effectively eliminating wavefield crosstalk and source position interference introduced by direct backpropagation of coupled records in traditional methods, and significantly improving the clarity and reliability of passive seismic reverse time migration imaging results. 2. The solution involves separating the P-waves and S-waves again after the reverse-time extrapolation of the recorded separation and extracting the pure wave field components. This ensures that the wave fields participating in the imaging are strictly corresponding in terms of physical type, eliminating the secondary crosstalk artifacts introduced by the residual coupled wave field in the imaging stage. At the same time, it avoids the polarity reversal and amplitude distortion problems caused by divergence-curl operations in traditional Helmholtz decomposition, ensuring the polarity correctness and physical interpretability of the imaging results. 3. By adopting amplitude correction normalization imaging conditions that preserve symbol information, reconstruct amplitude by multiplying by modulus and perform illumination normalization, the imaging amplitude attenuation and geometric diffusion caused by the angle between the vibration directions of P-wave and S-wave are effectively compensated while maintaining the interface polarity information. This overcomes the problem of insufficient amplitude fidelity in passive source converted wave imaging under traditional dot-matrix imaging conditions, making the imaging results more accurately reflect the true conversion response intensity of the subsurface interface.

[0100] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0101] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0102] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0103] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0104] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A passive seismic reverse-time migration imaging method based on vector wavefield decomposition, characterized in that, Includes the following steps: Step S1: Obtain multi-component passive seismic records and migration velocity models for the target area. The multi-component passive seismic records contain at least two orthogonal particle velocity components. Step S2: Separate the P-wave and S-wave of the multi-component passive seismic record while preserving vector attributes to obtain P-wave and S-wave records; Step S3: Calculate travel time information based on the offset velocity model, and extract direct P-wave and SP converted wave from the P-wave record and direct S-wave and PS converted wave from the S-wave record according to the travel time information; Step S4: Using the direct P-wave record, SP converted wave record, direct S-wave record, and PS converted wave record as boundary conditions, perform reverse time extrapolation of the elastic wave field to obtain the corresponding direct P-wave field, SP converted wave field, direct S-wave field, and PS converted wave field. Step S5: Perform PS imaging using the direct P-wave field and the PS converted wave field, and perform SP imaging using the SP converted wave field and the direct S-wave field to obtain an image of the underground structure of the target area.

2. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to claim 1, characterized in that, Step S4, after using the direct P-wave record, SP converted wave record, direct S-wave record, and PS converted wave record as boundary conditions for elastic wave field reverse-time extrapolation, further includes: Step S41: Separate the P-wave and S-wave in the direct P-wave field, and extract the pure direct P-wave field component from the separation result; Step S42: Separate the PS converted wave field from the P-wave to the S-wave, and extract the pure PS converted wave field component from the separation result; Step S43: Use the pure direct P-wave wave field component and the pure PS converted wave field component as inputs for PS imaging in step S5.

3. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to claim 2, characterized in that, After separating the P-wave and S-wave of the PS converted wave field in step S42 and extracting the pure PS converted wave field component from the separation result, the method further includes: Step S421: Perform P-wave and S-wave separation on the direct S-wave field, and extract the pure direct S-wave field component from the separation result; Step S422: Separate the SP converted wave field from the P-wave to the S-wave, and extract the pure SP converted wave field component from the separation result; Step S423: Use the pure SP converted wave field component and the pure direct S-wave field component as inputs for SP imaging in step S5.

4. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to claim 3, characterized in that, In step S5, PS imaging and SP imaging are performed, specifically using amplitude correction normalization imaging processing. The amplitude correction and normalization imaging processing includes: Step S51: Retain the sign information of the dot product result of the two wavefield components involved in imaging; Step S52: Reconstruct the amplitude information using the product of the magnitudes of the two wave field components; Step S53: Normalize the illumination terms of the corresponding direct wave wave field components in the two wave field components. Step S54: Based on the symbol information, the reconstructed amplitude information, and the normalization processing result, generate the imaging result.

5. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to claim 4, characterized in that, The illumination term corresponding to the direct wave field component in step S53 is specifically the autocorrelation of the corresponding direct wave field component.

6. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to claim 4, characterized in that, The reconstructed amplitude information in step S54 includes: Step S541: Based on the geometric relationship that the vibration directions of P-wave and S-wave are nearly orthogonal in an isotropic elastic medium, the vector inner product of the two wave field components is separated into a sign term and an amplitude term. Step S542: Extract the amplitude term and use the amplitude term as the amplitude information of the amplitude correction normalization imaging process.

7. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to claim 1, characterized in that, Step S3 involves extracting the direct P-wave, SP-converted wave, direct S-wave, and PS-converted wave based on the travel time information, including: Step S311: Based on the offset velocity model, calculate the theoretical arrival times of the direct P-wave, SP-converted wave, direct S-wave, and PS-converted wave respectively; Step S312: Based on the theoretical time, set corresponding extraction time windows for the P-wave record and the S-wave record respectively; Step S313: Within the extraction window, extract the direct P wave and the SP converted wave from the P wave record, and extract the direct S wave and the PS converted wave from the S wave record.

8. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to claim 7, characterized in that, Step S3 involves extracting the PS converted wave from the S-wave record, including: Step S321: In the S-wave record, a cutoff window is set for the first preset time period after the direct S-wave, the first preset time period corresponding to surface waves or other non-target strong interference wave trains; Step S322: Extract the PS converted wave during the wake period outside the cutoff window.

9. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to any one of claims 1-8, characterized in that, Step S4 involves using the direct P-wave record, SP converted wave record, direct S-wave record, and PS converted wave record as boundary conditions for reverse-time extrapolation of the elastic wave field, including: Step S401: Using the direct P-wave record and the PS converted wave record as the first imaging branch, perform elastic wave field reverse time extrapolation to obtain the direct P-wave wave field and the PS converted wave field; Step S402: Using the direct S-wave record and the SP converted wave record as the second imaging branch, perform elastic wave field reverse time extrapolation to obtain the direct S-wave wave field and the SP converted wave field; The first imaging branch is used for subsequent PS imaging, and the second imaging branch is used for subsequent SP imaging. The reverse time extrapolation of the two branches is performed independently.

10. The passive seismic reverse-time migration imaging method based on vector wavefield decomposition according to claim 9, characterized in that, Step S5 obtains an image of the underground structure of the target area, including: Step S501: For each seismic source, each event, or each time period, the PS imaging results obtained from the first imaging branch are superimposed to obtain a PS superimposed image; Step S502: For each seismic source, each event, or each time period, the SP imaging results obtained from the second imaging branch are superimposed to obtain an SP superimposed image; Step S503: Based on the PS overlay image and the SP overlay image, obtain the underground structure image of the target area.

Images