CCP imaging method, device and equipment based on receiving function and storage medium

Abstract

The application discloses a CCP imaging method and device based on a receiving function, equipment and a storage medium, and the method comprises the following steps: preprocessing seismic three-component data of a sedimentary basin area; generating a receiving function according to the seismic three-component data and filtering the receiving function; correcting the relative time of the receiving function by using a sedimentary layer bottom interface PdPpdS seismic phase; mapping the corresponding receiving function PpPs seismic phase amplitude to a virtual interface given a layered medium velocity model and the virtual interface; moving an imaging window on each virtual interface, superimposing the amplitude values in the imaging window as the amplitude at the center point of the imaging window, and taking the center point of the imaging window as a final imaging point; and calculating the imaging point to complete CCP migration imaging. The application eliminates multiple wave reverberation in a sedimentary basin, corrects the relative time of the receiving function, and can be used in a CCP migration imaging method of a normal crust model to obtain accurate CCP migration imaging.

Description

Technical Field

[0001] This invention relates to the field of crustal structure analysis, and in particular to CCP imaging methods, apparatus, devices, and storage media based on receiver functions. Background Technology

[0002] Receiver function (RF) imaging is a commonly used method for probing subsurface discontinuities and S-wave velocity structures. Common RRF imaging methods include H-κ stacking, common conversion point (CCP) migration imaging, and RRF inversion of one-dimensional S-wave velocity structures. Among these, CCP migration imaging can more intuitively characterize the depth variations of subsurface discontinuities, providing important evidence for regional geodynamics and tectonic genesis.

[0003] The inventors discovered through research that existing technologies using receiver functions for CCP imaging have at least the following drawbacks:

[0004] For receiver functions in sedimentary basins, the IASP91 reference model used in imaging does not take into account the low-velocity sedimentary layer on the surface. Therefore, the time difference between P-wave and Ps phases corresponds to a larger crustal thickness, resulting in inaccurate imaging results.

[0005] The information disclosed in this background section is intended only to enhance the understanding of the overall background of the invention and should not be construed as an admission or in any way implying that the information constitutes prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to improve the accuracy of CCP imaging in sedimentary basin areas.

[0007] This invention provides a CCP imaging method based on a receiver function, comprising the following steps:

[0008] S11. Based on seismic observation data, preprocess the three-component seismic data of sedimentary basin areas;

[0009] S12. Generate a receiving function based on the three-component seismic data, and filter the receiving function using a time-domain predictive deconvolution or a frequency-domain resonance filter to eliminate multiple reverberation in the sedimentary layer.

[0010] S13. Correcting the relative time of the receiver function using the PdPpdS seismic phase at the bottom interface of the sedimentary layer, including: when performing CCP migration imaging of the crustal structure using the PpPs seismic phase, correcting the relative 0 time in the receiver function from the P-wave seismic phase to the PdPpdS seismic phase.

[0011] S14. Given a layered medium velocity model, set the required virtual interface, and map the corresponding receiver function amplitude to the virtual interface using the time difference formula between P-wave and PpPs-wave in the normal crust model.

[0012] S15. Move the imaging window on each of the virtual interfaces, superimpose the amplitude values ​​within the imaging window, and use the amplitude at the center point of the imaging window as the final imaging point.

[0013] S16. Calculate all imaging points in the profile in sequence, use the two-dimensional profile to display the crustal structure below the sedimentary layer, and complete CCP migration imaging.

[0014] Preferably, in this invention, the preprocessing of the three-component seismic data of the sedimentary basin region based on seismic observation data includes:

[0015] Create an earthquake catalog within the time range of earthquake observation data; for P-wave receiver functions, select earthquake events with epicentral distances between 30° and 90° and magnitudes greater than 5.5.

[0016] The theoretical arrival time of P-waves is calculated using a global standard one-dimensional velocity model; seismic event data is extracted based on the theoretical arrival time; and after removing instrument response, the seismic event data is processed to remove spikes, mean, and tilt to generate the three-component seismic data.

[0017] Preferably, in this invention, generating the receiving function based on the three-component seismic data includes:

[0018] Rotating the horizontal component of the preprocessed three-component seismic data to the radial and tangential directions includes: the horizontal component refers to the amplitude of particle vibration in the north and east directions horizontally; the radial component refers to the direction of the great circle path of the epicenter and the observation station; and the tangential direction refers to the direction orthogonal to the radial direction horizontally. The formula for rotating the horizontal component of the seismic observation data to the radial-tangential coordinate system includes:

[0019]

[0020] Where R, T, E, and N represent the radial, tangential, eastward, and northward components, respectively, and Baz represents the reverse azimuth angle;

[0021] The receiver function is generated using a frequency-domain deconvolution method, comprising: the receiver function being the response of the subsurface medium beneath the observation station; the receiver function being based on the equivalent source assumption, equating the vertical component of the earthquake to an impulse function; the receiver function including a radial receiver function and a tangential receiver function; wherein, the radial receiver function is the deconvolution of the vertical component using the rotated radial component; the tangential receiver function is the deconvolution of the tangential component with respect to the vertical component; the frequency-domain calculation formula for the receiver function includes:

[0022]

[0023]

[0024] In the formula, E R (ω) and E T (ω) represents the spectrum of the radial and tangential receiver functions; D R (ω), D T (ω) and D V (ω) represents the radial, tangential, and vertical seismic event spectra, respectively; I(ω) and S(ω) represent the instrument response spectrum and the source function spectrum, respectively.

