PP and PS combined reverse time migration imaging method based on cosine similarity measurement

Abstract

The invention provides a PP and PS combined reverse time migration imaging method based on cosine similarity measurement, and the method comprises the steps: 1, obtaining observed seismic data, and extracting seismic wavelets; step 2, calculating a P-wave forward propagation wave field at a seismic source end, a P-wave backward propagation wave field at a detection point end and an S-wave backward propagation wave field at the detection point end; 3, calculating a PP imaging data volume and a PS imaging data volume; 4, calculating the cosine similarity of the PP and PS imaging data volumes; 5, filtering the PP imaging data volume and the PS imaging data volume by using cosine similarity; and step 6, superposing the filtered PP imaging data volume and PS imaging data volume corresponding to all the single shot records, and outputting the PP imaging data volume and the PS imaging data volume. According to the PP and PS combined reverse time migration imaging method based on cosine similarity measurement, the overall quality of PP and PS imaging data volumes is improved, calculation is simple, implementation is easy, high practicability is achieved, and the reliability of the imaging result is high.

Description

Technical Field

[0001] This invention relates to the field of seismic data elastic wave reverse time migration imaging technology, and in particular to a joint reverse time migration imaging method of PP and PS based on cosine similarity metric. Background Technology

[0002] Elastic wave reverse time migration imaging (RTM) using nodal-recorded multi-wave, multi-component data is an important method for acquiring P-wave (PP) and converted S-wave (PS) images of subsurface media. Since S-wave velocity is lower than P-wave velocity, the resolution of the phase axis in PS imaging profiles is higher than that in PP imaging profiles. This makes PS imaging profiles more effective for identifying thin interlayers. Furthermore, S-waves are insensitive to fluids. Therefore, PS imaging profiles offer higher image quality in gas cloud areas and fractured buried hills. While existing RTM methods can simultaneously acquire PP and PS imaging profiles, significant artifacts exist in both. This leads to a decrease in the overall quality of PP and PS imaging profiles and interferes with subsequent seismic data interpretation, hindering accurate identification of subsurface geological structures. Current techniques primarily suppress noise in P-wave migration imaging but cannot suppress noise in multi-wave, multi-component data migration. Therefore, there is an urgent need to develop noise suppression techniques for multi-wave, multi-component data migration.

[0003] Chinese patent application CN112379430B discloses an angle-domain multi-component migration imaging method. This method divides the seismic records of PP and PS wave co-shot points for each shot into multiple local windows, determining the beam center point. Local windows are then locally tilted and stacked to obtain local plane waves. The propagation angles of the P-waves originating from the source, the central ray of the P-waves originating from the beam center point, and the central ray of the S-waves originating from the beam center point are determined. The local plane waves are then reverse-extended to obtain the reverse-extended wavefield, and the P-wave incident angles at the imaging points are determined for both PP and PS waves. Finally, the imaging values ​​from all preset incident angle ranges within the extracted PP and PS wave angle-domain co-shot point gathers are superimposed to obtain PP and PS wave angle-domain pre-stack depth migration imaging profiles. This invention achieves efficient migration imaging of PP and PS wave seismic data in the angle domain, obtaining high-quality multi-component migration imaging profiles. The aforementioned patent is based on ray theory and requires calculating the incident angle of the P-wave and PS-wave at the imaging point. However, it does not consider the similarity between the PP and PS imaging profiles during the offset imaging process, leading to inconsistencies between them. Furthermore, the patent does not account for noise during the offset process, resulting in offset artifacts in the PP and PS imaging profiles and reducing overall image quality.

[0004] In Chinese patent application publication number CN115755175A, a method for reverse time migration imaging of elastic waves based on the Hilbert transform is disclosed. This method involves obtaining HβHt, H-βHt, and HβH-β terms through time and oblique direction Hilbert transforms, respectively; substituting these terms into the oblique direction Hilbert full-field separation formula to obtain four wavefield components separated from the original shot-receiver wavefield data along any β and its symmetric direction; obtaining the corresponding forward propagation composite P-wave (P-wave) wavefield and reverse propagation P-wave and S-wave (S-wave) composite wavefield based on the four wavefield components; and cross-correlating the forward propagation composite wavefield with the reverse propagation up-wave and the forward propagation down-wave with the reverse propagation composite wavefield to obtain the PP and PS elastic wave reverse time migration results that retain the steep structure. This invention utilizes the construction of a composite wave field to eliminate gyratory wave artifacts and reduce low-frequency noise during migration, while retaining the algorithm's advantage in imaging steep structures. Furthermore, compared to the traditional elastic wave reverse time migration algorithm, its computational complexity only increases by about double. This provides a new approach for using elastic wave reverse time migration algorithms to image complex, steep underground structures. The aforementioned patent requires Hilbert transform; moreover, compared to the traditional elastic wave reverse time migration algorithm, the patent's computational complexity is doubled. The patent also does not consider the similarity between the PP and PS imaging profiles during migration imaging, leading to inconsistencies between the PP and PS imaging profiles.

