Decoupling Poynting vector based on divergence and curl operator and application method and system of decoupling Poynting vector in elastic wave reverse time migration

Abstract

The invention belongs to the technical field of elastic wave reverse time migration, and particularly relates to a decoupling Poynting vector based on divergence and rotation operators and an application method of the decoupling Poynting vector in elastic wave reverse time migration, and the method comprises the following steps: P wave and S wave decoupling and Poynting vector calculation based on the divergence and rotation operators; the invention relates to an elastic wave RTM imaging condition based on decoupling Poynting vector constraint. According to the method, the Poynting vector calculation method compatible with divergence / rotation wave field decoupling is developed theoretically, the comprehensive advantages of the method in the aspects of improving the elastic wave RTM imaging quality (especially the reliability of PS wave imaging), effectively suppressing noise and improving the calculation efficiency are shown in the practical application level, and the method has important scientific value and wide application prospects and is worthy of popularization and application. And an effective new way is provided for further improving the practical level and the imaging precision of the elastic wave RTM under the existing computing resource condition.

Description

Technical Field

[0001] This invention belongs to, but is not limited to, the field of elastic wave reverse time migration technology, and particularly relates to a decoupled Poynting vector based on divergence and curl operators and its application method and system in elastic wave reverse time migration. Background Technology

[0002] Elastic wave reverse time migration (RTM) is an important method for acquiring high-precision images of complex geological structures. In elastic wave RTM, decoupling of the P and S waves before imaging is crucial to effectively avoid the impact of crosstalk between P and S waves on the imaging results. Simultaneously, using the Poynting vector to indicate the direction of wavefield energy propagation is key to suppressing migration noise and artifacts and improving image quality. However, existing technologies face the following core technical challenges in this area: 1. Ambiguity in the representation of P- and S-wave directions using traditional Poynting vectors: The Poynting vector calculated by traditional methods through the product of the stress tensor and the particle velocity vector is a mixture of P- and S-wave energy flows. When the P- and S-waves propagate in different directions, this mixed Poynting vector cannot accurately represent the propagation directions of both waves simultaneously. This limits its potential for application in precisely controlling the imaging process and effectively suppressing noise and artifacts associated with specific wave types.

[0003] 2. Limitations of existing Poynting vector calculation methods for P- and S-wave decoupling: (1) Some proposed methods for decoupling P and S waves in Poynting vector calculation either rely on modifying the original wave equation (such as introducing a decoupling wave equation form), which usually increases additional wave field variables and differential operators, thereby increasing computational complexity and memory requirements, and may make it difficult to apply boundary conditions such as free surfaces.

[0004] (2) Other methods may not be pure enough when separating the S-wave Poynting vector. The calculated S-wave energy flux density will still be contaminated by the residual P-wave component, affecting the accuracy of its application.

[0005] 3. Lack of a mature decoupling Poynting vector calculation approach compatible with efficient and concise wavefield separation methods: Methods for directly separating P- and S-wave modes in elastic wavefields based on divergence and curl operators still have significant practical value in elastic RTM due to their advantages of simple implementation, low computational requirements, and no change to the original wave equation form. However, there is currently a lack of a mature and accurate decoupling P- and S-wave Poynting vector calculation theory and practical procedure directly compatible with such efficient wavefield separation methods.

[0006] 4. The Challenge of S-Wave Polarity Correction in Multi-Component (Especially PS-Wave) RTM Imaging: PS-converted wave imaging is crucial for acquiring more comprehensive subsurface information, but the S-wave polarity reversal problem during its imaging process has always been a technical challenge affecting image quality. If the propagation direction information of the S-wave cannot be accurately obtained (which is closely related to the accuracy of the Poynting vector), traditional imaging conditions cannot effectively solve the S-wave polarity problem, easily leading to loss of effective reflected energy, phase axis misalignment, or blurring in the imaging results, thereby reducing the reliability of the imaging results.

[0007] The closest prior art to this invention is the "Source-Free P-SVConverted-WaveReverse-TimeMigrationUsingFirst-OrderVelocity-Dilatation-RotationEquations" proposed by He et al. in 2022. This method rewrites the wavefield extension stage into a first-order system of velocity-volume strain-rotation equations. During wavefield propagation, P / S modes are automatically decoupled, and the Poynting vectors of pure P-waves and pure S-waves can be directly calculated to distinguish uplink and downlink energies. Based on this, imaging conditions for source-receiver wavefield cross-correlation are constructed, thereby obtaining PP and PS reflection imaging results. Compared to the traditional approach of combining velocity-stress equations with divergence / curl separation, this scheme is representative and advanced in its integration of pure wave direction identification and converted wave imaging.

[0008] However, the aforementioned techniques still suffer from three prominent problems: First, the velocity-volume strain-rotation equation requires additional storage of strain and rotation variables, expanding the dimensions from 6 velocity-stress components to more than 9, further increasing the already high computational and memory overhead of elastic RTM, making it at least an order of magnitude higher than acoustic RTM in 2D / 3D scenarios; Second, the paper mainly uses zero-delay cross-correlation imaging without introducing explicit angle or direction constraints, resulting in insufficient suppression of multiple waves and bounce energy, and a tendency for low-frequency crosstalk and artifacts; Third, the authors acknowledge that when the structure is complex, the S-wave polarity is difficult to accurately correct, and local fragmentation remains significant, requiring subsequent manual or empirical processing to improve image continuity. These shortcomings are precisely the technical problems that this invention aims to solve by using a lightweight decoupling of the Poynting vector and angle threshold imaging conditions based on divergence / curl. Summary of the Invention