[0025] After obtaining the radial and tangential receiver function spectra, the spectrum is transformed to the time domain by inverse Fourier transform to generate the corresponding receiver function.

[0026] Preferably, in this invention, filtering the receiving function using time-domain predictive deconvolution includes:

[0027] The multiple reverberation of the deposition layer refers to the multiple waves formed by multiple reflections and propagation within the low-velocity deposition layer.

[0028] Time-domain predictive deconvolution utilizes the periodicity of the reverberation of the receiver function to eliminate periodic multiple waves while retaining the primary wave. The calculation formula includes:

[0029]

[0030] Where, r gg (t) represents the sequence of autocorrelation functions of the receiving function, m is the prediction filter factor length, a is the multiple wave period, and c is the filter factor for constructing the deconvolution filter; its coefficient matrix is ​​a Torbritz matrix, and the equation of the Torbritz matrix is ​​solved quickly using the Levinson recursive algorithm; pre-whitening processing is required in the process of solving this system of equations, and (1+b)r is used on the main diagonal of the Torbritz matrix. gg (0) replaces r gg (0), where b is the white noise coefficient, and its value range includes 0.01-0.1;

[0031]

[0032] Convolve the receiver function with multiple reverberation into the inverse filter factor d(t) to obtain the receiver function after removing multiple reverberation.

[0033] Preferably, in this invention, filtering the receiving function using a frequency domain resonant filter includes:

[0034] The reverberant receiving function is transformed to the frequency domain via Fourier transform using a frequency domain resonant filter. In the frequency domain, the series form of the multiple waves is converted to an exponential form with decay. Then, the multiple wave reverberation is eliminated in the frequency domain by division. The calculation formula includes:

[0035] F(iω)=H(iω)(1+r0e -iωΔt )

[0036] In the formula, F(iω) is the spectrum of the receiver function after multiple reverberation removal, H(iω) is the spectrum of the receiver function before multiple reverberation removal, and 1+r0e -iωΔt For a resonant filter, the parameter r0 in the resonant filter is the reflection intensity, which is defined as the ratio of the amplitude of the first trough to the amplitude of the first peak in the receiving function; the parameter Δt is the time difference between the first peak and the first trough; the parameter r0 and the parameter Δt are determined by the normalized autocorrelation function of the receiving function, and each receiving function corresponds to a set of (r0, Δt) values.

[0037] Preferably, in this invention, mapping the corresponding receiver function amplitude onto the virtual interface using the time difference formula between P-waves and PpPs-waves from a normal crustal model includes:

[0038] The formulas for calculating the time difference between PpPs waves and direct P waves in a normal crustal model include:

[0039]

[0040] In the formula, the P-wave refers to the direct P-wave from a distant earthquake with an epicentral distance between 30° and 90°, and the PpPs converted wave refers to the converted S-wave generated by the distant P-wave at the Moho surface. PpPs V represents the time difference between the P-wave and the PpPs converted wave. S and V P This refers to the S-wave and P-wave velocities, where h is the layer thickness and p is the ray parameter.

[0041] Preferably, in this invention, the given layered medium velocity model and the setting of the required virtual interface include:

[0042] The layered medium velocity model references the P-wave and S-wave velocities of the IASP91 one-dimensional reference velocity model; the imaging depth range of the virtual interface is set, including: the maximum depth of the virtual interface for the migration imaging study of the Moho surface is set to 100km, and a virtual discontinuity is set at every 0.5km depth interval, and the P-wave and S-wave velocity values ​​at each virtual discontinuity are obtained by interpolation from the IASP91 model.

[0043] Preferably, in this invention, the given layered medium velocity model includes:

[0044] For a horizontally layered medium in spherical coordinates, the formulas for calculating the ray incident angle and the central angle include:

[0045]

[0046]

[0047] In the formula, α k Let θ be the incident angle corresponding to the ray. k R0 is the central angle corresponding to the intralayer ray, where k represents the P-wave or S-wave phase, p is the ray parameter, R0 is the Earth's radius, and R0 is the radius of the Earth. j Let be the distance from the lower interface of the j-th layer to the center of the sphere;

[0048] For PpPs seismic phases, the point of penetration from the mantle through the Mohorovičić discontinuity can be represented as:

[0049]

[0050]

[0051]

[0052] In the formula λ represents latitude and longitude respectively, and Baz is the inverse azimuth. λ0 and λ0 represent the latitude and longitude of the observation station, respectively.

[0053] In another aspect of the invention, a CCP imaging apparatus based on a receiver function is also provided, comprising:

[0054] The three-component data generation unit is used to preprocess the three-component seismic data of sedimentary basin areas based on seismic observation data.

[0055] The reverberation filtering unit is used to generate a receiving function based on the three-component seismic data, and to filter the receiving function through a time-domain predictive deconvolution or a frequency-domain resonance filter to eliminate the reverberation of multiple waves in the sedimentary layer.

[0056] The relative time correction unit is used to correct the relative time of the receiver function using the PdPpdS seismic phase at the bottom interface of the sedimentary layer, including: when performing CCP migration imaging with PpPs seismic phase relative to the crustal structure, correcting the relative 0 time in the receiver function from the P-wave seismic phase to the PdPpdS seismic phase.

