A method for fast simulation of 3D scattered wave in acoustic remote sensing in transversely isotropic media

Abstract

The application discloses a method for rapidly simulating three-dimensional scattered waves in a sound wave far detection transversely isotropic stratum, and belongs to the field of seismic wave forward simulation. The application is used to solve the technical problem that it is difficult to efficiently and universally simulate three-dimensional scattered waves in single well imaging of the transversely isotropic stratum in the prior art. The application comprises the following steps: obtaining parameters required for simulation; deriving an analytical solution of a sound source far field wave under a transversely isotropic medium according to a steepest descent integral method, and obtaining an incident wave field, wherein saddle points of quasi-SV waves and quasi-P waves are solved by using a Newton-Raphson method; establishing an elastic wave velocity-stress staggered grid, introducing the incident wave field into a calculation space boundary by using a total field / scattered field technology, and performing numerical simulation on a local region containing a scatterer to obtain a near field scattered wave field; calculating a virtual source wave field based on the analytical solution; and converting the near field scattered wave field and the virtual source wave field into a scattered wave field at a position of a well receiver based on a fluid-solid coupling reciprocity theorem. The application is suitable for a transversely isotropic medium with arbitrary parameters.

Description

Technical Field

[0001] This invention belongs to the field of seismic wave forward modeling, specifically relating to a method for rapidly simulating three-dimensional scattered waves from distant acoustic probes in transversely isotropic strata. Background Technology

[0002] Single-well imaging is an advanced technology for detecting subsurface fractures and karst reservoir structures. Compared to traditional seismic exploration, single-well imaging places both the seismic source and receiver inside the well, thus avoiding surface noise interference and enabling the detection of deeper subsurface structures. Furthermore, single-well imaging utilizes reflected waves from outside the well, extending the radial detection distance from approximately one meter around the well to tens of meters, significantly improving the detection rate of oil and gas reservoirs.

[0003] However, actual sedimentary strata typically exhibit anisotropy. Anisotropy refers to the variation of the physical properties of a stratum (such as the velocity of sound waves) with direction. Most sedimentary rocks (such as shale and sandstone) have layered structures during their formation, resulting in properties that are essentially the same in the horizontal direction but differ in the vertical direction. This specific type of anisotropy is called a vertically transversely isotropic (VTI) medium. Since actual reservoirs often exhibit characteristics of VTI media, studying the sound wave propagation mechanism in single-well imaging within this medium is of great significance for improving the detection capabilities of oil and gas reservoirs.

[0004] Previous studies have addressed wellbore acoustics in VTI formations. For example, the dispersion curves of Stoneley waves and pseudo-Rayleigh waves excited by monopole sources were obtained through computational analysis, the full acoustic field in radially layered VTI formations was analyzed, and the generation mechanisms of critically refracted P-waves and S-waves and their impact on low-frequency direct shear wave logging were investigated. However, these studies primarily focused on guided waves propagating inside the wellbore or near the wellbore wall, with relatively little research on external radiation waves in VTI formations. In single-well imaging applications, a crucial step is acquiring reflected signals from external targets, which first requires an accurate understanding of the acoustic radiation process from the wellbore to the formation. However, the coupling between quasi-P-waves and quasi-SV-waves in the VTI medium makes the calculation of wellbore radiation waves quite complex.

[0005] For calculating far-field radiation waves in VTI (Vacuum-Induced Tilt) media, some researchers have used the steepest descent method to analyze the data. By establishing the correspondence between the slow-motion surface and the wavefront, they have derived the radiation pattern of VTI formations under specific parameters. However, this method relies on a tabular analysis of specific formation parameters and is difficult to generalize to VTI formations with arbitrary parameters. Therefore, it cannot meet the need for efficient and universal calculation of the radiation wave field in single-well imaging forward modeling. Rapid and accurate calculation of far-field radiation waves is precisely the foundation for conducting single-well imaging forward modeling research.

[0006] Furthermore, in practical single-well imaging scenarios, when scatterers such as caverns and fractures exist in the formation, analytical solutions for three-dimensional scattered waves are usually unavailable, necessitating numerical simulation methods for analysis. The Finite-Difference Time-Domain (FDTD) method is a commonly used approach for solving such problems. However, this method requires fine-grid discretization of the entire VTI formation, including the borehole and far-field scatterers, resulting in enormous computational demands and making it difficult to meet practical application requirements. Therefore, how to efficiently simulate three-dimensional elastic wave scattering in single-well imaging within VTI formations has become a pressing technical challenge.

[0007] In summary, to address the issue of high computational complexity in 3D scattering wave simulation for single-well imaging in VTI media, it is urgent to develop a new method that balances computational accuracy and efficiency, thus laying a theoretical foundation for the development of remote detection technology for unconventional oil and gas reservoirs. Summary of the Invention

[0008] To address the technical problem of the difficulty in efficiently and universally simulating three-dimensional scattered waves in single-well imaging of transversely isotropic formations in existing technologies, this invention provides a method for rapidly simulating three-dimensional scattered waves from acoustic remote detection in transversely isotropic formations.

[0009] The present invention discloses a method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, comprising the following steps:

[0010] S1. Obtain the formation parameters, wellbore acoustic source parameters, receiver parameters, and scatterer parameters required for the simulation;