[0005] Chinese patent application CN114814944B discloses a method for elastic P-wave and S-wave field separation and reverse-time migration imaging based on divergence and curl. The method includes: Step 1, obtaining the total vector seismic wave field U based on forward modeling of seismic waves; Step 2, solving for the Helmholtz potential P-wave Pcor and S-wave Scor after amplitude and phase correction; and solving for the vector P-wave UP and vector S-wave US; and Step 3, generating PP reflectance images and PS reflectance images using scalar imaging conditions and vector imaging conditions based on amplitude and sign, respectively. This invention solves the technical problem that existing reverse-time migration imaging overlay profiles cannot reflect true subsurface reflectance information. The aforementioned patent requires Helmholtz decomposition; furthermore, the P-waves and S-waves obtained from Helmholtz decomposition exhibit phase and amplitude distortion, resulting in amplitude-insecure PP and PS reflectances. The aforementioned patent does not consider the similarity between the PP and PS imaging profiles during migration imaging, leading to inconsistencies between the PP and PS imaging profiles.

[0006] The existing technologies described above are significantly different from this invention. A search reveals no literature in the XY category, indicating that this invention possesses novelty. Since no solution exists in the existing technology to address the technical problem we seek, we have invented a novel PP and PS joint reverse time migration imaging method based on cosine similarity measurement. Summary of the Invention

[0007] The purpose of this invention is to provide a joint reverse time migration imaging method based on cosine similarity metric, which is computationally simple, easy to implement, highly practical, and produces reliable imaging results.

[0008] The objective of this invention can be achieved through the following technical measures: a joint reverse time migration imaging method for PP and PS based on cosine similarity measurement, comprising:

[0009] Step 1: Acquire observed seismic data and extract seismic wavelets;

[0010] Step 2: Calculate the forward propagation wavefield of the P-wave at the source end, the reverse propagation wavefield of the P-wave at the receiver end, and the reverse propagation wavefield of the S-wave at the receiver end.

[0011] Step 3: Calculate the PP imaging data volume and the PS imaging data volume;

[0012] Step 4: Calculate the cosine similarity between the PP and PS imaging data volumes;

[0013] Step 5: Filter the PP imaging data volume and the PS imaging data volume using cosine similarity;

[0014] Step 6: Superimpose the filtered PP imaging data volume and PS imaging data volume corresponding to all single-shot records, and output the PP and PS imaging data volumes.

[0015] The objective of this invention can also be achieved through the following technical measures:

[0016] In step 1, the observed seismic data is acquired by placing a four-component nodal seismograph on the seabed and recording the artificially generated source wave field. The recorded data is used as the observed seismic data. For the theoretical model, the observed seismic data is obtained by forward modeling the elastic wave equation based on the true velocity model.

[0017] In step 1, after acquiring the observed seismic data, the observed seismic data is preprocessed to obtain preprocessed seismic data. The preprocessing includes noise suppression, multiple wave suppression, node attitude correction, node position correction, and removal of direct and refracted waves.

[0018] In step 1, the seismic wavelet is extracted using the preprocessed seismic data and a higher-order statistical method.

[0019] In step 2, using the given migration velocity model and the extracted seismic wavelet, the forward propagation wavefield of the P-wave at the source end is calculated. The calculation formula is as follows:

[0020]