[0009] This invention aims to solve the aforementioned technical problems. Its core lies in proposing a new Poynting vector calculation formula for P-waves and S-waves that is compatible with wavefield decoupling methods based on divergence and curl operators. This method allows for the accurate determination of the propagation directions of P-waves and S-waves without altering the original wave equations, and can be applied to the imaging conditions of elastic wave RTM. Its purpose is to more effectively suppress migration noise and imaging artifacts, improve the clarity and resolution of the imaging profile, and especially achieve accurate polarity correction for S-waves, thereby comprehensively improving the imaging quality and practicality of multi-component elastic wave RTM while maintaining high computational efficiency.

[0010] This invention is implemented as follows: a method for decoupling the Poynting vector based on divergence and curl operators and its application in the reverse-time migration of elastic waves, the method comprising: S1: Decoupling of P-waves and S-waves and Poynting vector calculation based on divergence and curl operators.

[0011] S2: Elastic wave RTM imaging conditions based on decoupled Poynting vector constraints.

[0012] Furthermore, S1 specifically includes: First, the wavefield extension in RTM is based on the first-order velocity-stress elastic wave equation:

[0013] in It is a velocity wave field The Quantity, These are stress components, and the point at the top represents the time derivative. and It is Lamé's constant. It is the density of the medium. It uses Kronecker notation, the physical term is omitted, and Einstein's summation convention is adopted.

[0014] According to elastic wave theory, P-waves in isotropic media are irrotational, with their polarization direction parallel to the propagation direction; S-waves are pressureless, with their polarization direction perpendicular to the propagation direction. Therefore, P-waves and S-waves can be decoupled by directly calculating the divergence and curl of the wave field, which is one of the most commonly used decoupling methods in active source RTM.

[0015] in It is the Nabla operator. Denotes divergence, Indicates curl. and These represent the decoupled scalar P-wave and vector S-wave (in the two-dimensional case), respectively. Degenerate into a scalar The main advantage of this decoupling method is that all the differential terms required to calculate divergence and curl have been pre-calculated when solving equation (1), so no additional differential calculation is required during decoupling, only simple addition and subtraction of the calculated terms are involved, and the amount of calculation is negligible compared with other decoupling methods. The decoupling Poynting vector associated with the traditional decoupling wave equation has the following form:

[0016] in and These represent the velocities of the P-wave and S-wave in the medium, respectively, and the volumetric strain is also present. It is a rotation vector. and These are the decoupled velocity vectors, and these physical quantities satisfy:

[0017] Since longitudinal waves are irrotational fields ,Right now Transverse waves are non-scattered. ,Right now Taking the time derivative of both sides of equation (3), and exchanging the time and space derivatives, we can simplify to obtain:

[0018] in and These are the decoupled accelerated wave field vectors. It can be seen that the physical quantities... , , and The equation (4) is satisfied, and , , and The equation (3) satisfies the same form; therefore, by analogy with equation (2), we can write the wave field after decoupling through divergence and curl. and Poynting vector:

[0019] Based on the relationship between the wave field velocity vector and the acceleration vector, and combined with equation (5a), we can obtain:

[0020] In an isotropic medium, volume strain It can be represented as:

[0021] In the three-dimensional case In two-dimensional case That is, volumetric strain. It can be directly obtained by summing all the normal stresses and then multiplying by . We obtain it. Then, according to formula (6), we perform... By calculating the gradient, the decoupled acceleration P-wave field can be obtained. This is the only place in this paper where the Poynting vector calculation requires difference calculation.

[0022] The accelerating wave field of the transverse wave can be directly obtained by subtracting the accelerating wave field of the longitudinal wave from the overall accelerating wave field; the overall accelerating wave field has also been obtained when solving the first-order velocity-stress equation (1a), and its components are... Written in vector form as The accelerated wave field of the transverse wave. It can be written as:

[0023] Will , Substituting into formula (6), the corresponding decoupling Poynting vector can be obtained.

[0024] Furthermore, S2 specifically includes: After completing the forward continuation of the source wavefield and the reverse continuation of the receiver wavefield, imaging is performed using the decoupled P-waves and S-waves obtained during the continuation process, applying appropriate imaging conditions. In conventional elastic wave RTM, imaging primarily relies on first-arrival reflection signals from the subsurface. During RTM, the reflected waves generated by the forward and reverse continuation wavefields propagate in opposite directions, therefore the angle between the two wavefields is typically obtuse. By constraining this angle during imaging, migration noise and artifacts can be effectively suppressed. Using pi / 2 as the constraint angle allows for the direct calculation of stable and efficient imaging decisions based on the sign of the Poynting vector dot product of the source and receiver wavefields, without the need to calculate an explicit inverse cosine function.

[0025] In the two-dimensional case, the normalized cross-correlation imaging condition based on angle constraints can be expressed as:

[0026] in P-wave from the source wave field P-waves of the receiving point wave field The result obtained is called the PP imaging result; P-wave from the source wave field S-wave of the receiving point wave field The result obtained is called the PS imaging result; yes and The angle between the directions of propagation yes and The angle between the propagation directions. The propagation direction of the wave field can be represented by the Poynting vector; therefore, combining the calculated decoupled P-wave and S-wave Poynting vectors, equation (10) can be equivalently written as:

[0027] in It is the P-wave Poynting vector of the source wave field. It is the P-wave Poynting vector of the received wave field. The S-wave Poynting vector of the received wave field can be calculated using equation (6).

[0028] On the other hand, the S-wave decoupled by the curl operator exhibits polarity reversal, which leads to inconsistent polarity of the imaging results at different locations when superimposing different single-shot PS imaging results, resulting in energy cancellation during the superposition process. Therefore, polarity correction of the S-wave is required during imaging. A commonly used polarity correction method in elastic wave RTM is to calculate the directional relationship between the incident and reflected wave fields using the Poynting vector, and then correct the S-wave polarity. For PS wave RTM results, the polarity correction method can be expressed as:

[0029] in This is the PS migration result after S-wave correction.

[0030] Based on the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solution to be protected by this invention are as follows: First, in conventional seismic imaging and wavefield visualization workflows, the Poynting vector is widely used to determine the wavefield propagation direction, assess wavelet illumination quality, and guide key steps such as reverse time migration (RTM). However, traditional algorithms typically perform cross-products or dot-products between the full vector velocity field and the stress tensor, leading to cross-contamination of P-wave and S-wave energy. The pseudo-energy generated by interface mode conversion is amplified in complex structures, introducing stripe noise, azimuth misjudgment, and amplitude distortion into cross-sectional imaging. This limits the resolution and reliability of industrial applications such as deep oil and gas exploration, mineral exploration, and urban active fault location.

[0031] This invention, within the framework of the first-order velocity-stress elastic wave equation, first uses divergence and curl operators to rigorously decompose the velocity field into P-wave and S-wave components. Then, it further sums the normal stress components to obtain the volumetric strain and calculates the gradient, constructing the P-wave acceleration field. Subtracting the P-wave acceleration from the overall acceleration yields the S-wave acceleration field. Using this approach, pure mode separation is achieved simultaneously at both the source and receiver ends. Furthermore, locally differentiable operations can yield physically clear and phase-consistent acceleration information, laying a pure foundation for the subsequent construction of the P-wave / S-wave Poynting vector in the form of "acceleration × velocity".

[0032] Thanks to the approach of decoupling the wavefield before calculating the flux, the generated P and S energy flux direction vectors no longer exhibit polarity reversal and crosstalk near the interface. The energy amplitude and incident angle maintain a strict physical correspondence, significantly reducing imaging artifacts caused by pseudo-scattering, boundary reflection, and numerical dispersion. In a large-scale GPU parallel environment, this algorithm only involves first-order spatial gradient and tensor summation. Compared to traditional methods that require additional Fourier filtering or higher-order pseudospectral calculations, it can save 20–30% of GPU memory and computation time, providing an efficient solution for industrial scenarios with extremely high real-time requirements, such as deep-water multi-component exploration and pre-drilling prediction of ultra-deep wells.

[0033] At the industry level, pure modal Poynting vectors not only improve the accuracy of illumination analysis in RTM, full waveform inversion (FWI), and spread spectrum wavelet design, but also provide more reliable quantitative indicators of energy flow for evaluating fracture strike, fracture density, and reservoir anisotropy. Combined with existing integrated seismic acquisition and processing platforms, this invention reduces subsurface structural interpretation errors by more than 10° and improves thin-layer resolution by more than 15%, directly resulting in significant economic benefits such as optimized drilling targets, reduced exploration risks, and savings in development costs. This demonstrates a substantial improvement and significant progress compared to existing technologies.

[0034] Compared with existing methods that mainly employ traditional coupled Poynting vectors or other complex decoupled Poynting vector methods that rely on specific decoupled wave equations, the decoupled Poynting vector based on divergence and curl operators proposed in this invention and its application in RTM have the following significant advantages: 1. Significantly improves image quality and signal-to-noise ratio, effectively suppressing noise and artifacts.

[0035] (1) By accurately separating the P-wave and S-wave Poynting vectors and applying them to imaging conditions, the energy propagation direction can be constrained more precisely, thereby more effectively suppressing low-frequency noise and various imaging artifacts generated during RTM.

[0036] (2) Both PP wave and PS wave imaging can obtain clearer geological interfaces, richer structural details and higher imaging resolution, especially in complex structural areas (such as steeply dipped interfaces), where the imaging effect is significantly improved, as shown in the test results of graben model and Marmousi model. Figure 6 , Figure 8 , Figure 9 For PP wave imaging, decoupled Poynting vectors can more effectively suppress low-frequency noise and offset artifacts compared to traditional coupled Poynting vectors.

[0037] 2. It efficiently and accurately solves the problem of S-wave polarity correction in PS converted wave imaging.

[0038] (1) The present invention utilizes the accurately calculated decoupled S-wave Poynting vector direction information to reliably correct the polarity of PS-wave reflection (as shown in Equation (12)), effectively avoiding imaging energy loss, phase axis misalignment and interface distortion caused by S-wave polarity reversal.

[0039] (2) This makes the PS wave imaging results more reliable and easier to interpret, which is of great significance for the comprehensive use of P-wave and S-wave information for reservoir prediction and lithology identification. Even after multi-shot data are superimposed, its advantages remain significant (e.g. Figure 6 Part (f) in the middle, Figure 9 (part (b) of the text).

[0040] 3. It is highly compatible with widely used efficient wave field separation methods, easy to implement, and does not change the original wave equation.

[0041] (1) The decoupled Poynting vector calculation proposed in this method is closely integrated with the P and S wave field separation method based on divergence and curl operators widely used in the field of geophysics, without the need to make fundamental modifications to the wave equation solver (such as the first-order velocity-stress equation) in the existing RTM process.

[0042] (2) It fills the gap in the lack of mature and accurate matching decoupling Poynting vector calculation for this type of efficient wave field separation method, and provides a practical technical approach.

[0043] 4. It has advantages in computing efficiency, while its memory requirements are comparable.

[0044] (1) Compared with some decoupled Poynting vector methods that require solving additional differential equations or introducing complex auxiliary variables (such as methods based on modified specific decoupled wave equations), the method of this invention is more computationally efficient because its Poynting vector calculation process mainly relies on algebraic operations and limited spatial differentiation (only required when calculating the volume strain gradient).

[0045] (2) Marmousi model tests show that this method improves computation speed by about 7% compared to an improved traditional decoupled Poynting vector method, while maintaining essentially the same computer memory requirements. Figure 10 ).