[0011] S2. The analytical solution of the far-field wave of the sound source under the transversely isotropic medium is derived according to the steepest descent integral method, and the incident wave field is obtained. The saddle points of the mutually coupled quasi-SV wave and quasi-P wave in the analytical solution of the far-field wave of the sound source are solved by the Newton-Raphson method, and the saddle points of the independently decoupled SH wave are solved by explicit method.

[0012] S3. Establish an elastic wave velocity-stress interleaved grid, discretize the computational space containing the target scatterer, use the total field / scattered field technique to introduce the incident wave field into the boundary of the computational space, perform numerical simulation on the local region containing the scatterer, obtain the near-field scattered wave field generated by the real source near the scatterer, and add a perfectly matched layer to the outer layer of the computational space.

[0013] S4. Based on the analytical solution, calculate the virtual source wave field generated by the virtual source near the scatterer at the receiver location;

[0014] S5. Based on the fluid-structure interaction reciprocity theorem, the near-field scattered wave field and the virtual source wave field are converted into the scattered wave field at the receiver position in the well, thus completing the simulation of three-dimensional scattered waves.

[0015] Preferably, the analytical solution of the far-field wave of the sound source in step S2 is expressed by the following far-field displacement expression:

[0016]

[0017]

[0018]

[0019] In the formula, These are the radial distance from the well axis, the vertical distance from the sound source, and the interface azimuth angle, respectively.

[0020] for directional far-field radiation displacement for directional far-field radiation displacement for Directional far-field radiation displacement;

[0021] and These are caused by quasi-P waves and quasi-SV waves, respectively. Directional displacement, and These are caused by quasi-P waves and quasi-SV waves, respectively. Directional displacement;

[0022] for Amplitude coefficient of directional far-field radiation; The radial wavenumber of the SH wave; The azimuth angle of the sound source; Angular frequency, This represents the straight-line distance from the field point to the sound source. Axial slowness, This is the saddle point of the SH wave;

[0023] It is a complex exponential function with axial slowness as the independent variable, obtained by integrating the displacement corresponding to the SH wave with the steepest descent. yes The second derivative of .

[0024] Preferably, quasi-P waves and quasi-SV waves cause Directional displacement and , caused by quasi-P waves and quasi-SV waves Directional displacement and Determined by the following formula:

[0025]

[0026]

[0027]

[0028]

[0029] In the formula, , for Amplitude coefficient of directional far-field radiation, , for Amplitude coefficient of directional far-field radiation, , These are the radial wave numbers of the quasi-P-wave and the quasi-SV-wave, respectively. The radius of the sound source; It is a unit imaginary number; As a quasi-P-wave saddle point, As a quasi-SV wave saddle point;

[0030] , These are complex exponential functions with axial slowness as the independent variable, obtained by integrating the displacements corresponding to the quasi-P-wave and quasi-SV-wave at the steepest descent, respectively. , They are respectively , The second derivative of .

[0031] Preferably, the amplitude coefficient is determined by the following formula:

[0032]

[0033]

[0034]

[0035]

[0036]

[0037] In the formula, , and These are the undetermined coefficients for quasi-P wave, quasi-SV wave, and SH wave, respectively. The order of the sound source. This represents the maximum value of the change in the volume of the sound source. The frequency function of the sound source. and These are the axial wave number and the fluid radial wave number, respectively. For quasi-P-wave coupling coefficients, The quasi-SV wave coupling coefficient.

[0038] Preferably, the undetermined coefficients of the quasi-P wave, quasi-SV wave, and SH wave are... , and The boundary conditions are determined by applying stress and displacement to the borehole wall.

[0039] Preferably, the SH wave saddle point in the analytical solution of the far-field wave of the sound source in step S2. Solve using the following explicit expression:

[0040]

[0041] In the formula, Let be the tilt angle of the radiated wave relative to the vertical axis of symmetry. The solid density of the formation medium. , The coefficients at the locations corresponding to the stiffness matrix of transversely isotropic strata. The shear elasticity coefficient is the coefficient of elasticity in the vertical plane. is the shear elasticity coefficient in the horizontal plane.

[0042] Preferably, the saddle points of the quasi-SV wave and quasi-P wave in the analytical solution of the far-field wave of the sound source in step S2 are obtained by solving the following implicit equations using the Newton-Raphson method:

[0043]

[0044] In the formula, It is a function for finding the saddle point of quasi-P-waves and quasi-SV-waves;

[0045] , The radial slowness of the quasi-P wave, , The radial slowness of the quasi-SV wave; The angle of inclination of the radiated wave relative to the vertical axis of symmetry;

[0046] Calculate using the following formula:

[0047]

[0048] In the formula, The solid density of the formation medium. , , The coefficients at the locations corresponding to the stiffness matrix of transversely isotropic strata. The elastic modulus in the horizontal direction is... The coupling elasticity coefficient, The elastic modulus is the coefficient of elasticity in the vertical direction.

[0049] Preferably, in step S3, the incident wave field is introduced into the boundary of the computational space and applied using the following discretization scheme within the elastic wave velocity-stress staggered grid framework:

[0050]

[0051]

[0052]

[0053]

[0054]

[0055]

[0056]

[0057]

[0058] In the formula, They are respectively Orientation grid index, The grid size is in the y-direction. The time step of FDTD They are respectively The velocity component in the direction; They are respectively Stress components in the direction of application;

[0059] superscript Indicates the time step index, superscript The positive and negative half-time step indices are indicated; the superscript ordinary FDTD indicates the physical quantity calculated using the traditional FDTD method, and the superscript inc indicates the incident wave field quantity calculated using the analytical expression described in step S2.