[0021] Where ρ represents the density of the medium, k represents the bulk modulus of the medium, μ represents the shear modulus of the medium, and v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x s v represents the transverse wave velocity component of a particle vibrating in the x-direction. z s v represents the transverse wave velocity component of a particle vibrating in the z-direction. x v represents the velocity component of a particle vibrating in the x-direction. z Let τ represent the velocity component of the particle vibrating in the z-direction, p represent the pressure field, and τ represent the velocity component of the particle vibrating in the z-direction. xx τ represents the normal stress field in the x-direction. zz τ represents the normal stress field in the z-direction. xz f represents the shear stress field in the xz plane. s This represents the seismic wavelet extracted in step three, where t represents time, x represents the horizontal spatial coordinate, and z represents the vertical spatial coordinate; v x p v z p Together with p, they constitute the forward propagation wave field of the P-wave at the source end.

[0022] In step 2, using the given migration velocity model and observed seismic data, the P-wave backpropagation wavefield at the receiver end and the S-wave backpropagation wavefield at the receiver end are calculated.

[0023] In step 3, using the forward propagation wavefield of the P-wave at the source end and the reverse propagation wavefield of the P-wave at the receiver end calculated in step 2, the PP imaging data volume is calculated under vector imaging conditions; the calculation formula is as follows:

[0024]

[0025] Among them, I pp Represents the PP imaging data volume, v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x p,b v represents the longitudinal wave velocity component in the x-direction vibration of a particle in the backpropagation wave field of the P-wave at the detector end. z p,b This represents the longitudinal wave velocity component in the z-direction vibration of a particle in the backpropagation wave field of the P-wave at the detector point.

[0026] In step 3, using the forward propagation wavefield of the P-wave at the source end and the reverse propagation wavefield of the S-wave at the receiver end calculated in step 2, the PS imaging data volume is calculated under vector imaging conditions; the calculation formula is as follows:

[0027]

[0028] Among them, I ps Represents the PS imaging data volume, v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x s,b v represents the transverse wave velocity component in the x-direction vibration of a particle in the S-wave backpropagation field at the detector end. z s,b This represents the transverse wave velocity component in the z-direction vibration of a particle in the S-wave backpropagation wave field at the detector point.

[0029] In step 4, the formula for calculating the cosine similarity between the PP and PS imaging data volumes is as follows:

[0030]

[0031] Where cos represents the cosine function, θ represents the angle between the PP and PS imaging data volumes, and I pp I represents a one-dimensional vector formed by each element in the PP imaging data volume. ps represents a one-dimensional vector composed of each element in the PS imaging data volume, · represents the dot product operation, and || represents calculating the magnitude of the vector.

[0032] In step 5, the formula for filtering the PP imaging data volume is as follows:

[0033] I' pp =I pp cos(θ), (5)

[0034] Among them, I' pp I represents the filtered PP imaging data volume. pp Let θ represent the PP imaging data volume, cos represent the cosine function, and θ represent the angle between the PP and PS imaging data volumes.

[0035] In step 5, the formula for filtering the PS imaging data volume is as follows:

[0036] I' ps =I ps cos(θ), (6)

[0037] Among them, I' ps I represents the filtered PS imaging data volume.ps Let θ represent the PS imaging data volume, cos represent the cosine function, and θ represent the angle between the PP and PS imaging data volumes.

[0038] The objective of this invention can also be achieved through the following technical measures: a PP and PS joint reverse time migration imaging system based on cosine similarity metric, wherein the PP and PS joint reverse time migration imaging system based on cosine similarity metric employs a PP and PS joint reverse time migration imaging method based on cosine similarity metric to perform elastic wave migration imaging.

[0039] The objective of this invention can also be achieved through the following technical measures: a computer-readable storage medium comprising a stored computer program; wherein, when the computer program is running, it controls the device containing the computer-readable storage medium to execute a PP and PS joint reverse time migration imaging method based on cosine similarity metric.

[0040] The cosine similarity-based combined PP and PS reverse time migration imaging method of this invention eliminates artifacts in PP and PS imaging profiles, improves the overall quality of PP and PS imaging profiles, and overcomes the problem of migration artifacts affecting the accurate identification of subsurface geological structures. Using the cosine similarity-based PP and PS combined reverse time migration imaging method described in this invention, the cosine similarity of PP and PS imaging data volumes is calculated to filter the PP and PS imaging data volumes, eliminating artifacts in the PP and PS imaging data volumes, improving the overall quality of PP and PS imaging data volumes, and overcoming the adverse effects of migration artifacts on the accurate identification of subsurface geological structures. This invention is computationally simple, easy to implement, highly practical, and produces highly reliable imaging results. Attached Figure Description