[0057] The function amplitude mapping unit, given a layered medium velocity model, sets the required virtual interface, and maps the corresponding received function amplitude to the virtual interface using the time difference formula between P-wave and PpPs-wave in the normal crust model.

[0058] An imaging point determination unit is used to move the imaging window on each of the virtual interfaces, superimpose the amplitude values ​​within the imaging window as the amplitude at the center point of the imaging window, and use the center point of the imaging window as the final imaging point.

[0059] The imaging point calculation unit is used to calculate all imaging points in the profile in sequence, and to display the crustal structure below the sedimentary layer using a two-dimensional profile, thus completing CCP migration imaging.

[0060] In another aspect of this invention, a CCP imaging device based on a receiving function is also provided, comprising:

[0061] Memory, used to store computer programs;

[0062] A processor for invoking and executing the computer program to implement the various steps of the CCP imaging method based on the receiver function as described in any of the preceding claims.

[0063] In another aspect of the present invention, a storage medium is also provided, on which a computer program is stored, which, when executed by a processor, implements the various steps of the CCP imaging method based on the receiving function as described in any of the preceding claims.

[0064] The CCP imaging device based on the receiving function includes a computer program stored on a medium. The computer program includes program instructions that, when executed by a computer, cause the computer to perform the methods described in the above aspects and achieve the same technical effects.

[0065] Compared with the prior art, the present invention has the following beneficial effects:

[0066] The inventors discovered that existing CCP migration imaging methods are inaccurate when applied to receiver functions in sedimentary basins because they do not consider the low-velocity sedimentary layers at the surface. Therefore, this invention, after generating the receiver function using frequency-domain deconvolution, employs time-domain predictive deconvolution or a frequency-domain resonance filter to filter the receiver function, thereby eliminating multiple reverberation from the sedimentary layers. Then, the relative time of the filtered receiver function is corrected using the PdPpdS phase correction of the sedimentary layer's bottom interface. Thus, after mapping the corresponding receiver function amplitude onto a virtual interface using the time difference formula between P-waves and PpPs waves from a normal crustal model, CCP migration imaging is completed by calculating the amplitudes of all imaging points within the profile.

[0067] Because this invention eliminates multiple reverberation in sedimentary basins and utilizes the relative time of the PdPpdS phase correction receiver function at the lower interface of the sedimentary layer, it can be used to obtain accurate CCP migration imaging using a normal crustal model CCP migration imaging method.

[0068] The above description is merely an overview of the technical solution of the present invention. In order to better understand the technical means of the present invention and to implement it according to the contents of the specification, and to make the above and other objects, technical features and advantages of the present invention easier to understand, one or more preferred embodiments are listed below and described in detail with reference to the accompanying drawings. Attached Figure Description

[0069] To more clearly illustrate the technical solution of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0070] Figure 1 This is a flowchart illustrating the steps of the CCP imaging method based on the receiver function described in this invention;

[0071] Figures 2 to 5 This is a schematic diagram of various ray paths of the seismic phase described in this invention;

[0072] Figure 6 This is a simulation diagram of the receiver function with the reverberation effect of the deposition layer as described in this invention;

[0073] Figure 7 and Figure 8 A schematic diagram of the result of eliminating multiple reverberation of the deposition layer using time-domain deconvolution in the receiving function described in this invention;

[0074] Figure 9This is a schematic diagram of the structure of the CCP imaging device based on the receiver function described in this invention.

[0075] Figure 10 This is a schematic diagram of the CCP imaging device based on the receiving function described in this invention. Detailed Implementation

[0076] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings, but it should be understood that the scope of protection of the present invention is not limited to the specific embodiments.

[0077] Unless otherwise expressly stated, throughout the specification and claims, the term "comprising" or its variations such as "including" or "comprises" shall be understood to include the stated elements or components without excluding other elements or other components.

[0078] In this document, the terms "first," "second," etc., are used to distinguish two different elements or parts, and are not used to define specific positions or relative relationships. In other words, in some embodiments, the terms "first," "second," etc., can also be used interchangeably.

[0079] Example 1

[0080] To improve the accuracy of CCP imaging in sedimentary basin areas, such as Figure 1 As shown, this embodiment of the invention provides a CCP imaging method based on a receiver function, including the following steps:

[0081] S11 preprocesses the three-component seismic data of sedimentary basin areas based on seismic observation data;

[0082] CCP is an abbreviation for Cmmom Conversion Point, which refers to wavefield backpropagation based on the kinematic properties of the ray path between underground phase conversion points and surface receiving stations. The underground conversion points are arranged in two or three dimensions, and the amplitude of the receiving function corresponding to the time at the location of the surface observation station is backpropagated to the underground conversion point, and imaging is performed in the depth domain.

[0083] In this embodiment of the invention, the preprocessed seismic three-component data of the sedimentary basin region refers to: firstly, creating a seismic catalog within the time range of seismic observation data; for P-wave receiver functions, typically selecting seismic events with epicentral distances between 30 and 90° and magnitudes greater than 5.5; then calculating the theoretical arrival time of P-waves using global standard one-dimensional velocity models (IASP91, AK135 model, etc.); extracting seismic event data based on the theoretical arrival time; removing instrument responses; and then performing despiking, mean removal, and tilt removal processing on the seismic event data.