[0060] Preferably, the perfectly matched layer in step S3 is used to absorb artificially reflected waves generated at the truncated boundary.

[0061] Preferably, the fluid-structure interaction reciprocity theorem described in step S5 is expressed by the following formula:

[0062]

[0063] In the formula, This represents the fluid displacement at the receiver location inside the well. Let be the unit vector along the dipole direction. The boundary of the closed surface domain surrounding the scatterer;

[0064] and These are the displacements and stresses generated near the scatterer by the real source obtained through the numerical simulation described in step S3;

[0065] and These represent the displacement and stress generated near the scatterer by the real source obtained through the numerical simulation described in step S4;

[0066] The distance between the two mutually interchangeable points; The density of the liquid;

[0067] For a small region on the integral surface, To be perpendicular to A unit vector in direction.

[0068] The beneficial effects of this invention are:

[0069] 1. It solves the problem of insufficient versatility of traditional analytical methods.

[0070] This invention introduces the Newton-Raphson method to solve the implicit equations for saddle points of quasi-P-waves and quasi-SV-waves, replacing the traditional steepest descent method which relies on tabular analysis for specific parameter solutions. This enables the invention to quickly and accurately calculate far-field radiation waves from wellbore acoustic sources in transversely isotropic media with arbitrary parameters, overcoming the limitation of traditional methods that are only applicable to specific formation parameters, and providing a universal analytical basis for single-well imaging forward modeling.

[0071] 2. Significantly reduced the computational load of numerical simulation.

[0072] This invention employs the Total Field / Scattered Field (TF / SF) technique to isolate the local region containing the scatterer from the entire computational domain for detailed simulation. Only in this local region is an elastic wave velocity-stress interleaved grid established, and a perfectly matched layer absorbing boundary is applied. Compared to traditional finite-difference time-domain methods, which require detailed grid discretization of the entire large-scale space including the borehole and the scatterer, this invention's local simulation strategy significantly reduces memory usage and CPU consumption.

[0073] 3. It achieves an organic combination of analytical and numerical solutions.

[0074] This invention utilizes the fluid-structure interaction reciprocity theorem to perform an integral transformation between the near-field scattered wave field obtained from numerical simulation of a local region and the virtual source wave field calculated based on analytical solutions, directly obtaining the scattered wave field at the receiver location within the wellbore. This method avoids the process of the simulated wave field propagating back from the scatterer to the wellbore, further improving computational efficiency while ensuring computational accuracy.

[0075] 4. Significantly improved computational efficiency

[0076] Numerical experimental results show that for large-scale 3D single-well imaging problems under VTI backgrounds, the computational efficiency of the method proposed in this invention is significantly improved compared with the traditional global finite-difference time-domain (FDTD) method. On the same computer equipped with an i9-12900K CPU and 64GB of memory, in the example containing cubic scatterers, the CPU time of the method proposed in this invention is 4053.52 seconds, while the CPU time of the global FDTD method is 85698.52 seconds. In the large-scale example containing two ellipsoidal scatterers, the memory usage of the method proposed in this invention is approximately 3.12GB, and the CPU time is approximately 3.86 hours. In contrast, if the traditional FDTD method is used to simulate the same problem, it is estimated that it will require approximately 464.91GB of memory and 519.40 hours of CPU time.

[0077] 5. Wide range of applications

[0078] The method of this invention is applicable to transversely isotropic media with arbitrary parameters, and can simulate scatterers of arbitrary shapes (including geological anomalies such as caves and fractures). It supports multiple sound source types such as monopoles and dipoles, and can effectively simulate complex physical phenomena such as transverse wave splitting, laying a theoretical foundation for the development of remote detection technology for unconventional oil and gas reservoirs. Attached Figure Description

[0079] Figure 1 This is a flowchart of the method for simulating three-dimensional scattered wave imaging in a transversely isotropic medium according to the present invention;

[0080] Figure 2 This is a schematic diagram of an arbitrary transversely isotropic medium model corresponding to the formula derived in this invention;

[0081] Figure 3 A schematic diagram of a single-well imaging simulation model containing a cubic scatterer provided by the present invention;

[0082] Figure 4 The figure shows a comparison between the analytical and numerical solutions used in this invention to calculate the far-field wave. The red solid line represents the result calculated by the real axis integration method (the exact numerical solution), and the black dashed line represents the result calculated by the steepest descent method used in this invention (the approximate analytical solution). Figure 4 (a), (b), and (c) show the correlation results corresponding to stratum 1, where, Figure 4 (a) is directional far-field radiation displacement Waveform comparison chart, Figure 4 (b) is directional far-field radiation displacement Waveform comparison chart, Figure 4 (c) is directional far-field radiation displacement Waveform comparison chart; Figure 4 (d), (e), and (f) are the correlation results corresponding to stratum 2, where: Figure 4 (d) is directional far-field radiation displacement Waveform comparison chart, Figure 4 (e) is directional far-field radiation displacement Waveform comparison chart, Figure 4 (f) is directional far-field radiation displacement Waveform comparison chart;

[0083] Figure 5 This is a schematic diagram of the time-domain finite-difference grid and total field / scattered field technique used in this invention, wherein, Figure 5 (a) is a schematic diagram of the elastic wave velocity-stress interleaved grid used in this invention; Figure 5 (b) is a schematic diagram of the TF / SF technology introduced in this invention on the finite-difference time-domain (FDTD) framework;