[0041] Figure 1 This is a flowchart of a specific embodiment of the PP and PS joint reverse time migration imaging method based on cosine similarity measurement of the present invention;

[0042] Figure 2 This is a PP imaging profile calculated using a conventional reverse time migration imaging method in a specific embodiment of the present invention;

[0043] Figure 3 This is a PS imaging profile calculated using a conventional reverse time migration imaging method in a specific embodiment of the present invention;

[0044] Figure 4 This is a PP imaging profile calculated in a specific embodiment of the present invention;

[0045] Figure 5 This is a PS imaging profile calculated in a specific embodiment of the present invention. Detailed Implementation

[0046] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0047] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments of the present invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, and / or combinations thereof.

[0048] The cosine similarity metric-based PP and PS joint reverse time migration imaging method of the present invention includes:

[0049] Step 1: Acquire observed seismic data and extract seismic wavelets;

[0050] Step 2: Calculate the forward propagation wavefield of the P-wave at the source end, the reverse propagation wavefield of the P-wave at the receiver end, and the reverse propagation wavefield of the S-wave at the receiver end.

[0051] Step 3: Calculate the PP imaging data volume and the PS imaging data volume;

[0052] Step 4: Calculate the cosine similarity between the PP and PS imaging data volumes;

[0053] Step 5: Filter the PP imaging data volume and the PS imaging data volume using cosine similarity;

[0054] Step 6: Superimpose the filtered PP imaging data volume and PS imaging data volume corresponding to all single-shot records, and output the PP and PS imaging data volumes.

[0055] The following are several specific embodiments of the application of the present invention.

[0056] Example 1

[0057] In a specific embodiment 1 of the present invention, such as Figure 1 As shown, Figure 1 This is a flowchart of the PP and PS joint reverse time migration imaging method based on cosine similarity metric of the present invention. The PP and PS joint reverse time migration imaging method based on cosine similarity metric includes:

[0058] Step 1: Acquire the observed seismic data: Place a four-component nodal seismograph on the seabed and record the artificially generated source wave field. Use the recorded data as the observed seismic data. For the theoretical model, the observed seismic data is obtained by forward modeling the elastic wave equation with second-order time and 12th-order spatial accuracy based on the real velocity model.

[0059] Step 2: Preprocess the observed seismic data to obtain preprocessed seismic data. Preprocessing includes noise suppression, multiple suppression, node attitude correction, node position correction, and removal of direct and refracted waves.

[0060] Step 3: Using the preprocessed seismic data, extract the seismic wavelet using a higher-order statistical method;

[0061] Step 4: Using the given migration velocity model and the extracted seismic wavelet, calculate the forward propagation wavefield of the P-wave at the source end;

[0062] Step 5: Using the given migration velocity model and observed seismic data, calculate the P-wave backpropagation field at the receiver end;

[0063] Step 6: Using the given migration velocity model and observed seismic data, calculate the S-wave backpropagation field at the receiver end;

[0064] Step 7: Using the forward propagation wavefield of the P-wave at the source end calculated in Step 4 and the reverse propagation wavefield of the P-wave at the receiver end calculated in Step 5, calculate the PP imaging data volume under vector imaging conditions.

[0065] Step 8: Using the forward propagation wavefield of the P-wave at the source end calculated in Step 4 and the reverse propagation wavefield of the S-wave at the receiver end calculated in Step 6, calculate the PS imaging data volume under vector imaging conditions.

[0066] Step 9: Using the PP imaging data volume calculated in Step 7 and the PS imaging data volume calculated in Step 8, calculate the cosine similarity between the PP and PS imaging data volumes.

[0067] Step 10: Filter the PP imaging data using the cosine similarity calculated in Step 9;

[0068] Step 11: Filter the PS imaging data using the cosine similarity calculated in Step 9;

[0069] Step 12: Overlay the filtered PP imaging data volumes corresponding to all single-shot records to obtain the PP imaging data volume of the subsurface medium;

[0070] Step 13: Overlay the filtered PS imaging data volumes corresponding to all single-shot records to obtain the PS imaging data volume of the subsurface medium;

[0071] Step Fourteen: Output the superimposed PP and PS imaging data volumes.

[0072] Furthermore, in step four, the formula for calculating the forward propagation wavefield of the P-wave at the source end is as follows:

[0073]