[0046] 5. Enhance the imaging capabilities and reliability of complex geological structures.

[0047] For geological targets such as steeply dipping structures, faults, and thin interbedded layers, this method, through more precise energy direction constraints and S-wave polarity correction, can provide more accurate, focused, and reliable imaging results, which is beneficial for detailed geological interpretation and subsequent oil and gas exploration and development decisions.

[0048] 6. It has good potential for three-dimensional expansion.

[0049] The decoupled Poynting vector calculation formula and imaging method proposed in this invention are not limited to two-dimensional media. Their theoretical basis and calculation process can be directly applied to three-dimensional elastic wave RTM, and are expected to play a positive role in the scalarization processing of three-dimensional vector S-waves.

[0050] This invention not only theoretically develops a Poynting vector calculation method compatible with divergence / curl wave field decoupling, but also demonstrates its comprehensive advantages in improving the imaging quality of elastic wave RTM (especially the reliability of PS wave imaging), effectively suppressing noise, and improving computational efficiency in practical applications. It has important scientific value and broad application prospects, and provides an effective new approach to further improve the practicality and imaging accuracy of elastic wave RTM under existing computing resource conditions.

[0051] Second, the technical solution of this invention fills the technical gap in the field of elastic wave RTM, which requires accurate and efficient decoupling Poynting vector calculation that is compatible with the most widely used divergence and curl wave field separation methods.

[0052] In elastic wave real-time measurement (RTM), the separation of P-waves and S-waves using divergence and curl operators is widely adopted due to its efficiency, simplicity, and the fact that it does not require modification of the wave equation. However, a long-standing technical bottleneck is the lack of a mature decoupling Poynting vector computation theory that can seamlessly integrate with this separation method.

[0053] Existing technologies either use a mixed Poynting vector that cannot accurately distinguish the propagation directions of P and S waves, limiting its role in suppressing imaging artifacts; or employ other complex decoupling methods that require modification of the wave equation or the introduction of a large amount of additional computation, sacrificing computational efficiency and practicality. Therefore, how to accurately obtain the energy propagation directions of P and S waves while maintaining the efficiency of the divergence / curl method has become an unresolved technological gap.

[0054] This invention presents, for the first time, an innovative Poynting vector calculation formula compatible with divergence and curl operator decoupling methods. Through ingenious mathematical analogy and derivation of the extension results of the first-order velocity-stress equation, it successfully calculates the pure P-wave and S-wave Poynting vectors without adding the burden of solving additional wave equations. This effectively fills the gap in the industry's efficient wavefield separation technology route, which lacks a corresponding high-precision, low-cost Poynting vector calculation method, providing a new and practical technical approach to improving RTM imaging quality.

[0055] The technical solution of this invention solves the technical problem of S-wave polarity reversal correction, which has long existed and seriously affected the imaging quality in multi-component (especially PS converted wave) elastic wave reverse time migration, based on the decoupled Poynting vector of divergence and curl operators.

[0056] PS-converted wave imaging is crucial for acquiring comprehensive subsurface information. However, S-waves undergo complex polarity reversals during reflection, causing reflected energy from the same geological interface to cancel each other out due to polarity inconsistencies during conventional RTM imaging and multi-shot stacking. This results in discontinuities in the effective reflection phase axis, energy loss, interface blurring, and even distortion in the imaging results. Accurate polarity correction of S-waves has been a long-standing technical challenge for geophysicists, one they have yet to fully resolve. The fundamental reason is that accurate polarity correction requires precise S-wave propagation direction information, and traditional techniques struggle to provide pure, reliable S-wave Poynting vectors, leading to poor correction results.