[0084] Figure 6 This figure shows a comparison of the results of calculating the far-field wave of a borehole acoustic source using the hybrid method (analytical solution + finite-difference time domain) of this invention and the real-axis integration method. The red solid line represents the result calculated using the real-axis integration method (the accurate numerical solution), and the black dashed line represents the result calculated using the hybrid method (the method of this invention). Figure 6 (a) is directional far-field radiation displacement Waveform comparison chart, Figure 6 (b) is directional far-field radiation displacement Waveform comparison Figure 6 (c) is directional far-field radiation displacement Waveform comparison;

[0085] Figure 7 This is a schematic diagram illustrating the reciprocity relationship of the single-well imaging model used in this invention to process in-well wave fields, wherein... Figure 7 (a) is a schematic diagram of the actual source wave field. Figure 7 (b) is a schematic diagram of the virtual source wave field;

[0086] Figure 8This is a comparison of the results of the proposed method and the traditional method under a cubic scattering model. The black solid line represents the calculation result of the reference method, and the red dashed line represents the calculation result of the hybrid method.

[0087] Figure 9 A schematic diagram of a large-scale model containing two far-field scatterers provided by the present invention;

[0088] Figure 10 The figure shows the result of calculating a large-scale model containing two far-field scatterers using the method proposed in this invention. Figure 10 (a) is a scattering waveform recorded by six receivers. Figure 10 (b) is a waveform diagram of the scattered wave velocity;

[0089] Figure 11 A schematic diagram of a scattering model for simulating shear wave splitting provided by the present invention;

[0090] Figure 12 This is a simulation result of shear wave splitting according to the present invention, wherein, Figure 12 (a) shows the simulation results received by the six receivers in this transverse wave splitting simulation. Figure 12 (b) shows the reception results of receiver 1 in the transverse wave splitting simulation;

[0091] Figure 13 This is a slowness surface diagram of three types of waves in a transversely isotropic background medium under the transverse wave splitting simulation of this invention. Detailed Implementation

[0092] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0093] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0094] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, but this is not intended to limit the scope of the invention.

[0095] Specific Implementation Method 1: The following is combined with... Figures 1 to 13 This embodiment describes a method for rapidly simulating three-dimensional scattered acoustic waves at a distance in transversely isotropic strata, comprising the following steps:

[0096] Step 1: Parameter Preprocessing

[0097] First, based on the simulated formation conditions, input various parameters for this simulation example, including formation medium information, wellbore sound source parameters, receiver parameters, and the location and geometric dimensions of the scatterer.

[0098] In this specific embodiment, with Figure 3 The single-well imaging model containing cubic scatterers shown is used as an example for illustration. The monopole sound source is located at (0,0,0.5) m, with a center frequency of 3 kHz and a half-bandwidth of 2 kHz. Seven receivers are evenly distributed between (0,0,2.5) m and (0,0,3.7) m. VTI formation medium information is shown in Table 1, where formation 1 and formation 2 correspond to different combinations of elastic parameters.

[0099] Table 1 Formation Parameters

[0100]

[0101] , , , , The coefficients at the locations corresponding to the stiffness matrix of transversely isotropic strata. The elastic modulus in the horizontal direction is... The coupling elasticity coefficient, The elastic modulus is in the vertical direction. The shear elasticity coefficient is the coefficient of elasticity in the vertical plane. is the shear elasticity coefficient in the horizontal plane.

[0102] Step 2: Derive the analytical solution for far-field waves

[0103] This step aims to provide the necessary incident wave field expression for the subsequent total field / scattered field boundary, and is specifically divided into the following four parts:

[0104] Step 2.1: Expression of the potential function

[0105] In cylindrical coordinates, the elastic wave field excited by a sound source in a VTI medium can be represented by three potential functions: quasi-P-wave potential function. Quasi-SV wave potential function and SH wave potential function Their expressions are as follows:

[0106] (1)

[0107] (2)

[0108] (3)

[0109] In the formula, These are the radial distance from the well axis, the vertical distance from the sound source, and the interface azimuth angle, respectively. Angular frequency, and These are the axial wave number and the fluid radial wave number, respectively. and These are the sound source radius and azimuth angle, respectively. This represents the maximum value of the change in the volume of the sound source. The frequency function of the sound source. Category II The modified Bessel function of order 1 , , These are the radial wavenumbers of the SH wave, quasi-P wave, and quasi-SV wave, respectively. It is a unit imaginary number. It is the natural logarithm.

[0110] , and These are the undetermined coefficients for quasi-P-waves, quasi-SV-waves, and SH-waves, determined by applying stress and displacement boundary conditions on the wellbore wall. Specifically, at the wellbore wall, the conditions of continuous fluid radial displacement and solid radial displacement, continuous fluid pressure and solid radial stress, and zero solid shear stress must be met, thereby obtaining the coefficients for... , and The three undetermined coefficients can be determined by solving the system of linear equations.

[0111] It is worth noting that, according to equation (1-3), quasi-P waves and quasi-SV waves are coupled to each other in VTI formations, while SH waves are independently decoupled. For quasi-P-wave coupling coefficients, The quasi-SV wave coupling coefficient reflects the mutual coupling relationship between the quasi-P wave and the quasi-SV wave.

[0112] The radial wavenumber and intermediate parameters are determined by the following formula:

[0113] (4a)

[0114] , (4b)

[0115] (4c)

[0116] (4d)

[0117] (4e)

[0118] (4f)

[0119] (4g)