[0074] Where ρ represents the density of the medium, k represents the bulk modulus of the medium, μ represents the shear modulus of the medium, and v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x s v represents the transverse wave velocity component of a particle vibrating in the x-direction. z s v represents the transverse wave velocity component of a particle vibrating in the z-direction. x v represents the velocity component of a particle vibrating in the x-direction. z Let τ represent the velocity component of the particle vibrating in the z-direction, p represent the pressure field, and τ represent the velocity component of the particle vibrating in the z-direction. xx τ represents the normal stress field in the x-direction. zz τ represents the normal stress field in the z-direction. xz f represents the shear stress field in the xz plane. s This represents the seismic wavelet extracted in step three, where t represents time, x represents the horizontal spatial coordinate, and z represents the vertical spatial coordinate; v x p v z p Together with p, they constitute the forward propagation wave field of the P-wave at the source end.

[0075] Furthermore, in step seven, the calculation formula for the PP imaging data volume is as follows:

[0076]

[0077] Among them, I pp Represents the PP imaging data volume, v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x p,b v represents the longitudinal wave velocity component in the x-direction vibration of a particle in the backpropagation wave field of the P-wave at the detector end. z p,b This represents the longitudinal wave velocity component in the z-direction vibration of a particle in the backpropagation wave field of the P-wave at the detector point.

[0078] Furthermore, in step eight, the calculation formula for the PS imaging data volume is as follows:

[0079]

[0080] Among them, I ps Represents the PS imaging data volume, v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x s,b v represents the transverse wave velocity component in the x-direction vibration of a particle in the S-wave backpropagation field at the detector end. z s,b This represents the transverse wave velocity component in the z-direction vibration of a particle in the S-wave backpropagation wave field at the detector point.

[0081] Furthermore, in step nine, the formula for calculating the cosine similarity of the PP and PS imaging data volumes is as follows:

[0082]

[0083] Where cos represents the cosine function, θ represents the angle between the PP and PS imaging data volumes, and I pp I represents a one-dimensional vector formed by each element in the PP imaging data volume. ps represents a one-dimensional vector composed of each element in the PS imaging data volume, · represents the dot product operation, and || represents calculating the magnitude of the vector.

[0084] Furthermore, in step ten, the calculation formula for filtering the PP imaging data volume is as follows:

[0085] I' pp =I pp cos(θ), (5)

[0086] Among them, I' pp I represents the filtered PP imaging data volume. pp Let θ represent the PP imaging data volume, cos represent the cosine function, and θ represent the angle between the PP and PS imaging data volumes.

[0087] Furthermore, in step eleven, the calculation formula for filtering the PS imaging data volume is as follows:

[0088] I' ps =I ps cos(θ), (6)

[0089] Among them, I' ps I represents the filtered PS imaging data volume. psLet θ represent the PS imaging data volume, cos represent the cosine function, and θ represent the angle between the PP and PS imaging data volumes.

[0090] Example 2

[0091] In a specific embodiment 2 of the present invention, in order to further illustrate the implementation idea and process of the present invention and prove the effectiveness of the method, a two-layer uniform velocity model was used for testing, and the results were compared with those calculated by the conventional reverse time migration imaging method.

[0092] S1: Construct a two-layer uniform velocity model as the true velocity model. The true velocity model has a width of 7.5 km and a depth of 2 km. The first layer of the true velocity model has a thickness of 1 km, a P-wave velocity of 2 km / s, a S-wave velocity of 1.3 km / s, and a density of 2.5 g / cm³. 3 In the real velocity model, the thickness of the second medium is 1 km, the longitudinal wave velocity is 3 km / s, the transverse wave velocity is 1.8 km / s, and the density is 2.6 g / cm³. 3 The sample is discretized using a square grid with a grid size of 5m.

[0093] S2: Observation System: Seismic sources are evenly distributed at 50m intervals from 0km in the X direction to the far right of the model. Nodes are evenly distributed at 10m intervals from 0km in the X direction to the far right of the model. The depth of both seismic sources and nodes is set to 5m. The seismic data sampling time is 2s, and the sampling interval is 0.5ms.

[0094] S3: Using a real velocity model and a Ricker wavelet with a dominant frequency of 15Hz, multi-wave multi-component seismic data is obtained through forward modeling of regular grid elastic wave equations with second-order time and 12th-order spatial accuracy, and with fully matched layer boundary conditions. This data is then used as the observed seismic data.