[0057] This invention, through its unique decoupled Poynting vector calculation method, can capture the true propagation direction of S-waves with unprecedented accuracy. Based on this, this invention designs a novel polarity correction method based on the sign of the decoupled Poynting vector cross product. This method can reliably determine the reflection polarity of the PS wave and accurately correct it during imaging. Figure 6 As shown in part (f), the PS-wave imaging profile corrected and superimposed by the method of this invention shows a continuous and clear geological interface, while the results of traditional methods ( Figure 6The (e) portion exhibits significant energy loss and discontinuity. Therefore, this invention successfully overcomes the industry-recognized technical challenge of S-wave polarity correction in PS-wave RTM imaging. Attached Figure Description

[0058] Figure 1 This is a flowchart of the method for decoupling the Poynting vector based on divergence and curl operators and its application in the reverse time migration of elastic waves, provided by an embodiment of the present invention.

[0059] Figure 2 This is a schematic diagram of a double-layer graben model provided in an embodiment of the present invention. Figure 2 Part (a) in the figure is used for the longitudinal and transverse wave velocities of the real medium model in forward modeling; Figure 2 Part (b) is used for the longitudinal and transverse wave velocities of the smooth medium model for wavefield inverse time delay topology and imaging; Figure 2 Part (c) describes the observation system for the five shot points used in this test.

[0060] Figure 3 This is a snapshot of the wavefield before and after decoupling the wavefield at the third shot point in the double-layer graben model provided in this embodiment of the invention. The wavefield energy has been normalized, and the P and S waves have been completely decoupled. Figure 3 (a) A snapshot of the vx component of the wavefield before decoupling in part; Figure 3 (b) is a snapshot of the vz component of the decoupled wavefield. Figure 3 A snapshot of the P-wave field after decoupling in part (c); Figure 3 A snapshot of the S-wave field after decoupling part (d) in the image.

[0061] Figure 4 This is a snapshot of the Poynting vector of the third shot point in the double-graben model provided in this embodiment of the invention before and after decoupling. After decoupling, the Poynting vectors of the P-wave and S-wave are completely separated, and each accurately represents the propagation direction of the corresponding wave. Figure 4 (a) The x-component of the Poynting vector before decoupling; Figure 4 The x-component of the Poynting vector of the P-wave after decoupling in part (b); Figure 4 The x-component of the S-wave Poynting vector after decoupling in part (c); Figure 4 The z-component of the Poynting vector before decoupling in part (d); Figure 4 The z-component of the Poynting vector of the P-wave after decoupling of part (e) in the figure; Figure 4 The z-component of the S-wave Poynting vector after decoupling part (f) in the figure.

[0062] Figure 5The image shows the single-shot RTM imaging result of the third shot point in the double-layer graben model provided in this embodiment of the invention. Without Poynting vector or traditional method constraints, obvious imaging artifacts are present at the red circles and arrows. After applying Poynting vector constraints, the imaging artifacts are significantly suppressed. Figure 5 Part (a) shows PP wave imaging results that did not use Poynting vector constraints; Figure 5 Part (b) uses PP wave imaging results with conventional Poynting vector constraints; Figure 5 Part (c) uses the PP wave imaging results of the decoupled P-wave and S-wave Poynting vector constraints in this invention; Figure 5 Part (d) in the figure shows PS wave imaging results that do not use Poynting vector constraints; Figure 5 Part (e) uses PS-wave imaging results with conventional Poynting vector constraints; Figure 5 Part (f) uses the PS-wave imaging results of the decoupled P-wave and S-wave Poynting vector constraints in this invention.

[0063] Figure 6 The above is the RTM stacking result of five shot points in the double-layer graben model provided in this embodiment of the invention. After applying the Poynting vector constraint proposed in this invention, the polarity of the S-wave is accurately corrected, the offset noise is suppressed, and the imaging quality of the stacked offset profile is significantly improved. Figure 6 Part (a) shows PP wave imaging results that did not use Poynting vector constraints; Figure 6 Part (b) uses PP wave imaging results with conventional Poynting vector constraints; Figure 6 Part (c) uses the PP wave imaging results of the decoupled P-wave and S-wave Poynting vector constraints in this invention; Figure 6 Part (d) in the figure shows PS wave imaging results that do not use Poynting vector constraints; Figure 6 Part (e) uses PS-wave imaging results with conventional Poynting vector constraints; Figure 6 Part (f) uses the PS-wave imaging results of the decoupled P-wave and S-wave Poynting vector constraints in this invention.

[0064] Figure 7 This is a schematic diagram of the Marmousi model provided in an embodiment of the present invention. Figure 7 Part (a) in the text is used for the accurate P-wave velocity model in numerical simulation; Figure 7 Part (b) is used for the smooth longitudinal wave velocity model of wavefield inverse time extension and RTM.

[0065] Figure 8The image shows the PP wave migration imaging results of the Marmousi model provided in this embodiment of the invention. The red circle indicates that the method of this invention displays a clearer imaging interface in the steep tilt region. Figure 8 Part (a) in the image uses the superimposed imaging results of traditional Poynting vector constraints; Figure 8 Part (b) uses the superimposed imaging results of the decoupled P-wave and S-wave Poynting vector constraints in this invention.

[0066] Figure 9 The image shows the PS-wave migration imaging results of the Marmousi model provided in this embodiment of the invention. The PS-wave migration results with decoupled Poynting vector constraints proposed in this invention show significantly clearer imaging of shallow reflective interfaces (indicated by the red circle). Figure 9 Part (a) in the image uses the superimposed imaging results of traditional Poynting vector constraints; Figure 9 Part (b) uses the superimposed imaging results of the decoupled P-wave and S-wave Poynting vector constraints in this invention.

[0067] Figure 10 This is a comparison of the calculation time of RTM imaging provided by the embodiment of the present invention with that of traditional methods. Detailed Implementation

[0068] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0069] This invention provides a method for decoupling the Poynting vector based on divergence and curl operators and its application in the reverse-time migration of elastic waves. The method includes: S1: Decoupling of P-waves and S-waves and Poynting vector calculation based on divergence and curl operators.

[0070] S2: Elastic wave RTM imaging conditions based on decoupled Poynting vector constraints.

[0071] S1 specifically includes: First, the wavefield extension in RTM is based on the first-order velocity-stress elastic wave equation:

[0072] in It is a velocity wave field The Quantity, These are stress components, and the point at the top represents the time derivative. and It is Lamé's constant. It is the density of the medium. It uses Kronecker notation, the physical term is omitted, and Einstein's summation convention is adopted.

[0073] According to elastic wave theory, P-waves in isotropic media are irrotational, with their polarization direction parallel to the propagation direction; S-waves are pressureless, with their polarization direction perpendicular to the propagation direction. Therefore, P-waves and S-waves can be decoupled by directly calculating the divergence and curl of the wave field, which is one of the most commonly used decoupling methods in active source RTM.

[0074] in It is the Nabla operator. Denotes divergence, Indicates curl. and These represent the decoupled scalar P-wave and vector S-wave (in the two-dimensional case), respectively. Degenerate into a scalar The main advantage of this decoupling method is that all the differential terms required to calculate divergence and curl have been pre-calculated when solving equation (1), so no additional differential calculation is required during decoupling, only simple addition and subtraction of the calculated terms are involved, and the amount of calculation is negligible compared with other decoupling methods.

[0075] The decoupling Poynting vector associated with the traditional decoupling wave equation has the following form:

[0076] in and These represent the velocities of the P-wave and S-wave in the medium, respectively, and the volumetric strain is also present. It is a rotation vector. and These are the decoupled velocity vectors, and these physical quantities satisfy:

[0077] Since longitudinal waves are irrotational fields ,Right now Transverse waves are non-scattered. ,Right now Taking the time derivative of both sides of equation (3), and exchanging the time and space derivatives, we can simplify to obtain:

[0078] in and These are the decoupled accelerated wave field vectors. It can be seen that the physical quantities... , , and The equation (4) is satisfied, and , , and The equation (3) satisfies the same form; therefore, by analogy with equation (2), we can write the wave field after decoupling through divergence and curl. and Poynting vector:

[0079] Based on the relationship between the wave field velocity vector and the acceleration vector, and combined with equation (5a), we can obtain:

[0080] In an isotropic medium, volume strain It can be represented as:

[0081] In the three-dimensional case In two-dimensional case That is, volumetric strain. It can be directly obtained by summing all the normal stresses and then multiplying by . We obtain it. Then, according to formula (6), we perform... By calculating the gradient, the decoupled acceleration P-wave field can be obtained. This is the only place in this paper where the Poynting vector calculation requires difference calculation.

[0082] The accelerating wave field of the transverse wave can be directly obtained by subtracting the accelerating wave field of the longitudinal wave from the overall accelerating wave field; the overall accelerating wave field has also been obtained when solving the first-order velocity-stress equation (1a), and its components are... Written in vector form as The accelerated wave field of the transverse wave. It can be written as:

[0083] Will , Substituting into formula (6), the corresponding decoupling Poynting vector can be obtained.

[0084] S2 specifically includes: After completing the forward continuation of the source wavefield and the reverse continuation of the receiver wavefield, imaging is performed using the decoupled P-waves and S-waves obtained during the continuation process, applying appropriate imaging conditions. In conventional elastic wave RTM, imaging primarily relies on first-arrival reflection signals from the subsurface. During RTM, the reflected waves generated by the forward and reverse continuation wavefields propagate in opposite directions, therefore the angle between the two wavefields is typically obtuse. By constraining this angle during imaging, migration noise and artifacts can be effectively suppressed. Using pi / 2 as the constraint angle allows for the direct calculation of stable and efficient imaging decisions based on the sign of the Poynting vector dot product of the source and receiver wavefields, without the need to calculate an explicit inverse cosine function.

[0085] In the two-dimensional case, the normalized cross-correlation imaging condition based on angle constraints can be expressed as:

[0086] in P-wave from the source wave field P-waves of the receiving point wave field The result obtained is called the PP imaging result; P-wave from the source wave field S-wave of the receiving point wave field The result obtained is called the PS imaging result; yes and The angle between the directions of propagation yes and The angle between the propagation directions. The propagation direction of the wave field can be represented by the Poynting vector; therefore, combining the calculated decoupled P-wave and S-wave Poynting vectors, equation (10) can be equivalently written as:

[0087] in It is the P-wave Poynting vector of the source wave field. It is the P-wave Poynting vector of the received wave field. The S-wave Poynting vector of the received wave field can be calculated using equation (6).

[0088] On the other hand, the S-wave decoupled by the curl operator exhibits polarity reversal, which leads to inconsistent polarity of the imaging results at different locations when superimposing different single-shot PS imaging results, resulting in energy cancellation during the superposition process. Therefore, polarity correction of the S-wave is required during imaging. A commonly used polarity correction method in elastic wave RTM is to calculate the directional relationship between the incident and reflected wave fields using the Poynting vector, and then correct the S-wave polarity. For PS wave RTM results, the polarity correction method can be expressed as:

[0089] in This is the PS migration result after S-wave correction.

[0090] Figure 1 This is a flowchart of the method for decoupling the Poynting vector based on divergence and curl operators and its application in the reverse time migration of elastic waves, provided by an embodiment of the present invention.

[0091] Figure 2 This is a schematic diagram of a double-layer graben model provided in an embodiment of the present invention. Figure 2 Part (a) in the figure is used for the longitudinal and transverse wave velocities of the real medium model in forward modeling; Figure 2 Part (b) is used for the longitudinal and transverse wave velocities of the smooth medium model for wavefield inverse time delay topology and imaging; Figure 2 Part (c) describes the observation system for the five shot points used in this test.

[0092] Figure 3 This is a snapshot of the wavefield before and after decoupling the wavefield at the third shot point in the double-layer graben model provided in this embodiment of the invention. The wavefield energy has been normalized, and the P and S waves have been completely decoupled. Figure 3 (a) A snapshot of the vx component of the wavefield before decoupling in part; Figure 3 (b) is a snapshot of the vz component of the decoupled wavefield. Figure 3 A snapshot of the P-wave field after decoupling in part (c); Figure 3 A snapshot of the S-wave field after decoupling part (d) in the image.

[0093] Figure 4 This is a snapshot of the Poynting vector of the third shot point in the double-graben model provided in this embodiment of the invention before and after decoupling. After decoupling, the Poynting vectors of the P-wave and S-wave are completely separated, and each accurately represents the propagation direction of the corresponding wave. Figure 4 (a) The x-component of the Poynting vector before decoupling; Figure 4 (b) The x-component of the Poynting vector of the P-wave after decoupling; Figure 4 (c) The x-component of the Poynting vector of the S-wave after decoupling; Figure 4(d) The z-component of the Poynting vector before decoupling; Figure 4 (e) The z-component of the Poynting vector of the P-wave after decoupling; Figure 4 (f) is the z-component of the decoupled S-wave Poynting vector.