[0120] In the formula, and These are the fluid density and the solid density, respectively. , , These are intermediate parameters introduced for the convenience of computational expression.

[0121] Substituting equation (1-3) into the Helmholtz decomposition, we obtain the frequency domain shift expression:

[0122] (5)

[0123] (6)

[0124] (7)

[0125] In the formula, for directional far-field radiation displacement for directional far-field radiation displacement for Directional far-field radiation displacement; for Amplitude coefficient of directional far-field radiation; The second The modified Bessel function of order 1

[0126] , for Amplitude coefficient of directional far-field radiation, , for The amplitude coefficients of the directional far-field radiation are as follows:

[0127] (8a)

[0128] (8b)

[0129] (8c)

[0130] (8d)

[0131] (8f)

[0132] In the formula, n is the order of the sound source, n=0 for a monopole source, n=1 for a dipole source, and n=2 for a quadrupole source.

[0133] Step 2.2: Calculate the far-field displacement using the steepest descent integral method

[0134] The real-axis integration method is a commonly used direct method for calculating radiation displacement (Equation 5-7). However, obtaining the spatial domain wave field through wavenumber integration is computationally expensive. When the radial distance is much larger than the dominant wavelength of the radiation wave, the following conditions must be met: , and Under the condition of correcting the far-field asymptote of the Bessel function for:

[0135] (9)

[0136] Substituting equation (9) into equation (5-7) and using the steepest descent integral method, the far-field displacement can be derived:

[0137] (10)

[0138] (11)

[0139] (12)

[0140] in

[0141] (13a)

[0142] (13b)

[0143] (13c)

[0144] (13d)

[0145] The complex exponential function is defined as:

[0146] (14a)

[0147] (14b)

[0148] (14c)

[0149] Axial slowness and radial slowness satisfy:

[0150] , , , , (15)

[0151] In the formula, This represents the straight-line distance from the field point to the sound source. Axial slowness, , and These are the radial slowness of the quasi-P wave, quasi-SV wave, and SH wave, respectively. , and These are the axial slow saddle points for quasi-P-wave, quasi-SV-wave, and SH-wave, respectively. and These are caused by quasi-P waves and quasi-SV waves, respectively. Directional displacement, and These are caused by quasi-P waves and quasi-SV waves, respectively. Directional displacement, , and These are complex exponential functions with axial slowness as the independent variable, obtained by integrating the displacements corresponding to the quasi-P wave, quasi-SV wave, and SH wave at the steepest descent. , They are respectively , The second derivative, The angle of inclination of the radiated wave relative to the vertical axis of symmetry.

[0152] Step 2.3: Solving for saddle points

[0153] For the azimuth displacement component composed of pure SH waves, the following equation can be obtained: Explicit expression:

[0154] (16)

[0155] for The first derivative.

[0156] Substituting equations (4a, 14c, 15) into equation (16), we can obtain the saddle point corresponding to SH:

[0157] (17)