[0095] S4: Preprocess the observed seismic data, removing direct waves and refracted waves to obtain preprocessed seismic data.

[0096] S5: Use the 15Hz Reck wavelet as the seismic wavelet.

[0097] S6: Using the longitudinal wave velocity, transverse wave velocity and density of the first layer of medium in the real velocity model, the forward propagation wave field of the P-wave at the source end is calculated using formula (1).

[0098] S7: Calculate the P-wave backpropagation wave field at the receiver end using the longitudinal wave velocity, transverse wave velocity, and density of the first layer of medium in the real velocity model.

[0099] S8: Calculate the S-wave backpropagation wave field at the receiver end using the longitudinal wave velocity, transverse wave velocity, and density of the first layer of medium in the real velocity model.

[0100] S9: Using the forward propagation wave field of the P-wave at the source end calculated in step six and the reverse propagation wave field of the P-wave at the receiver end calculated in step seven, the PP imaging data volume is calculated using formula (2).

[0101] S10: Using the forward propagation wave field of the P-wave at the source end calculated in step six and the reverse propagation wave field of the S-wave at the receiver end calculated in step eight, the PS imaging data volume is calculated using formula (3).

[0102] S11: Using the PP imaging data volume calculated in step nine and the PS imaging data volume calculated in step ten, calculate the cosine similarity between the PP and PS imaging data volumes using formula (4).

[0103] S12: Using the cosine similarity calculated in step eleven, filter the PP imaging data volume calculated in step nine using formula (5) to obtain the filtered PP imaging data volume.

[0104] S13: Using the cosine similarity calculated in step eleven, filter the PS imaging data volume calculated in step ten using formula (6) to obtain the filtered PS imaging data volume.

[0105] S14: Superimpose the filtered PP imaging data volumes corresponding to all single-shot records to obtain the PP imaging data volume of the subsurface medium.

[0106] S15: Overlay the filtered PS imaging data volumes corresponding to all single-shot records to obtain the PS imaging data volume of the subsurface medium.

[0107] S16: Output the superimposed PP and PS imaging data volumes.

[0108] Figure 2 and Figure 3 These are the PP and PS imaging profiles calculated using the conventional reverse time migration imaging method. Figure 4 and Figure 5 These are the PP and PS imaging profiles calculated by this invention. From Figure 2 and Figure 3 As can be seen, beyond the velocity interface at a depth of 1 km, there are visible in-phase axes. The depth of these in-phase axes is not equal to 1 km; they are illusions. Figure 4 and Figure 5 In the study, in-phase axes were found only at the velocity interface at a depth of 1 km, and not in other regions. This indicates that the present invention can effectively eliminate artifacts calculated by conventional reverse time migration imaging methods, improve the overall quality of PP and PS imaging data volumes, and overcome the adverse effects of imaging artifacts on subsequent seismic data interpretation.

[0109] Finally, it should be noted that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0110] Except for the technical features described in the specification, all other technologies are known to those skilled in the art.

Claims

1. A joint reverse time migration imaging method for PP and PS based on cosine similarity metric, characterized in that, This joint reverse time migration imaging method based on cosine similarity metric (PP) and PS includes: Step 1: Acquire observed seismic data and extract seismic wavelets; Step 2: Calculate the forward propagation wavefield of the P-wave at the source end, the reverse propagation wavefield of the P-wave at the receiver end, and the reverse propagation wavefield of the S-wave at the receiver end. Step 3: Calculate the PP imaging data volume and the PS imaging data volume; Step 4: Calculate the cosine similarity between the PP and PS imaging data volumes; Step 5: Filter the PP imaging data volume and the PS imaging data volume using cosine similarity; Step 6: Superimpose the filtered PP imaging data volume and PS imaging data volume corresponding to all single-shot records, and output the PP and PS imaging data volumes.

2. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 1, characterized in that, In step 1, the observed seismic data is acquired by placing a four-component nodal seismograph on the seabed and recording the artificially generated source wave field. The recorded data is used as the observed seismic data. For the theoretical model, the observed seismic data is obtained by forward modeling the elastic wave equation based on the true velocity model.

3. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 2, characterized in that, In step 1, after acquiring the observed seismic data, the observed seismic data is preprocessed to obtain preprocessed seismic data. The preprocessing includes noise suppression, multiple wave suppression, node attitude correction, node position correction, and removal of direct and refracted waves.

4. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 3, characterized in that, In step 1, the seismic wavelet is extracted using the preprocessed seismic data and a higher-order statistical method.

5. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 1, characterized in that, In step 2, using the given migration velocity model and the extracted seismic wavelet, the forward propagation wavefield of the P-wave at the source end is calculated. The calculation formula is as follows: Where ρ represents the density of the medium, k represents the bulk modulus of the medium, μ represents the shear modulus of the medium, and v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x s v represents the transverse wave velocity component of a particle vibrating in the x-direction. z s v represents the transverse wave velocity component of a particle vibrating in the z-direction. x v represents the velocity component of a particle vibrating in the x-direction. z Let τ represent the velocity component of the particle vibrating in the z-direction, p represent the pressure field, and τ represent the velocity component of the particle vibrating in the z-direction. xx τ represents the normal stress field in the x-direction. zz τ represents the normal stress field in the z-direction. xz f represents the shear stress field in the xz plane. s This represents the seismic wavelet extracted in step three, where t represents time, x represents the horizontal spatial coordinate, and z represents the vertical spatial coordinate; v x p v z p Together with p, they constitute the forward propagation wave field of the P-wave at the source end.

6. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 1, characterized in that, In step 2, using the given migration velocity model and observed seismic data, the P-wave backpropagation wavefield at the receiver end and the S-wave backpropagation wavefield at the receiver end are calculated.

7. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 1, characterized in that, In step 3, using the forward propagation wavefield of the P-wave at the source end and the reverse propagation wavefield of the P-wave at the receiver end calculated in step 2, the PP imaging data volume is calculated under vector imaging conditions. The calculation formula is as follows: Among them, I pp Represents the PP imaging data volume, v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x p,b v represents the longitudinal wave velocity component in the x-direction vibration of a particle in the backpropagation wave field of the P-wave at the detector end. z p,b This represents the longitudinal wave velocity component in the z-direction vibration of a particle in the backpropagation wave field of the P-wave at the detector point.

8. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 1, characterized in that, In step 3, using the forward propagation wavefield of the P-wave at the source end and the reverse propagation wavefield of the S-wave at the receiver end calculated in step 2, the PS imaging data volume is calculated under vector imaging conditions; the calculation formula is as follows: Among them, I ps Represents the PS imaging data volume, v x p v represents the longitudinal wave velocity component of a particle vibrating in the x-direction. z p V represents the longitudinal wave velocity component of a particle vibrating in the z-direction. x s,b v represents the transverse wave velocity component in the x-direction vibration of a particle in the S-wave backpropagation field at the detector end. z s,b This represents the transverse wave velocity component in the z-direction vibration of a particle in the S-wave backpropagation wave field at the detector point.

9. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 1, characterized in that, In step 4, the formula for calculating the cosine similarity between the PP and PS imaging data volumes is as follows: Where cos represents the cosine function, θ represents the angle between the PP and PS imaging data volumes, and I pp I represents a one-dimensional vector formed by each element in the PP imaging data volume. ps represents a one-dimensional vector composed of each element in the PS imaging data volume, · represents the dot product operation, and || represents calculating the magnitude of the vector.

10. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 1, characterized in that, In step 5, the formula for filtering the PP imaging data volume is as follows: I' pp =I pp cos(θ), (5) Among them, I' pp I represents the filtered PP imaging data volume. pp Let θ represent the PP imaging data volume, cos represent the cosine function, and θ represent the angle between the PP and PS imaging data volumes.

11. The PP and PS joint reverse time migration imaging method based on cosine similarity metric according to claim 1, characterized in that, In step 5, the formula for filtering the PS imaging data volume is as follows: I' ps =I ps cos(θ), (6) Among them, I' ps I represents the filtered PS imaging data volume. ps Let θ represent the PS imaging data volume, cos represent the cosine function, and θ represent the angle between the PP and PS imaging data volumes.

12. A joint reverse time migration imaging system for PP and PS based on cosine similarity metric, characterized in that, The PP and PS joint reverse time migration imaging system based on cosine similarity metric is used for elastic wave migration imaging according to any one of claims 1-11.

13. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored computer program; wherein, when the computer program is executed, it controls the device on which the computer-readable storage medium is located to perform a PP and PS joint reverse time migration imaging method based on cosine similarity metric as described in any one of claims 1-11.

Images