[0084] S12. Generate a receiving function based on the three-component seismic data, and filter the receiving function using a time-domain predictive deconvolution or a frequency-domain resonance filter to eliminate multiple reverberation in the sedimentary layer.

[0085] In this step, the horizontal component of the preprocessed seismic three-component data is rotated to the radial and tangential directions, which may specifically include:

[0086] The horizontal component refers to the amplitude of particle vibration in the north and east directions horizontally, the radial component refers to the direction of the great circle path of the epicenter and the observation station, and the tangential component refers to the direction orthogonal to the radial direction horizontally. The formulas for rotating the horizontal component of seismic observation data to the radial-tangential coordinate system include:

[0087]

[0088] Where R, T, E, and N represent the radial, tangential, eastward, and northward components, respectively, and Baz represents the reverse azimuth angle.

[0089] Next, the receiver function is generated using the frequency domain deconvolution method, which may specifically include:

[0090] The receiver function refers to the response of the subsurface medium beneath the observation station. Based on the equivalent source assumption, the receiver function equates the vertical component of the earthquake to an impulse function. The receiver function includes radial and tangential receiver functions. The radial receiver function is obtained by deconvolving the vertical component with the rotated radial component; the tangential receiver function is obtained by deconvolving the tangential component with respect to the vertical component. The formula for calculating the receiver function in the frequency domain includes:

[0091]

[0092]

[0093] In the formula, E R (ω) and E T (ω) represents the spectrum of the radial and tangential receiver functions; D R (ω), D T (ω) and D V (ω) represents the radial, tangential, and vertical seismic event spectra, respectively; I(ω) and S(ω) represent the instrument response spectrum and the source function spectrum, respectively.

[0094] After obtaining the radial and tangential receiver function spectra, the corresponding receiver function can be generated by performing an inverse Fourier transform to the time domain.

[0095] In this embodiment of the invention, multiple reverberation of the deposition layer refers to multiple waves formed by multiple reflections and propagation within the low-velocity deposition layer;

[0096] In this embodiment of the invention, filtering the receiver function to eliminate multiple reverberation of the deposition layer can be done in two ways: one is to filter the receiver function through time-domain predictive deconvolution; the other is to filter the receiver function through a frequency-domain resonant filter.

[0097] Filtering the receiver function using time-domain predictive deconvolution can specifically include: Time-domain predictive deconvolution utilizes the periodicity of the receiver function's reverberation to eliminate periodic multiple waves while retaining the primary wave. Specific calculation formulas include:

[0098]

[0099] Where, r gg (t) represents the sequence of autocorrelation functions of the receiving function, m is the prediction filter factor length, a is the multiple period, and c is the filter factor used to construct the deconvolution filter; its coefficient matrix is ​​a Tobrizoid matrix, and the equation of the Tobrizoid matrix is ​​solved quickly using the Levinson recursive algorithm; it should be noted that pre-whitening processing is required during the solution of this system of equations, and (1+b)r is used on the main diagonal of the Tobrizoid matrix. gg (0) replaces r gg (0), where b is the white noise coefficient. In this embodiment of the invention, the white noise coefficient is a very small positive number, and its value range preferably includes 0.01-0.1.

[0100]

[0101] By convolving the receiver function with multiple reverberation with the inverse filter factor d(t), the receiver function after removing multiple reverberation can be obtained.

[0102] In addition, the received function is filtered by a frequency domain resonant filter, specifically including:

[0103] The reverberant receiver function is transformed to the frequency domain via Fourier transform using a frequency domain resonant filter. In the frequency domain, the series form of the multiple waves is converted to an exponential form with decay. Then, the multiple reverberation is eliminated in the frequency domain by division. Specific calculation formulas may include:

[0104] F(iω)=H(iω)(1+r0e -iωΔt )

[0105] In the formula, F(iω) is the spectrum of the receiver function after multiple reverberation removal, H(iω) is the spectrum of the receiver function before multiple reverberation removal, and 1+r0e -iωΔtFor a resonant filter, the parameter r0 in the resonant filter is the reflection intensity, which is defined as the ratio of the amplitude of the first trough to the amplitude of the first peak in the receiving function; the parameter Δt is the time difference between the first peak and the first trough; the parameter r0 and the parameter Δt are determined by the normalized autocorrelation function of the receiving function, and each receiving function corresponds to a set of (r0, Δt) values.

[0106] S13. Correcting the relative time of the receiver function using the PdPpdS seismic phase at the bottom interface of the sedimentary layer, including: when performing CCP migration imaging of the crustal structure using the PpPs seismic phase, correcting the relative 0 time in the receiver function from the P-wave seismic phase to the PdPpdS seismic phase.

[0107] S14. Given a layered medium velocity model, set the required virtual interface, and map the corresponding receiver function amplitude to the virtual interface using the time difference formula between P-wave and PpPs-wave in the normal crust model.

[0108] In this step, the formula for calculating the time difference between the PpPs wave and the teleseismic direct P wave in the normal crustal model may specifically include:

[0109]

[0110] In the formula, the P-wave refers to the direct P-wave from a distant earthquake with an epicentral distance between 30° and 90°, and the PpPs converted wave refers to the converted S-wave generated by the distant P-wave at the Moho surface. PpPs V represents the time difference between the P-wave and the PpPs converted wave. S and V P This refers to the S-wave and P-wave velocities, where h is the layer thickness and p is the ray parameter.

