VTI medium seismic wave numerical simulation method based on pre-processing iteration

Abstract

This invention discloses a numerical simulation method for seismic waves in VTI media based on preprocessing iteration, belonging to the field of deep earth engineering geophysical exploration technology. Addressing the problems of easy divergence and unstable wavefield calculations in traditional Born series in strongly disturbed VTI media, this invention first decomposes the velocity and Thomsen anisotropy parameters of the VTI medium into background and disturbance components. A virtual component is introduced into the background vertical wavenumber to construct a complex wavenumber equation. Preprocessing of the equation is performed using a pre-solution factor to obtain a renormalized Born series that converges under strong scattering. Thin-plate partitioning is used to decompose the global calculation into local layered extensions, sequentially performing downward forward scattering extensions and upward backward scattering backpropagation, ultimately obtaining the full wavefield seismic record. This invention effectively improves the stability and accuracy of seismic wave simulation in strongly anisotropic and high-velocity-contrast media, while ensuring computational efficiency through thin-plate layered calculations, making it suitable for forward modeling of seismic waves in deep, complex structures.

Description

Technical Field

[0001] This invention belongs to the field of deep earth engineering geophysical exploration technology, specifically relating to a numerical simulation method for seismic waves in VTI media based on preprocessing iteration. Background Technology

[0002] As deep resource extraction and underground space development extend into the Earth's depths, deep-earth engineering faces a complex geological environment characterized by "high ground stress, high temperature, high osmotic pressure," as well as "strong disturbance and strong anisotropy." Strong anisotropy stems from the directional arrangement of minerals, the predominance of fractures, and the bedding structure in deep strata, resulting in a significant direction dependence of seismic wave velocity and attenuation characteristics. Strong disturbance is triggered by rock excavation and seismic source activation, disrupting the original stress balance and causing the dynamic evolution of fracture networks, leading to strong discontinuities and strong scattering characteristics in the medium.

[0003] Traditional numerical simulation methods for seismic waves are mostly based on isotropic and weak scattering assumptions, making it difficult to accurately characterize the propagation of wavefields in deep, complex VTI (vertically and laterally isotropic) media. Among existing forward modeling methods for VTI media, grid-based methods such as the finite difference method and pseudospectral method have high computational costs and require high numerical stability of anisotropic parameters; while scattering methods based on Born series are prone to series divergence under conditions of strong velocity contrast and strong anisotropy, failing to meet the requirements of high-precision simulation.

[0004] Preprocessing iteration is a renormalized scattering series method. It constructs the complex wavenumber Lippmann-Schwinger (LS) equation by introducing a small imaginary component into the background wavenumber, and then preprocesses the equation using a pre-solution factor, ensuring the Born series remains convergent under strong scattering conditions. In existing research, some industry scholars have applied this method to single-frequency wavefield calculations in acoustic media. Subsequently, Chinese invention patent CN120873329A discloses a numerical simulation method for scattering waves in viscous acoustic media based on preprocessing iteration, utilizing this method to perform numerical simulations of seismic waves in viscous acoustic media. However, the application of this technique in VTI anisotropic media has not yet been reported. Because VTI media exhibit multi-parameter coupling between velocity and two anisotropic parameters, and the wavenumber dispersion relation and scattering potential operator differ fundamentally from those of isotropic viscoacoustic media, the absorption attenuation of viscoacoustic media leads to seismic wave amplitude attenuation and velocity dispersion, resulting in relatively weak amplitude of deep structures in migration imaging results. Furthermore, the anisotropy of VTI media affects the travel time of seismic waves, leading to numerous instances of diffracted waves not returning to their original positions and structures deviating from their actual locations in the migration results. Therefore, existing preprocessing frameworks cannot be directly applied between different media, necessitating the development of preprocessing iterative numerical simulation methods adapted to VTI media to provide theoretical support for the detection of deep concealed structures and early warning of disaster sources. Summary of the Invention