[0094] Figure 5 The image shows the single-shot RTM imaging result of the third shot point in the double-layer graben model provided in this embodiment of the invention. Without Poynting vector or traditional method constraints, obvious imaging artifacts are present at the red circles and arrows. After applying Poynting vector constraints, the imaging artifacts are significantly suppressed. Figure 5 (a) PP wave imaging results without using Poynting vector constraints; Figure 5 (b) shows the PP wave imaging results using traditional Poynting vector constraints; Figure 5 (c) shows the PP wave imaging results using the decoupled P-wave and S-wave Poynting vector constraints of this invention; Figure 5 (d) PS wave imaging results without using Poynting vector constraints; Figure 5 (e) shows PS-wave imaging results using traditional Poynting vector constraints; Figure 5 (f) shows the PS-wave imaging results using the decoupled P-wave and S-wave Poynting vector constraints of this invention.

[0095] Figure 6 The above is the RTM stacking result of five shot points in the double-layer graben model provided in this embodiment of the invention. After applying the Poynting vector constraint proposed in this invention, the polarity of the S-wave is accurately corrected, the offset noise is suppressed, and the imaging quality of the stacked offset profile is significantly improved. Figure 6 (a) PP wave imaging results without using Poynting vector constraints; Figure 6 (b) shows the PP wave imaging results using traditional Poynting vector constraints; Figure 6 (c) shows the PP wave imaging results using the decoupled P-wave and S-wave Poynting vector constraints of this invention; Figure 6 (d) PS wave imaging results without using Poynting vector constraints; Figure 6 (e) shows PS-wave imaging results using traditional Poynting vector constraints; Figure 6 (f) shows the PS-wave imaging results using the decoupled P-wave and S-wave Poynting vector constraints of this invention.

[0096] Figure 7 This is a schematic diagram of the Marmousi model provided in an embodiment of the present invention. Figure 7(a) in the figure is the accurate P-wave velocity model used in numerical simulation; Figure 7 (b) is used for the smooth longitudinal wave velocity model of wavefield inverse time extension and RTM.

[0097] Figure 8 The image shows the PP wave migration imaging results of the Marmousi model provided in this embodiment of the invention. The red circle indicates that the method of this invention displays a clearer imaging interface in the steep tilt region. Figure 8 (a) shows the superimposed imaging results using traditional Poynting vector constraints; Figure 8 (b) shows the superimposed imaging results using the decoupled P-wave and S-wave Poynting vector constraints of the present invention.

[0098] Figure 9 The image shows the PS-wave migration imaging results of the Marmousi model provided in this embodiment of the invention. The PS-wave migration results with decoupled Poynting vector constraints proposed in this invention show significantly clearer imaging of shallow reflective interfaces (indicated by the red circle). Figure 9 (a) shows the superimposed imaging results using traditional Poynting vector constraints; Figure 9 (b) shows the superimposed imaging results using the decoupled P-wave and S-wave Poynting vector constraints of the present invention.

[0099] Figure 10 The present invention provides a comparison of the computation time of the proposed method with that of the traditional method for RTM imaging. The computation time of the proposed method is approximately 92.92% of that of the improved traditional decoupling method, which is more efficient.

[0100] Based on the decoupled Poynting vector calculation method proposed in this invention and its application in elastic wave reverse time migration, it can be widely applied to the following fields or related products.

[0101] Oil and gas exploration and development: As a core application area, this invention can be used to process seismic data collected on land and at sea, and to perform high-precision imaging of complex geological structures (such as faults, salt domes, small-scale traps, steeply dipping strata, etc.), effectively improving the accuracy of reservoir prediction and oil and gas reservoir description.

[0102] Mineral resource exploration: It can be applied to the detailed exploration of solid minerals such as metallic and non-metallic minerals. Through clear structural imaging, it helps to determine the location, occurrence and boundaries of ore bodies.

[0103] Geothermal resource development: Precise imaging of underground thermal reservoir structures helps identify drilling targets and assess the development potential of geothermal resources.

[0104] Engineering geology and disaster assessment: It can be used for site stability assessment of major projects (such as tunnels, dams, and nuclear power plant site selection), as well as for internal structure detection and monitoring of potential geological hazards such as landslides and ground subsidence.

[0105] Basic research in Earth sciences: can provide high-quality geophysical imaging data for studying scientific issues such as the fine structure of the Earth's crust and the activity of fault zones.

[0106] Related products: The method of this invention can be integrated into commercial or open-source geophysical data processing and interpretation software platforms as the core algorithm of their elastic wave reverse time migration imaging module, providing users with higher quality imaging services.

[0107] The present invention has been numerically simulated and tested using a standard double graben model and a complex Marmousi model, and directly compared with existing technologies, obtaining sufficient evidence of technical effectiveness and demonstrating the significant progress of the present invention.

[0108] 1. It has an excellent effect on suppressing imaging artifacts and low-frequency noise.