[0111] S15. Move the imaging window on each of the virtual interfaces, superimpose the amplitude values ​​within the imaging window, and use the amplitude at the center point of the imaging window as the final imaging point.

[0112] In this step, the layered medium velocity model can refer to the P-wave and S-wave velocities of the IASP91 one-dimensional reference velocity model; the imaging depth range of the virtual interface (i.e., the underground virtual interface) is set, including: the maximum depth of the virtual interface for the migration imaging study of the Moho surface is set to 100km, and a virtual discontinuity is set at every 0.5km depth interval, and the P-wave and S-wave velocity values ​​at each virtual discontinuity are obtained by interpolation from the IASP91 model.

[0113] In practical applications, the layered medium velocity model in this embodiment of the invention may specifically include:

[0114] For a horizontally layered medium in spherical coordinates, the formulas for calculating the ray incident angle and the central angle include:

[0115]

[0116]

[0117] In the formula, α k Let θ be the incident angle corresponding to the ray. k R0 is the central angle corresponding to the intralayer ray, where k represents the P-wave or S-wave phase, p is the ray parameter, R0 is the Earth's radius, and R0 is the radius of the Earth. j Let be the distance from the lower interface of the j-th layer to the center of the sphere;

[0118] For PpPs seismic phases, the point of penetration from the mantle through the Mohorovičić discontinuity can be represented as:

[0119]

[0120]

[0121]

[0122] In the formula λ represents latitude and longitude respectively, and Baz is the inverse azimuth. λ0 and λ0 represent the latitude and longitude of the observation station, respectively.

[0123] S16. Calculate all imaging points in the profile in sequence, use the two-dimensional profile to display the crustal structure below the sedimentary layer, and complete CCP migration imaging.

[0124] The following specific example illustrates the working principle and technical effects of the embodiments of the present invention;

[0125] Figures 2 to 5 These are schematic diagrams of various ray paths of the seismic phases in the above example, used to illustrate that the receiver function of a sedimentary basin can be corrected relative to time 0 and normal PpP. S The time difference formula is used to calculate the crustal structure beneath the sedimentary layer. Among other things, Figure 2 The path of P-wave rays passing through the Earth's crust and sedimentary layers; Figure 3 The ray path for the PdPpdS phase; Figure 4 The ray path for multiple reverberation of the sedimentary layer. Figure 5 The ray path for the PpPs phase.

[0126] In the above figure, solid lines represent P-wave paths, dashed lines represent S-wave paths, black lines represent propagation paths within the Earth's crust, red lines represent propagation paths within sedimentary layers, and orange lines represent the superposition of all seismic phases. Figure 2 and Figure 3The time difference represents the time difference between the direct P-wave and the PdPpdS-wave, that is, the correction value used in this embodiment of the invention to correct the relative time of the receiving function. Figure 2 and Figure 4 The ray path represents the time difference of the converted wave within the crust, which is the PpPs seismic back-propagation time difference required for imaging.

[0127] Then, by means of Figures 6 to 8 The illustration demonstrates the effectiveness of the correction method proposed in this embodiment of the invention by utilizing a simulated reception function. Figure 6 The receiver function is used to simulate the reverberation effect of the sediment layer. The vertical axis represents the sediment layer thickness, varying from 0 km to 1 km in 0.1 km intervals. The P-wave velocity of the sediment layer is 2.1 km / s, the S-wave velocity is 0.7 km / s, and the density is 2500 kg / m³. 3 The Earth's crust has a P-wave velocity of 6.1 km / s, an S-wave velocity of 3.49 km / s, and a density of 2700 kg / m³. 3 The thickness is 32 km; the P-wave velocity of the lower mantle is 8.0 km / s, the S-wave velocity is 4.5 km / s, and the density is 3300 kg / m³. 3 The deconvolution method uses frequency domain deconvolution, and the Gaussian coefficient is selected as 5.

[0128] Figure 7 for Figure 6 The receiver function utilizes time-domain deconvolution to eliminate multiple reverberation in the sediment layer. Figure 7 The five dashed lines represent the theoretical time difference lines for crustal P-waves, sedimentary bottom interface conversion waves and multiples, Ps waves, and PpPs multiples, respectively. It is clear that... Figure 7 The time difference of the PpPs wave does not correspond to the corresponding PpPs wave phase in the receiver function, therefore it cannot be used in CCP migration imaging.

[0129] Figure 8 The result is the result after correction using the method of the embodiments of the present invention. Figure 8 The first dashed line on the left represents the location of the PdPpdS phase, and the second dashed line represents the time difference between the theoretical P wave and PdPpdS calculated based on the crustal model. It can be seen that the theoretical time difference corresponds well with the PpPs wave, indicating that the method proposed in the embodiments of the invention can accurately estimate the crustal structure of sedimentary basin areas.

[0130] In summary, after calculating and generating the receiver function using the frequency domain deconvolution method, this embodiment of the invention further employs time domain predictive deconvolution or a frequency domain resonance filter to filter the receiver function, thereby eliminating multiple reverberation of the sedimentary layer. Then, the relative time of the filtered receiver function is corrected using the PdPpdS phase correction of the sedimentary layer's bottom interface. Thus, after mapping the corresponding receiver function amplitude onto the virtual interface using the time difference formula between P-waves and PpPs waves from a normal crustal model, CCP migration imaging is completed by calculating the amplitude of all imaging points within the profile.