[0158] For the z-component and r-component displacements composed of quasi-P-waves and quasi-SV-waves, the saddle point and Determined by the following formula ( The saddle point corresponding to the quasi-P wave, (Saddle point corresponding to quasi-SV wave)

[0159] (18)

[0160] for The first derivative.

[0161] Substituting equations (4b-4e, 14a-14b, 15) into equation (18), we obtain the implicit conditions for saddle points:

[0162] (19)

[0163] It is a quasi-P wave ( ) and quasi-SV waves ( Find the function needed to find the saddle point.

[0164] in:

[0165] (20a)

[0166] (20b)

[0167] (20c)

[0168] (20d)

[0169] in , , Intermediate parameters introduced into the calculation process.

[0170] in Corresponding to the quasi-P wave, Corresponding to the quasi-SV wave. Since equation (19) cannot be solved explicitly, that is, the implicit nature of equations (20a-20d) makes it impossible to solve directly. The root of the problem. Based on the slow surface equation, previous researchers used the tabular method to solve for the saddle point of VTI medium. However, this method is only applicable to specific parameters and lacks universality. To obtain the saddle point of any VTI formation, this paper uses the Newton-Raphson method to solve equation (19), and the detailed process is shown in Table 2. This method makes the analytical solution applicable to VTI formations with arbitrary parameters, thus overcoming the limitation of traditional methods that rely on tabular analysis.

[0171] Table 2 Newton-Raphson method framework for determining the slowness saddle points of quasi-P-waves and quasi-SV-waves in VTI formations

[0172]

[0173] Step 2.4: Far-field stress and coordinate system transformation

[0174] After obtaining the saddle points corresponding to the three waves, the constitutive relation of the VTI medium is combined and neglected. From this, we can obtain the expression for the far-field stress:

[0175] (21a)

[0176] (21b)

[0177] (21c)

[0178] (21d)

[0179] (21e)

[0180] (21f)

[0181] In the formula, , and These represent the normal stresses in the directions corresponding to the subscripts. , and The subscript represents the shear stress in the direction of application.

[0182] To incorporate borehole radiation waves into the FDTD algorithm, the physical quantity expressions need to be transformed to Cartesian coordinates using the following formula. :

[0183] (22a)

[0184] (22b)

[0185] In the formula, , and This represents the displacement in the direction corresponding to the subscript. , This represents the normal stress in the direction corresponding to the subscript. , and This represents the shear stress in the direction corresponding to the subscript.

[0186] Substituting equations (10-12) and (21a-21f) into equation (22a-22b) and using inverse Fourier transform, we can obtain the time-domain displacement and stress of the VTI medium generated by a borehole sound source of any order, i.e., the required incident wave field.

[0187] To verify the accuracy of the derived analytical solution, the results calculated using the steepest descent method are compared with the numerical solution calculated using the real axis integration method, such as... Figure 4 As shown, the results from both methods agree well in all six subplots, verifying the accuracy of the steepest descent method. Approximate solutions can be used instead of numerical solutions in wavefield calculations.

[0188] Step 3: Establish a local FDTD simulation framework based on TF / SF technology

[0189] Step 3.1: Mesh Generation

[0190] Using an elastic wave velocity-stress interleaved mesh, such as Figure 5 As shown in (a), the velocity field is defined at half-integer points of the grid, and the stress field is defined at integer points of the grid. They are updated alternately over time to ensure second-order accuracy.

[0191] Step 3.2: TF / SF Region Partitioning

[0192] To reduce the simulation space and thus the computational load, TF / SF technology is introduced into the FDTD framework, such as... Figure 5 As shown in (b), the entire simulation domain is divided into two regions:

[0193] Region 1 (Total Field Region): Contains scatterers, calculates the total field (incident field + scattered field).

[0194] Region 2 (scattering field region): does not contain scatterers, only the scattering field is calculated.

[0195] Step 3.3: Apply boundary incident wave

[0196] An analytical incident wave calculated in step 2 is applied to the connecting boundary between regions 1 and 2 to ensure the continuity of the velocity and stress fields. The discretization scheme for the left boundary is as follows:

[0197] (23a)

[0198] (23b)

[0199] (23c)

[0200] (23d)

[0201] (23e)

[0202] (23f)

[0203] (23g)

[0204] (23h)

[0205] In the formula, They are respectively Orientation grid index, The grid size is in the y-direction. The time step of FDTD They are respectively The velocity component in the direction; They are respectively Stress components in the direction of application;

[0206] superscript Indicates the time step index, superscript The positive and negative half-time step indices are indicated; the superscript "ordinary FDTD" indicates the physical quantity calculated using the traditional FDTD method, and the superscript "inc" indicates the incident wave field quantity calculated using the analytical expression described in step S2 (obtained through equations 10-13 and 21-22). The equations for calculating the physical quantities of the remaining surfaces of the incident boundary can be derived using a method similar to that used in equations (23a-23h).

[0207] Step 3.4: Absorption Boundary

[0208] A perfect matching layer (PML) is placed on the outer layer of the entire simulation area to absorb artificial reflections caused by the truncation boundary and improve simulation accuracy.

[0209] Step 4: Method Validation

[0210] Before the formal simulation, the accuracy of the established FDTD algorithm needs to be verified. In VTI formations where scatterers are not considered, numerical and theoretical solutions are calculated and compared simultaneously, such as... Figure 6 As shown, the comparison results of the three displacement components show that the waveforms of the hybrid method (black dashed line) and the real axis integration method (red solid line) are in high agreement, verifying the accuracy of the hybrid method in calculating far-field radiation waves. If discrepancies occur, the difference scheme needs to be optimized or the mesh scale reduced to improve simulation accuracy.

[0211] The results of the FDTD simulation space established in this example are accurate. If the source wave field simulated by the hybrid method does not match the real axis integration method very well, it is necessary to optimize the difference and improve the simulation accuracy by reducing the grid scale to ensure the accuracy of the hybrid method.

[0212] Step 5: Wavefield Calculation

[0213] Step 5.1: Calculation of the real source wave field

[0214] The required scatterer is added to the established TF / SF-FDTD framework. The near-field scattered wave field generated by the real source near the scatterer is calculated using the numerical method obtained in step 3, including displacement. and stress .

[0215] Step 5.2: Virtual Source Wavefield Calculation

[0216] A virtual monopole sound source is placed at each receiver location. The virtual source wave field generated near the scatterer by the virtual source is calculated using the analytical solution of the VTI medium sound source obtained in step 2, including displacement. and stress .

[0217] Step 6: Wavefield transformation based on reciprocity

[0218] After calculating the two wave fields, the local scattered wave is converted into a scattered wave field at the receiver location within the well based on the fluid-structure interaction reciprocity theorem. For the single-well imaging model, states A and B represent the real source wave field and the virtual source wave field, respectively, as follows: Figure 7 As shown.

[0219] Reciprocity is expressed by the following formula:

[0220]

[0221] In the formula, This represents the fluid displacement at the receiver location inside the well. Let be the unit vector along the dipole direction. The boundary of the closed surface domain surrounding the scatterer;

[0222] and These are the displacements and stresses generated near the scatterer by the real source obtained through the numerical simulation described in step S3;

[0223] and These represent the displacement and stress generated near the scatterer by the real source obtained through the numerical simulation described in step S4;

[0224] The distance between the two mutually interchangeable points; Liquid density

[0225] For a small region on the integral surface, To be perpendicular to A unit vector in direction.

[0226] By using surface integral, the local FDTD calculation results are directly converted into the scattered wave field at the receiver location in the well, avoiding the process of simulated waves returning from the scatterer to the wellbore, and further improving the calculation efficiency.

[0227] Step 7: Results Output and Analysis

[0228] After completing the above steps, extract the calculated data and plot the waveform. Figure 8The y-component velocity waveforms obtained by the hybrid method (red dashed line) and the reference method (black solid line) were compared. The results of the two methods are in good agreement, which verifies the accuracy of the method proposed in this invention.

[0229] On the same computer equipped with an i9-12900K CPU and 64GB of memory, the CPU time of the hybrid method in this example is 4053.52 seconds, while the CPU time of the global FDTD method is 85698.52 seconds, indicating that the method of the present invention has extremely high computational efficiency compared with the traditional FDTD method.

[0230] The method proposed in this invention can perform efficient simulations of any VTI medium (including any scatterer). To demonstrate the versatility of this invention, an application example for a large-scale model containing two far-field scatterers is provided below:

[0231] Application Example 1: Large-scale model with two far-field scatterers

[0232] To demonstrate the versatility of this invention, Figure 9 The simulation is performed using a large-scale model with two far-field scatterers as an example. The entire computational domain measures 60m × 60m × 60m. The three-axis lengths of the two ellipsoids are (a = 2m, b = 1m, c = 1m), and their centers are located at (x = 30m, y = 0m, z = 0m) and (x = 16m, y = 0m, z = 30m), respectively. The x-polarized dipole source is located at (z = 0m, r = a = 0.1m, Six receivers were evenly distributed within the range of z=2m to z=12m. The background strata and ellipsoidal strata used the parameters of strata 1 and strata 2, respectively.

[0233] Figure 10 The result of calculating the model using the method of the present invention, wherein Figure 10 (a) shows the scattering waveforms recorded by the six receivers. Figure 10 (b) shows the velocity waveform of the scattered wave. The interaction between the scattered wave and multiple curved scatterers results in a greater number and more complex wave groups, which are significantly different from the wave groups generated by a single cubic scatterer. Meanwhile, the hybrid method requires approximately 3.12 GB of memory and 3.86 hours of CPU time; while the traditional FDTD method is expected to require approximately 464.91 GB of memory and 519.40 hours of CPU time to simulate the same problem. This example demonstrates that the hybrid method can efficiently solve large-scale single-well imaging model problems with multiple scatterers in a VTI background.

[0234] Application Example 2: Single-well imaging model with shear wave splitting phenomenon

[0235] To verify the effectiveness of this method in simulating shear wave splitting, we used... Figure 11The model shown is used as an example for simulation. The scatterer is a 1m×1m×1m cuboid rotated 45° along the z-axis, located at (x=13m, y=0m, z=5m). The borehole axis coincides with the z-axis, and the borehole radius is 0.1m. The x-polarized dipole sound source is located at (z=0m, r=a=0.1m, The center frequency is 3kHz, the half-bandwidth is 2kHz, and six receivers are evenly distributed within the range of z=2m to z=12m. The total computational domain size is 30m×30m×30m.

[0236] Traditional FDTD methods require significant computational resources. The hybrid method of this invention requires only approximately 0.46 GB of memory and 1.01 hours of CPU time to obtain all waveforms recorded by six receivers, such as... Figure 12 As shown in (a). Figure 12 (b) The x-component velocity waveform recorded by receiver 1 is shown, and four distinct wave groups can be clearly observed: quasi-PP wave, quasi-P-SV / SV-P wave, quasi-SV-SV wave and SH-SH wave, indicating that shear wave splitting still occurs in VTI formations even in the presence of a single target body.

[0237] To explain this phenomenon, the slowness vectors of three wave modes in the VTI background medium were calculated by solving the Kelvin-Christopher equation, and their slowness surfaces were displayed, as shown below. Figure 13 As shown in the figure. The results revealed that the quasi-SV wave velocity, which propagates at an angle relative to the borehole axis, differs from the SH wave velocity, leading to shear wave splitting. However, in isotropic formations, the SV and SH waves have the same velocity in any direction, thus this phenomenon does not occur. Furthermore, compared to other types of scattered waves, the SH-SH wave has the largest amplitude, indicating that the SH wave generated by the dipole source has an advantage in far-field detection of VTI formations.

[0238] While the invention has been described herein with reference to specific embodiments, it should be understood that these embodiments are merely examples of the principles and applications of the invention. Therefore, it should be understood that many modifications can be made to the exemplary embodiments, and other arrangements can be designed without departing from the spirit and scope of the invention as defined by the appended claims. It should be understood that different dependent claims and features described herein can be combined in ways different from those described in the original claims. It is also understood that features described in conjunction with individual embodiments can be used in other described embodiments.

Claims

1. A method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, characterized in that, Includes the following steps: S1. Obtain the formation parameters, wellbore acoustic source parameters, receiver parameters, and scatterer parameters required for the simulation; S2. The analytical solution of the far-field wave of the sound source under the transversely isotropic medium is derived according to the steepest descent integral method, and the incident wave field is obtained. The saddle points of the mutually coupled quasi-SV wave and quasi-P wave in the analytical solution of the far-field wave of the sound source are solved by the Newton-Raphson method, and the saddle points of the independently decoupled SH wave are solved by explicit method. S3. Establish an elastic wave velocity-stress interleaved grid, discretize the computational space containing the target scatterer, use the total field / scattered field technique to introduce the incident wave field into the boundary of the computational space, perform numerical simulation on the local region containing the scatterer, obtain the near-field scattered wave field generated by the real source near the scatterer, and add a perfectly matched layer to the outer layer of the computational space. S4. Based on the analytical solution, calculate the virtual source wave field generated by the virtual source near the scatterer at the receiver location; S5. Based on the fluid-structure interaction reciprocity theorem, the near-field scattered wave field and the virtual source wave field are converted into the scattered wave field at the receiver position in the well, thus completing the simulation of three-dimensional scattered waves.

2. The method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, as described in claim 1, is characterized in that... The analytical solution of the far-field wave of the sound source described in step S2 is expressed by the following far-field displacement expression: In the formula, These are the radial distance from the well axis, the vertical distance from the sound source, and the interface azimuth angle, respectively. for directional far-field radiation displacement for directional far-field radiation displacement for Directional far-field radiation displacement; and These are caused by quasi-P waves and quasi-SV waves, respectively. Directional displacement, and These are caused by quasi-P waves and quasi-SV waves, respectively. Directional displacement; for Amplitude coefficient of directional far-field radiation; The radial wavenumber of the SH wave; The azimuth angle of the sound source; Angular frequency, This represents the straight-line distance from the field point to the sound source. Axial slowness, This is the saddle point of the SH wave; It is a complex exponential function with axial slowness as the independent variable, obtained by integrating the displacement corresponding to the SH wave with the steepest descent. yes The second derivative of .

3. The method for rapidly simulating three-dimensional scattering waves from transversely isotropic strata using acoustic long-range detection, as described in claim 2, is characterized in that... Quasi-P waves and quasi-SV waves caused Directional displacement and , caused by quasi-P waves and quasi-SV waves Directional displacement and Determined by the following formula: In the formula, , for Amplitude coefficient of directional far-field radiation, , for Amplitude coefficient of directional far-field radiation, , These are the radial wave numbers of the quasi-P-wave and the quasi-SV-wave, respectively. The radius of the sound source; It is a unit imaginary number; As a quasi-P-wave saddle point, As a quasi-SV wave saddle point; , These are complex exponential functions with axial slowness as the independent variable, obtained by integrating the displacements corresponding to the quasi-P-wave and quasi-SV-wave at the steepest descent, respectively. , They are respectively , The second derivative of .

4. The method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, as described in claim 3, is characterized in that... The amplitude coefficient is determined by the following formula: In the formula, , and These are the undetermined coefficients for quasi-P wave, quasi-SV wave, and SH wave, respectively. The order of the sound source. This represents the maximum value of the change in the volume of the sound source. The frequency function of the sound source. and These are the axial wave number and the fluid radial wave number, respectively. For quasi-P-wave coupling coefficients, The quasi-SV wave coupling coefficient.

5. The method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, as described in claim 4, is characterized in that... The undetermined coefficients of the quasi-P wave, quasi-SV wave, and SH wave , and The boundary conditions are determined by applying stress and displacement to the borehole wall.

6. The method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, as described in claim 2, is characterized in that... The SH wave saddle point in the analytical solution of the far-field wave of the sound source described in step S2 Solve using the following explicit expression: In the formula, Let be the tilt angle of the radiated wave relative to the vertical axis of symmetry. The density of the solid medium in the formation. , The coefficients at the locations corresponding to the stiffness matrix of transversely isotropic strata. The shear elasticity coefficient is the coefficient of elasticity in the vertical plane. is the shear elasticity coefficient in the horizontal plane.

7. The method for rapidly simulating three-dimensional scattering waves from transversely isotropic strata using acoustic long-range detection, as described in claim 2, is characterized in that... The saddle points of the quasi-SV wave and quasi-P wave in the analytical solution of the far-field wave of the sound source described in step S2 are obtained by solving the following implicit equations using the Newton-Raphson method: In the formula, It is a function for finding the saddle point of quasi-P-waves and quasi-SV-waves; , The radial slowness of the quasi-P wave, , The radial slowness of the quasi-SV wave; The angle of inclination of the radiated wave relative to the vertical axis of symmetry; Calculate using the following formula: In the formula, The density of the solid medium in the formation. , , The coefficients at the locations corresponding to the stiffness matrix of transversely isotropic strata. The elastic modulus in the horizontal direction is... The coupling elasticity coefficient, The elastic modulus is the coefficient of elasticity in the vertical direction.

8. The method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, as described in claim 7, is characterized in that... In step S3, the incident wave field is introduced into the boundary of the computational space and applied using the following discretization scheme within the elastic wave velocity-stress staggered mesh framework: In the formula, They are respectively Orientation grid index, The grid size is in the y-direction. The time step of FDTD They are respectively The velocity component in the direction; They are respectively Stress components in the direction of application; superscript Indicates the time step index, superscript Indicates the positive and negative half-time step indices; The superscript "ordinary FDTD" indicates a physical quantity calculated using the traditional FDTD method, and the superscript "inc" indicates an incident wave field quantity calculated using the analytical expression described in step S2.

9. The method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, as described in claim 1, is characterized in that... The perfectly matched layer described in step S3 is used to absorb artificially reflected waves generated at the truncated boundary.

10. A method for rapidly simulating three-dimensional scattered waves from transversely isotropic strata using acoustic long-range detection, as described in claim 8, characterized in that... The fluid-structure interaction reciprocity theorem described in step S5 is expressed by the following formula: In the formula, This represents the fluid displacement at the receiver location inside the well. Let be the unit vector along the dipole direction. The boundary of the closed surface domain surrounding the scatterer; and These are the displacements and stresses generated near the scatterer by the real source obtained through the numerical simulation described in step S3; and These represent the displacement and stress generated near the scatterer by the real source obtained through the numerical simulation described in step S4; The distance between the two mutually interchangeable points; The density of the liquid; For a small region on the integral surface, To be perpendicular to A unit vector in direction.

Images