[0005] The purpose of this invention is to address the shortcomings of existing technologies by providing a numerical simulation method for seismic waves in VTI media based on preprocessing iteration. This method solves the problems of easy divergence of traditional Born series and unstable wavefield calculation in strongly disturbed and strongly anisotropic media. At the same time, it improves computational efficiency through thin plate partitioning, thereby achieving high-precision and high-stability numerical simulation of seismic waves in complex VTI media.

[0006] To achieve the above objectives, the technical solution adopted by this invention is: a numerical simulation method for VTI medium seismic waves based on preprocessing iteration, comprising the following steps: S1. Model Parameter Settings: Construct the VTI medium velocity model and set the model coordinate parameters. The number of horizontal grids is The number of vertical grids is The horizontal grid spacing is The vertical grid spacing is Determine the longitudinal wave velocity Spatial distribution of Thomsen anisotropy parameters (i.e., anisotropy parameters) , ).

[0007] S2. Observation System Layout: Establish a coordinate system with the horizontal axis as the X-axis and the vertical axis as the Z-axis. Arrange geophones horizontally along the ground surface and set the channel spacing. Determine the source location, source wavelet parameters, sampling interval, and total sampling duration. The corresponding number of sampling points at each time is as follows: .

[0008] S3. Thin Plate Partitioning: The model is layered using thin plate partitioning, dividing the model along the depth direction into sections. Each thin plate has the same thickness as the longitudinal grid spacing. Each thin plate corresponds to a set of medium parameters, including velocity. Anisotropy parameters and .

[0009] S4. Downward and Forward Extension: The source wavelet is transformed to the frequency domain using a Fourier transform to obtain the frequency domain source wavelet. ;by As initial boundary conditions, the wave field is extended downwards along the depth direction, layer by layer, until the bottom of the model is reached, to obtain the forward scattering wave field and the first backscattering wave field of each depth layer.

[0010] S5. Upward and Backward Backpropagation: Using the wave field at the bottom of the model as the initial value, the backscattered wave field is extended upwards through thin plates. During the extension process, the backscattered wave field of the corresponding depth layer is superimposed until the surface, and the frequency domain result of the full wave field at the surface is obtained.

[0011] S6. Time-frequency conversion: The full-wave field in the frequency domain of the Earth's surface is converted to the time domain by inverse Fourier transform to obtain the final seismic record.

[0012] Furthermore, the specific steps of wavefield extension down the thinner plate by thinner plate in step S4 are as follows: S41, Parameter Decomposition: For the first... The dielectric parameters within a thin plate are decomposed into background and disturbance components.

[0013] S42. Preprocessing parameter calculation: Solving the... The imaginary part of the complex wave number corresponding to a thin plate With pre-solution factor .

[0014] S43. Background wave field extension: The phase shift of the input wave field is calculated in the frequency-wavenumber domain to obtain the extended background wave field.

[0015] S44. Disturbance Correction Calculation: Based on the dispersion characteristics of VTI media, the appropriate background complex vertical wavenumber expression is derived, and three independent wavenumber domain correction operators for velocity disturbance, ε parameter disturbance, and δ parameter disturbance are calculated to accurately characterize the influence of anisotropy on wavefield propagation.

[0016] S45. Scattering field calculation: The background wave field and correction terms are transformed into the frequency-space domain, and multiple forward scattering wave fields and one backscattering wave field are calculated in combination with the pre-solution factor; the equation preprocessing is completed by the pre-solution factor to realize the renormalization of the Born series.

[0017] S46. Use the forward-scattered wave field as the input boundary condition for the next thin plate and save the first backscattered wave field of the current layer. Repeat steps S41 to S45 until the downward extension of all thin plates is completed. This extension mode adopts the method of "forward multiple scattering + backscattering once". When extending downward, multiple forward scattering and one reflection at the interface are calculated. When propagating upward, the reflected waves of each layer are superimposed to balance simulation accuracy and computational efficiency.