[0131] Because the embodiments of the present invention eliminate multiple reverberation in sedimentary basins and utilize the relative time of the PdPpdS phase correction receiver function at the lower interface of the sedimentary layer, the CCP migration imaging method for normal crustal models can be used to obtain accurate CCP migration imaging.

[0132] Example 2

[0133] Corresponding to the method embodiment, another aspect of the present invention provides a CCP imaging device based on a receiver function. Figure 2 This diagram illustrates the structure of a CCP imaging device based on a reception function according to an embodiment of the present invention. The CCP imaging device based on a reception function is... Figure 1 The device corresponding to the CCP imaging method based on the receiver function described in the corresponding embodiment is implemented through a virtual device. Figure 1 In the corresponding embodiment of the CCP imaging method based on the receiver function, the various virtual modules constituting the CCP imaging device based on the receiver function can be executed by electronic devices, such as network devices, terminal devices, or servers. Specifically, the CCP imaging device based on the receiver function in this embodiment of the invention includes:

[0134] Three-component data generation unit 01 is used to preprocess the three-component seismic data of sedimentary basin areas based on seismic observation data.

[0135] The reverberation filtering unit 02 is used to generate a receiving function based on the three-component seismic data, and to filter the receiving function through a time-domain predictive deconvolution or a frequency-domain resonance filter to eliminate multiple reverberation of the sedimentary layer.

[0136] The relative time correction unit 03 is used to correct the relative time of the receiver function using the PdPpdS seismic phase at the bottom interface of the sedimentary layer, including: when performing CCP migration imaging of the crustal structure using the PpPs seismic phase, correcting the relative 0 time in the receiver function from the P-wave seismic phase to the PdPpdS seismic phase.

[0137] Function amplitude mapping unit 04, given a layered medium velocity model, sets the required virtual interface, and maps the corresponding received function amplitude to the virtual interface through the time difference formula between P wave and PpPs wave in the normal crust model.

[0138] The imaging point determination unit 05 is used to move the imaging window on each of the virtual interfaces, superimpose the amplitude values ​​within the imaging window as the amplitude at the center point of the imaging window, and use the center point of the imaging window as the final imaging point.

[0139] Imaging point calculation unit 06 is used to calculate all imaging points in the profile in sequence, and to display the crustal structure below the sedimentary layer using a two-dimensional profile, thus completing CCP migration imaging.

[0140] It should be noted that the specific implementation and technical effects of the CCP imaging device based on the receiver function in the embodiments of the present invention can be found by referring to... Figure 1 The corresponding CCP imaging methods based on the receiver function will not be elaborated here.

[0141] Example 3

[0142] Corresponding to the method embodiments, this invention also provides a CCP imaging device based on a receiving function, such as a terminal or a server. The server can be a standalone physical server, a server cluster or distributed system composed of multiple physical servers, or a cloud server providing basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communication, middleware services, domain name services, security services, CDN, and big data and artificial intelligence platforms. The terminal can be a smartphone, tablet, laptop, desktop computer, etc., but is not limited to these.

[0143] An example diagram of the hardware structure block diagram of the CCP imaging device based on the receiving function provided in this application is shown below. Figure 9 As shown, it may include:

[0144] Processor 1, communication interface 2, memory 3, and communication bus 4;

[0145] The processor 1, communication interface 2, and memory 3 communicate with each other via communication bus 4.

[0146] Optionally, communication interface 2 can be an interface of a communication module, such as the interface of a GSM module;

[0147] Processor 1 may be a central processing unit (CPU), an application-specific integrated circuit (ASIC), or one or more integrated circuits configured to implement the embodiments of this application.

[0148] Memory 3 may include high-speed RAM memory, and may also include non-volatile memory, such as at least one disk storage device.

[0149] Specifically, processor 1 is used to execute the computer program stored in memory 3 to perform the following steps:

[0150] S11. Process the three-component seismic data of sedimentary basin areas based on seismic observation data;

[0151] S12. Generate a receiving function based on the three-component seismic data, and filter the receiving function using a time-domain predictive deconvolution or a frequency-domain resonance filter to eliminate multiple reverberation in the sedimentary layer.

[0152] S13. Correcting the relative time of the receiver function using the PdPpdS seismic phase at the bottom interface of the sedimentary layer, including: when performing CCP migration imaging of the crustal structure using the PpPs seismic phase, correcting the relative 0 time in the receiver function from the P-wave seismic phase to the PdPpdS seismic phase.

[0153] S14. Given a layered medium velocity model, set the required virtual interface, and map the corresponding receiver function amplitude to the virtual interface using the time difference formula between P-wave and PpPs-wave in the normal crust model.

[0154] S15. Move the imaging window on each of the virtual interfaces, superimpose the amplitude values ​​within the imaging window, and use the amplitude at the center point of the imaging window as the final imaging point.

[0155] S16. Calculate all imaging points in the profile in sequence, use the two-dimensional profile to display the crustal structure below the sedimentary layer, and complete CCP migration imaging.

[0156] The above-described product can execute the method provided in the embodiments of the present invention, and has the corresponding functional modules and beneficial effects for executing the method. Technical details not described in detail in this embodiment can be found in the CCP imaging method based on the receiver function provided in the embodiments of the present invention.