[0109] like Figure 5 The comparison of single-shot imaging results for the double-layer graben model shown shows the method without using Poynting vector constraints. Figure 5 (a) in the middle, Figure 5 (d) contains numerous arcuate artifacts caused by multiple waves and gyratory waves (circled by red dashed lines). Using the traditional hybrid Poynting vector ( Figure 5 (b) in the middle Figure 5 While (e) shows improvement, residual noise remains. However, after applying the decoupled Poynting vector constraint proposed in this invention (…),… Figure 5 (c) in the middle Figure 5 In (f) of the model, whether it is PP wave or PS wave imaging, imaging artifacts are significantly suppressed, the profile is clean, and the signal-to-noise ratio is greatly improved.

[0110] 2. Accurate polarity correction of the S-wave was successfully achieved, greatly improving the quality of the superimposed imaging.

[0111] like Figure 6 The multi-shot stacking results shown are in the uncorrected ( Figure 6 (d) or correction using traditional methods ( Figure 6 In the PS imaging profile (e) of the graben, due to polarity issues, the energy of the reflective interface at the bottom of the graben is weak and blurry. However, after polarity correction using the method of this invention ( Figure 6 In (f)), the S-wave energy is correctly superimposed, the graben interface is clear, continuous and has strong energy, and the imaging quality has achieved a qualitative leap. This directly proves that the present invention has solved the problem of S-wave polarity correction.

[0112] 3. It has a stronger imaging capability for complex geological structures.

[0113] In internationally recognized tests of complex Marmousi models, such as Figure 8 and Figure 9 As shown, comparing the imaging results of traditional methods with those of the present invention, in the steeply dipping faults and shallow complex structural areas marked by red dashed circles, the stratigraphic interfaces obtained by the present invention are clearer, more continuous, and more accurately focused, demonstrating its superiority in dealing with complex geological problems.

[0114] 4. It has higher computational efficiency.

[0115] like Figure 10 The computation time comparison shown indicates that, when completing the same Marmousi model imaging task, the method of this invention takes 4814.04 seconds, which is approximately 7% more efficient than an improved traditional decoupling method (5180.76 seconds). This demonstrates that the present invention significantly improves imaging quality while also possessing the advantage of high computational efficiency.

[0116] In summary, a large amount of model test data irrefutably demonstrates the significant advantages and beneficial effects of this invention compared to existing technologies from multiple dimensions, including noise suppression, S-wave polarity correction, complex structure imaging, and computational efficiency.

[0117] It should be noted that embodiments of the present invention can be implemented in hardware, software, or a combination of both. The hardware portion can be implemented using dedicated logic; the software portion can be stored in memory and executed by a suitable instruction execution system, such as a microprocessor or dedicated-design hardware. Those skilled in the art will understand that the above-described devices and methods can be implemented using computer-executable instructions and / or included in processor control code, for example, such code provided on a carrier medium such as a disk, CD, or DVD-ROM, a programmable memory such as read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The devices and modules of the present invention can be implemented by hardware circuitry such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field-programmable gate arrays, programmable logic devices, etc., or by software executed by various types of processors, or by a combination of the above-described hardware circuitry and software, such as firmware.

[0118] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions, and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention, and within the spirit and principles of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for calculating the Poynting vector based on the decoupling of the divergence operator and the curl operator in an elastic wave field, characterized in that, Includes the following steps: a) Time-domain extension of the source wavefield or receiver wavefield using the first-order velocity-stress elastic wave equation; b) Perform divergence calculation on the obtained velocity wave field to obtain the longitudinal wave component, and perform curl calculation on the velocity wave field to obtain the transverse wave component; c) The volumetric strain is obtained by summing the normal stress components and the gradient of the strain is obtained to obtain the longitudinal wave acceleration field. The transverse wave acceleration field is obtained by subtracting the longitudinal wave acceleration field from the overall acceleration field. d) Construct the longitudinal wave Poynting vector using the longitudinal wave acceleration wave field and longitudinal wave components, and construct the transverse wave Poynting vector using the transverse wave acceleration wave field and transverse wave components.

2. The method according to claim 1, characterized in that, The volumetric strain is obtained by multiplying the sum of the three normal stress components by a scaling factor γ calculated based on the first Lamé constant and the second Lamé constant.

3. The method according to claim 1, characterized in that, Divergence and curl operations use the same finite difference operators as the wave equation discretization and reuse the already calculated spatial gradient data, thereby reducing the computational cost of decoupling.

4. A method for elastic wave reverse-time migration imaging based on the Poynting vector according to any one of claims 1 to 3, characterized in that, include: a) Perform forward extension of the source wavefield and reverse extension of the received wavefield; b) At each moment, the angle between the propagation directions of the source wavefield and the received wavefield is determined using the Poynting vectors of the P-wave and the S-wave. If the dot product of the P-wave Poynting vector is negative, P-wave-P-wave normalized cross-correlation imaging is performed. If the dot product of the P-wave Poynting vector and the S-wave Poynting vector is negative, P-wave-S-wave normalized cross-correlation imaging is performed. Otherwise, imaging energy is not accumulated at the corresponding position. c) The polarity of the P-wave-S-wave imaging results is corrected according to the sign of the cross product of the P-wave Poynting vector and the S-wave Poynting vector.

5. The imaging method according to claim 4, characterized in that, The normalized cross-correlation numerator is the sum of the products of the corresponding components of the source wavefield and the received wavefield at each time step, and the denominator is the sum of the squares of the source wavefield components at each time step.

6. The imaging method according to claim 4, characterized in that, The included angle threshold is set to 90 degrees, and the direction constraint is achieved without the need for inverse cosine function calculation by using the Poynting vector dot product symbol.

7. A device for calculating the reverse time migration of elastic waves, comprising a memory and a processor, characterized in that, The processor is configured to perform the method according to any one of claims 1 to 6.

8. The computing device according to claim 7, characterized in that, The processor employs a shared memory caching strategy to reuse gradient data during divergence calculation, curl calculation, and stress-strain conversion, thereby reducing the number of memory accesses.

9. A computer-readable storage medium having a computer program stored thereon, the program being executed by a processor to cause the processor to perform the method according to any one of claims 1 to 6.

10. The computer-readable storage medium according to claim 9, characterized in that, The computer program includes a wavefield extension module, a wavefield decoupling module, a Poynting vector calculation module, a direction-constrained imaging module, and a polarity correction module. Each module is called in a preset order to complete the elastic wave reverse time migration imaging.

Images