[0018] Furthermore, the specific method of parameter decomposition in step S41 is as follows: velocity decomposition is... ,in For the first The transverse average background velocity of a thin plate Let be the velocity disturbance component, where For the first The depth coordinates of a thin plate.

[0019] Anisotropic parameters are decomposed as follows:

[0020] in , For the first Background anisotropy parameters of a thin plate , This corresponds to the disturbance component.

[0021] Furthermore, in step S42, the imaginary component... With pre-solution factor The calculation formula is:

[0022] in, It is the angular frequency; For the velocity perturbation scattering potential, the following conditions are satisfied:

[0023] in, , Depth position ( The true complex velocity parameters and background complex velocity parameters, and the imaginary components. Take 1.05 to 1.2 times the maximum value of the perturbation term modulus to ensure the convergence of the scattering series.

[0024] Furthermore, the specific method of background wave field extension in step S43 is as follows: [The text abruptly ends here, likely due to an incomplete sentence or a formatting error.] The input wave field of the layer is transformed to the frequency-wavenumber domain by Fourier transform, and the phase shift operator is used to obtain the first wave. Background wave field of the layer :

[0025] in For transverse wavenumber, The expression for the complex vertical wavenumber of the VTI medium background after introducing the imaginary component is as follows:

[0026] Complex vertical wavenumber is denoted as , actual part With the imaginary part satisfy:

[0027] In the formula , These are the real and imaginary parts of the expression within the square root, respectively.

[0028] Furthermore, the calculation method for the three types of correction terms in step S44 is as follows:

[0029] in For velocity disturbance correction term, for Parameter disturbance correction term for Parameter perturbation correction term; correction operator , , They are respectively:

[0030]

[0031]

[0032] In the formula For background P-wave number, It is the square of the normalized transverse wavenumber.

[0033] Furthermore, the scattered field in step S45 is calculated as follows: the background wave field and the three types of correction terms are transformed into the frequency-spatial domain via inverse Fourier transform, resulting in... , , , ; Calculate the forward scattered wave field multiple times using the pre-solution factor:

[0034] The first backscattered wave field is:

[0035] in,

[0036]

[0037]

[0038]

[0039]

[0040]

[0041] In the above formula, For the background wave field, For the spatial domain velocity perturbation term, and Corresponding to anisotropic parameters , The disturbance term, Let be the pre-solution factor for the j-th thin plate; and These are the Fourier transform and the inverse transform, respectively. , and These are velocity perturbation and anisotropy, respectively. Parameter perturbation and anisotropy Parameter perturbation; i is the imaginary part. Thin plate spacers.

[0042] Furthermore, the specific method of upward backscattering continuation in step S5 is as follows: the imaginary part of the background vertical wavenumber is taken as negative, and the backscattered wave field at the bottom of the model is used as the initial value. Phase shift continuation, perturbation correction, and scattering field calculation are performed plate by plate upward; every time the continuation reaches the first... When using a thin plate, a pre-stored backscattered wave field is superimposed on the layer. This serves as the input boundary condition for the next thinner plate, and so on until the surface.

[0043] Furthermore, in step S1, for sedimentary rock VTI media, the anisotropy parameters are calculated using empirical velocity formulas:

[0044] in , These represent the maximum and minimum values ​​of the longitudinal wave velocity within the model, respectively.

[0045] Furthermore, the source wavelet adopts the Rick wavelet, and the source is located at the midpoint of the horizontal direction on the ground surface; the detectors are arranged along the entire aperture of the ground surface, and the channel spacing is consistent with the horizontal grid spacing.

[0046] Compared with the prior art, the present invention has the following advantages: 1. This invention introduces preprocessing iterative renormalization technology into VTI media for the first time. By renormalizing the Born series through imaginary components and pre-solution factors, it effectively overcomes the problem of series divergence under strong velocity contrast and strong anisotropy conditions, and significantly improves the stability and simulation accuracy of wave field calculation.

[0047] 2. For the Thomsen parameter system of VTI media, this invention re-derives the complex vertical wavenumber expression and three types of perturbation correction operators, realizing the accurate characterization of anisotropic wave characteristics. It is different from the existing preprocessing framework for isotropic and viscous acoustic media and is adapted to the multi-parameter coupling characteristics of VTI media.

[0048] 3. This invention adopts the thin plate partitioning concept to transform global calculations into local layered operators, avoiding large-scale matrix operations; at the same time, it achieves high-precision scattering simulation with a small number of iterations through preprocessing iterations, effectively reducing the amount of computation while ensuring wave field accuracy.

[0049] 4. This invention, through a two-way process of downward forward extension and upward backward propagation, retains multiple forward scattering and interface reflection information, and the simulation results are more consistent with the actual wave field propagation law. Attached Figure Description

[0050] Figure 1 This is a schematic diagram of the overall process of the method of the present invention.

[0051] Figure 2 This is a schematic diagram of model parameters according to an embodiment of the present invention, wherein, Figure 2 In this context, 'a' represents the P-wave velocity model. Figure 2 In this context, b represents the distribution of the ε anisotropy parameter. Figure 2 In this context, c represents the distribution of the δ anisotropic parameter.

[0052] Figure 3 This is a schematic diagram of a numerical simulation earthquake record according to an embodiment of the present invention. Detailed Implementation

[0053] The present invention will be further described below.

[0054] Example 1: This example provides a numerical simulation method for seismic waves in VTI media based on preprocessing iteration. The specific implementation steps are as follows: S1. Model Parameter Settings: Construct a two-dimensional VTI medium model with a horizontal range of 0~13.6km and a depth range of 0~3.2km; horizontal grid number... Vertical grid number Horizontal and vertical grid spacing The P-wave velocity range within the model is 2000–4000 m / s, encompassing layered interfaces and local anomalies. Anisotropy parameters are calculated using empirical formulas to account for the characteristics of sedimentary rock VTI media.

[0055] in , The spatial distributions of the ε and δ parameters are obtained as follows: Figure 2 As shown.

[0056] S2. Observation System Layout: The coordinate origin (0,0) is set at the upper left corner of the Earth's surface. The X-axis extends horizontally to the right, and the Z-axis extends downwards along the depth. 1360 geophones are deployed along the entire aperture of the surface, with channel spacing... The first detector is located at (0,0). The seismic source uses a Ricker wavelet with a dominant frequency of... The epicenter was located at the midpoint of the horizontal axis of the Earth's surface (6800m, 0m). Time sampling intervals were set. The total sampling time is 1 second, and the corresponding number of sampling points is... .

[0057] S3. Thin Plate Division: Divide the model into 320 thin plates along the depth direction, with each thin plate having a thickness of... The spacing is consistent with the vertical grid spacing. The depth coordinates of the thin plate are The velocity at that depth is preserved within each thin plate. Anisotropy parameters and .

[0058] S4, Downward and Forward Wavefield Extension: For the Reichsson Wavelet Perform a Fourier transform to obtain the frequency domain wavelet. It covers the effective frequency band of 0~80Hz. Using the surface seismic source wavelet as the initial boundary condition, the process is as follows: Starting from the first thin plate, the overlay is extended downwards plate by plate. S41, Parameter Decomposition: For the first... The background parameters are obtained by lateral averaging of the dielectric parameters within the thin plate.

[0059]

[0060]

[0061] This leads to the perturbation components of each parameter:

[0062]

[0063]

[0064] S42. Preprocessing parameter calculation: First, calculate the velocity perturbation scattering potential:

[0065] Calculate the imaginary part of the complex wave number The value is taken as 1.1 times the maximum magnitude of the total perturbation terms to ensure series convergence:

[0066] Calculate the pre-solution factor:

[0067] S43, Background Wavefield Phase Shift Extension: [The sentence is incomplete and requires further context.] The input wave field of the layer is transformed to the frequency-wavenumber domain by Fourier transform, and then extended to the 1st layer by the VTI medium phase shift operator. Layer, to obtain the background wave field:

[0068] The background complex vertical wavenumber Satisfying VTI medium dispersion relation:

[0069] After introducing the imaginary component For a complex number, denoted as The real and imaginary parts are respectively:

[0070] In the formula , These are the real and imaginary parts of the expression within the square root of the dispersion relation, respectively.

[0071] S44. Disturbance Compensation Calculation: Calculate wavenumber domain compensation terms for velocity disturbances, ε disturbances, and δ disturbances respectively.

[0072]

[0073]

[0074] The three types of correction operators are as follows:

[0075]

[0076]

[0077] In the formula For background P-wave number, The normalized transverse wavenumber squared.

[0078] S45. Scattered Field Calculation: The background wave field and three types of correction terms are transformed to the frequency-space domain using inverse Fourier transform, and the forward scattered wave field is calculated multiple times using the pre-solution factor.

[0079] Calculation of the first backscattered wave field based on the interface reflection coefficient:

[0080] S46. Layer-by-layer recursion: The forward-scattered wave field... As the first Input boundary conditions for thin plates, storing the first backscattered wave field of the current layer. Repeat the above steps of parameter decomposition to scattering field calculation, and sequentially complete the downward extension of all 320 thin plates until the bottom of the model.

[0081] S5. Upward Backscattering Propagation: Using the wavefield at the bottom of the model as the initial value, backscattering extension is performed thinly by thin plate upwards along the depth. During upward extension, the imaginary part of the background vertical wavenumber is negative, and the phase shift operator propagates upwards accordingly. Each extension to the [missing value]... When using a thin plate, a pre-stored backscattered wave field is superimposed on the layer. This serves as the input boundary condition for the next thinner plate. This process continues until the plate is extended to the Earth's surface, yielding the full frequency domain field at the surface.

[0082] S6. Time-Frequency Conversion and Result Output: Perform an inverse Fourier transform on the frequency domain wavefield of each geophone location on the ground surface to convert it into a time domain wavefield, obtaining the seismic records of all gathers, such as... Figure 3 As shown.

Claims

1. A numerical simulation method for seismic waves in VTI media based on preprocessing iteration, characterized in that, Includes the following steps: S1. Construct a VTI medium velocity model, set model coordinate parameters and mesh parameters; determine the longitudinal wave velocity and Thomsen anisotropy parameters at each mesh point in the model; S2. Establish a coordinate system and set the detector arrangement, source location, wavelet parameters, and time sampling parameters; S3. The model is divided into several thin plates with a thickness equal to the longitudinal grid spacing along the depth direction using a thin plate partitioning method, and the global calculation is decomposed into local layered calculations. S4. Transform the source wavelet to the frequency domain, extend the wave field downwards plate by plate along the depth direction, and complete parameter decomposition, preprocessing parameter calculation, background wave field phase shift, disturbance correction and forward scattering field calculation for each plate. Save the backscattering wave field of each layer until the bottom of the model. S5. Using the wave field at the bottom of the model as the initial condition, the backscattered wave field is extended layer by layer along the depth upwards. During the extension process, the backscattered wave field of the corresponding depth layer is superimposed until the surface, and the frequency domain result of the full wave field of the surface is obtained. S6. Inverse transform the surface frequency domain wavefield to the time domain and output the final seismic record.

2. The method according to claim 1, characterized in that, The specific steps of wavefield extension down the thinner plate by thinner plate in step S4 are as follows: S41, Parameter Decomposition: For the first... The dielectric parameters within a thin plate are decomposed into background and disturbance components. S42. Preprocessing parameter calculation: Solving the... The imaginary part of the complex wave number corresponding to a thin plate With pre-solution factor ; S43, Background wave field extension: The phase shift of the input wave field is calculated in the frequency-wavenumber domain to obtain the extended background wave field; S44. Disturbance Correction Calculation: Calculate velocity disturbance, Parameter perturbation, Three types of wavenumber domain correction terms corresponding to parameter perturbations; S45. Scattered field calculation: Transform the background wave field and correction term to the frequency-space domain, and calculate multiple forward scattered wave fields and one backscattered wave field in combination with the pre-solution factor. S46. Use the forward-scattered wave field as the input boundary condition for the next thin plate, and save the backscattered wave field of the current layer; repeat steps S41 to S45 until the downward extension of all thin plates is completed.

3. The method according to claim 2, characterized in that, The specific method of parameter decomposition in step S41 is as follows: velocity decomposed into ,in For the first The transverse average background velocity of a thin plate Let be the velocity disturbance component, where For the first The depth coordinates of the thin plate; Anisotropic parameters are decomposed as follows: in , For the first Background anisotropy parameters of a thin plate , This corresponds to the disturbance component.

4. The method according to claim 3, characterized in that, In step S42, the imaginary component With pre-solution factor The calculation formula is: in, It is the angular frequency; For the velocity perturbation scattering potential, the following conditions are satisfied: in, , Depth position ( The true complex velocity parameters and background complex velocity parameters, and the imaginary components. Take 1.05 to 1.2 times the maximum value of the perturbation term modulus to ensure the convergence of the scattering series.

5. The method according to claim 4, characterized in that, The specific method for background wave field extension in step S43 is as follows: [The text abruptly ends here, likely due to an incomplete sentence or a missing section.] The input wave field of the layer is transformed to the frequency-wavenumber domain by Fourier transform, and the phase shift operator is used to obtain the first wave. Background wave field of the layer : in For the thickness of a single sheet, For transverse wavenumber, The expression for the complex vertical wavenumber of the VTI medium background after introducing the imaginary component is as follows: Complex vertical wavenumber is denoted as , actual part With the imaginary part satisfy: In the formula , These are the real and imaginary parts of the expression within the square root, respectively.

6. The method according to claim 5, characterized in that, The calculation method for the three types of correction terms in step S44 is as follows: in For velocity disturbance correction term, for Parameter disturbance correction term for Parameter perturbation correction term; correction operator , , They are respectively: In the formula For background P-wave number, It is the square of the normalized transverse wavenumber.

7. The method according to claim 6, characterized in that, The scattered field in step S45 is calculated as follows: the background wave field and the three types of correction terms are transformed into the frequency-space domain via inverse Fourier transform, resulting in... , , , ; Calculate the forward scattering wave field multiple times using the pre-solution factor : First backscattered wave field for: in, For the background wave field, For the spatial domain velocity perturbation term, and Corresponding to anisotropic parameters , The disturbance term, Let be the pre-solution factor for the j-th thin plate.

8. The method according to claim 1, characterized in that, The specific method for upward backscattering continuation in step S5 is as follows: the imaginary part of the background vertical wavenumber is taken as negative, and the backscattered wave field at the bottom of the model is used as the initial value. Phase shift continuation, perturbation correction, and scattering field calculation are performed plate by plate upward; every time the continuation reaches the first... When using a thin plate, a pre-stored backscattered wave field is superimposed on the layer. This serves as the input boundary condition for the next thinner plate, and so on until the surface.

9. The method according to claim 1, characterized in that, In step S1, for sedimentary rock VTI media, the anisotropy parameters are calculated using empirical velocity formulas: in, For the P-wave velocity at the grid points, , The Thomsen anisotropy parameters for the grid points are... , These represent the maximum and minimum values ​​of the longitudinal wave velocity within the model, respectively.

10. The method according to claim 1, characterized in that, The source wavelet is the Rick wavelet, and the source is located at the midpoint of the horizontal direction on the ground surface; the geophones are arranged along the entire aperture of the ground surface, and the channel spacing is consistent with the horizontal grid spacing.

Images