[0157] Example 4

[0158] In this embodiment of the invention, a storage medium is also provided, which can store a program suitable for execution by a processor, the program being used for:

[0159] S11. Process the three-component seismic data of sedimentary basin areas based on seismic observation data;

[0160] S12. Generate a receiving function based on the three-component seismic data, and filter the receiving function using a time-domain predictive deconvolution or a frequency-domain resonance filter to eliminate multiple reverberation in the sedimentary layer.

[0161] S13. Correcting the relative time of the receiver function using the PdPpdS seismic phase at the bottom interface of the sedimentary layer, including: when performing CCP migration imaging of the crustal structure using the PpPs seismic phase, correcting the relative 0 time in the receiver function from the P-wave seismic phase to the PdPpdS seismic phase.

[0162] S14. Given a layered medium velocity model, set the required virtual interface, and map the corresponding receiver function amplitude to the virtual interface using the time difference formula between P-wave and PpPs-wave in the normal crust model.

[0163] S15. Move the imaging window on each of the virtual interfaces, superimpose the amplitude values ​​within the imaging window, and use the amplitude at the center point of the imaging window as the final imaging point.

[0164] S16. Calculate all imaging points in the profile in sequence, use the two-dimensional profile to display the crustal structure below the sedimentary layer, and complete CCP migration imaging.

[0165] Optionally, the refined and extended functions of the program can be found in the description above.

[0166] The above-described product can execute the methods provided in the embodiments of the present invention, and has the corresponding functional modules and beneficial effects for executing the methods. Technical details not described in detail in this embodiment can be found in the methods provided in other embodiments of the present invention.

[0167] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0168] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. Furthermore, the couplings or direct couplings or communication connections shown or discussed may be indirect couplings or communication connections through interfaces, devices, or units, and may be electrical, mechanical, or other forms.

[0169] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0170] In addition, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0171] It should be understood that in the embodiments of this application, the claims, various embodiments, and features can be combined with each other to solve the aforementioned technical problems.

[0172] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0173] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A CCP imaging method based on a receiver function, characterized in that, Including the following steps: S11. Based on seismic observation data, preprocess the three-component seismic data of sedimentary basin areas; S12. Generate a receiving function based on the three-component seismic data, and filter the receiving function using a time-domain predictive deconvolution or a frequency-domain resonance filter to eliminate multiple reverberation in the sedimentary layer. S13, Utilizing the bottom interface of the sedimentary layer Phase correction of the relative time of the received function includes: using When performing CCP migration imaging of the seismic relative crustal structure, the relative time 0 in the receiver function is corrected from P-wave seismic relative to... Zhen and his companions; S14. Given a layered medium velocity model, set the required virtual interface, and use P-waves from a normal crustal model. The time difference formula between waves maps the corresponding amplitude of the receiving function to the virtual interface; S15. Move the imaging window on each of the virtual interfaces, superimpose the amplitude values ​​within the imaging window, and use the amplitude at the center point of the imaging window as the final imaging point. S16. Calculate all imaging points in the profile in sequence, use the two-dimensional profile to display the crustal structure below the sedimentary layer, and complete CCP migration imaging.

2. The CCP imaging method based on the receiver function according to claim 1, characterized in that, The preprocessing of seismic three-component data in sedimentary basin areas based on seismic observation data includes: Create an earthquake catalog within the time frame of the earthquake observation data; for The wave receiver function is selected for earthquake events with an epicentral distance between 30° and 90° and a magnitude greater than 5.

5. Calculation using the global standard one-dimensional velocity model Theoretical arrival time of waves is determined; seismic event data is extracted based on the theoretical arrival time; and after removing instrument response, the seismic event data is processed to remove spikes, mean, and tilt to generate the three-component seismic data.

3. The CCP imaging method based on the receiver function according to claim 1, characterized in that, The step of generating a receiving function based on the three-component seismic data includes: Rotating the horizontal component of the preprocessed three-component seismic data to the radial and tangential directions includes: the horizontal component refers to the amplitude of particle vibration in the north and east directions horizontally; the radial component refers to the direction of the great circle path of the epicenter and the observation station; and the tangential direction refers to the direction orthogonal to the radial direction horizontally. The formula for rotating the horizontal component of the seismic observation data to the radial-tangential coordinate system includes: ； in, These represent the radial, tangential, eastward, and northward components, respectively. Indicates the reverse azimuth angle; The receiver function is generated using a frequency-domain deconvolution method, comprising: the receiver function being the response of the subsurface medium beneath the observation station; the receiver function being based on the equivalent source assumption, equating the vertical component of the earthquake to an impulse function; the receiver function including a radial receiver function and a tangential receiver function; wherein, the radial receiver function is the deconvolution of the vertical component using the rotated radial component; the tangential receiver function is the deconvolution of the tangential component with respect to the vertical component; the frequency-domain calculation formula for the receiver function includes: ； In the formula, and The spectrum of the radial and tangential receiver functions; , and These are the seismic event spectra in the radial, tangential, and vertical directions, respectively. and These are the instrument response spectrum and the source function spectrum, respectively. After obtaining the radial and tangential receiver function spectra, the spectrum is transformed to the time domain by inverse Fourier transform to generate the corresponding receiver function.

4. The CCP imaging method based on the receiver function according to claim 2, characterized in that, The filtering of the receiving function by predicting deconvolution in the time domain includes: The multiple reverberation of the deposition layer refers to the multiple waves formed by multiple reflections and propagation within the low-velocity deposition layer. Time-domain predictive deconvolution utilizes the periodicity of the reverberation of the receiver function to eliminate periodic multiple waves while retaining the primary wave. The calculation formula includes: ； in, For the sequence of autocorrelation functions of the receiving function, To predict the filter factor length, For multiple wave cycles, To construct the filtering factor of the deconvolution filter; its coefficient matrix is ​​a Tobrizoid matrix, the equation of which is solved quickly using the Levinson recursive algorithm; pre-whitening processing is required during the solution of this system of equations, and the main diagonal of the Tobrizoid matrix is ​​represented by... replace ,in This is the white noise coefficient, and its value ranges from 0.01 to 0.

1. ； Convolve the receiver function with multiple reverberation to the inverse filter factor Obtain the receiver function after removing multiple reverberation.

5. The CCP imaging method based on the receiver function according to claim 1, characterized in that, The filtering of the received function using a frequency domain resonant filter includes: The reverberant receiver function is transformed to the frequency domain via a Fourier transform using a frequency domain resonant filter. In the frequency domain, the series form of the multiple waves is converted into attenuated forms. The exponential form, followed by division in the frequency domain to eliminate multiple reverberation, is calculated using the following formulas: ； In the formula, To eliminate the spectrum of the receiver function after multiple reverberations, To eliminate the spectrum of the receiver function before multiple reverberations, This is a resonant filter, and the parameters in the resonant filter are... The reflection intensity is defined as the ratio of the amplitude of the first trough to the amplitude of the first peak in the receiver function; parameter The time difference between the first peak and the first trough; this parameter is determined using the normalized autocorrelation function of the receiver function. and parameters Each receiving function corresponds to a set value.

6. The CCP imaging method based on the receiver function according to claim 2, characterized in that, The passage through the normal crustal model P wave and The time difference formula between waves maps the corresponding amplitude of the receiving function onto the virtual interface, including: Normal crustal model Wave and direct P The formulas for calculating the time difference between waves include: ； In the formula, the aforementioned P Waves refer to distal earthquakes with epicentral distances between 30° and 90°. P Wave, Converted waves refer to teleseismic waves P The converted S-wave generated by the wave at the Moho surface for P wave and Time difference between converted waves and It means wave and P wave velocity, For layer thickness, These are the ray parameters.

7. The CCP imaging method based on the receiver function according to claim 1, characterized in that, Given a layered medium velocity model, the required virtual interface is set, including: The layered medium velocity model references the IASP91 one-dimensional reference velocity model. P wave and Wave velocity; setting the imaging depth range of the virtual interface, including: for the migration imaging study of the Moho surface, the maximum depth of the virtual interface is set to 100km, and a virtual discontinuity is set at every 0.5km depth interval, and at each virtual discontinuity... P wave and The wave velocity values ​​were obtained by interpolation using the IASP91 model.

8. The CCP imaging method based on the receiver function according to claim 7, characterized in that, The layered medium velocity model includes: For a horizontally layered medium in spherical coordinates, the formulas for calculating the ray incident angle and the central angle include: ； In the formula, The incident angle corresponding to the ray. The central angle corresponding to the ray within the layer is denoted as , where represent P wave or Wave phase, For ray parameters, For the Earth's radius, For the first The distance from the lower interface to the center of the sphere; for The seismic phase, whose point of penetration from the mantle through the Mohorovičić discontinuity, can be represented as: ； In the formula They represent latitude and longitude, respectively. For the opposite azimuth angle, These represent the latitude and longitude of the observation station, respectively.

9. A CCP imaging device based on a receiver function, characterized in that, include: The three-component data generation unit is used to preprocess the three-component seismic data of sedimentary basin areas based on seismic observation data. The reverberation filtering unit is used to generate a receiving function based on the three-component seismic data, and to filter the receiving function through a time-domain predictive deconvolution or a frequency-domain resonance filter to eliminate the reverberation of multiple waves in the sedimentary layer. Relative time correction unit, used to utilize the bottom interface of the sedimentary layer Phase correction of the relative time of the received function includes: using When performing CCP migration imaging of the seismic relative crustal structure, the relative time 0 in the receiver function is corrected from P-wave seismic relative to... Zhen and his companions; Function amplitude mapping element, given a layered medium velocity model, sets the required virtual interface, and passes through the P-wave and normal crust model. The time difference formula between waves maps the corresponding amplitude of the receiving function to the virtual interface; An imaging point determination unit is used to move the imaging window on each of the virtual interfaces, superimpose the amplitude values ​​within the imaging window as the amplitude at the center point of the imaging window, and use the center point of the imaging window as the final imaging point. The imaging point calculation unit is used to calculate all imaging points in the profile in sequence, and to display the crustal structure below the sedimentary layer using a two-dimensional profile, thus completing CCP migration imaging.

10. A CCP imaging device based on a receiver function, characterized in that, include: Memory, used to store computer programs; A processor for invoking and executing the computer program to implement the steps of the CCP imaging method based on the receiver function as described in any one of claims 1-8.

11. A storage medium, characterized in that, Includes a software program adapted by a processor to perform the steps of the CCP imaging method based on the receiver function as described in any one of claims 1-8.

Images