Advanced tunnel detection method and system for viscous-acoustic medium reverse-time migration of high-order staggered grids

Abstract

The invention relates to a high-order staggered grid viscous sound medium reverse time migration tunnel advanced detection method and system, and belongs to the field of tunnel engineering geophysical prospecting. According to the method, a viscoelastic anisotropic TTI medium anisotropy wave equation expressed by a fractional Laplace operator is deduced, and a viscoelastic anisotropic TTI medium first-order speed-stress elastic wave equation is solved by adopting a high-order staggered grid finite difference method. According to the invention, the wave field information at any moment is obtained, so that the propagation process of seismic waves in an underground medium is simulated more accurately. In order to solve the problem of reflection possibly caused by artificial boundaries, a PML (perfectly matched layer) absorption boundary processing method of a new attenuation function is adopted. According to the method, seismic waves can be effectively absorbed, reflection is avoided, and therefore simulation stability and accuracy are guaranteed. Generally speaking, the invention provides an accurate and stable tunnel advanced detection forward modeling method, and a reverse time migration algorithm is carried out in an attenuation TTI medium.

Description

Technical Field

[0001] The present invention relates to a reverse time migration tunnel advanced detection method and system for visco-acoustic media with high-order staggered grids, belonging to the field of geophysical prospecting for tunnel engineering. Background Art

[0002] As an important part of the construction of modern transportation, water conservancy, mines and other infrastructure, the construction safety and efficiency of tunnels are directly related to the overall quality and progress of the project. However, during the tunnel construction process, complex geological conditions are often encountered, such as faults, fracture zones, karst caves, groundwater and other unfavorable geological bodies. These geological anomalies may cause disasters such as collapses and water inrushes, seriously threatening construction safety. Therefore, accurately detecting the geological conditions ahead and predicting potential risks before tunnel excavation has become an indispensable key link in tunnel engineering. With the development of geophysical exploration technologies, seismic wave detection methods have gradually become one of the mainstream technologies for tunnel advanced detection due to their advantages such as large detection depth, high resolution, and wide application range. The seismic wave detection method artificially generates seismic waves and uses the law of seismic wave propagation in underground media to image geological structures. However, in practical applications, seismic waves will be affected by the absorption and attenuation of the medium during propagation, resulting in a decrease in wave field energy and phase distortion. This phenomenon is particularly significant in visco-acoustic media (such as aquifers and soft rock formations). This attenuation effect will reduce the resolution and accuracy of seismic wave imaging, and thus affect the ability to identify geological anomaly bodies.

[0003] Reverse Time Migration (RTM) is a high-precision seismic imaging algorithm. It can achieve high-resolution imaging of complex geological structures through the matching of two-way wave field propagation and imaging conditions. Combining the high-order staggered grid finite difference method with reverse time migration can significantly improve the accuracy and reliability of seismic wave imaging in visco-acoustic media. It provides a scientific basis and safety guarantee for tunnel construction and has important engineering application value. Summary of the Invention

[0004] The technical problem solved by the present invention is: The present invention provides a reverse time migration tunnel advanced detection method and system for visco-acoustic media with high-order staggered grids to solve the influence of amplitude loss and phase dispersion during seismic wave propagation on seismic imaging caused by the attenuation and anisotropy characteristics of real earth media, and the reflection problems that may be caused by artificial boundaries. The present invention provides an accurate and stable forward simulation method for tunnel advanced detection, realizes the reverse time migration algorithm in attenuating TTI media, and the method of the present invention can effectively absorb seismic waves and avoid the generation of reflections, thus ensuring the stability and accuracy of the simulation.

[0005] The technical solution of the present invention is: A reverse time migration tunnel advanced detection method for visco-acoustic media with high-order staggered grids, the method comprising:

[0006] Step 1: Construct a velocity model, set the spatial and temporal sampling intervals, and then discretize the velocity model to prepare for subsequent numerical simulations.

[0007] Step 2: Based on the SLS model (Standard Linear Solid), a viscoacoustic wave equation is established to characterize the propagation characteristics of seismic waves in attenuating media; the SLS theory is combined with the pure qP wave dispersion relation to derive the pure viscous acoustic wave TTI wave equation represented by memory variables;

[0008] Step 3: Use the derived wave equation to simulate the propagation of seismic waves inside the velocity model to obtain the wave field data of each grid point at each time; calculate the wave field values ​​of the grid points on the boundary of the velocity model at each time through the wave field absorption boundary equation and the wave field absorption coefficient, so as to obtain the boundary wave field data;

[0009] Step 4: Combine the boundary wave field data, the internal wave field data and the corner wave field data to obtain complete wave field data, output the final wave field snapshot of the velocity model, and perform reverse time migration using source normalized correlation imaging conditions.

[0010] Furthermore, the Step 2 includes:

[0011] Based on the constitutive relationship of the SLS model, the viscoacoustic wave equation is derived by considering the motion equation and constitutive relationship of the medium, and combined with the pure qP wave dispersion relationship, and then the memory variable is introduced to simplify the solution of the acoustic wave equation, and finally the pure viscous acoustic wave TTI (transversely isotropic) wave equation represented by the memory variable is obtained; there are many kinds of relationships between stress and strain in viscoelasticity, and there are different models according to different assumed relationships. The SLS model usually refers to the Standard Linear Solid Model, also known as the Zener model. It is a classic viscoelastic material model consisting of the following two parts: the spring represents the elastic properties of the material and follows Hooke's law. The spring-damper series represents the viscous properties of the material, which is used to describe the characteristics of the material that exhibits both elastic and viscous behavior when subjected to force. Specifically include:

[0012] Step 2.1. Based on the SLS model, establish the viscoacoustic wave equation to characterize the propagation characteristics of seismic waves in the attenuation medium;

[0013] Step 2.2, based on the established viscoacoustic wave equation, by combining it with the pure qP wave dispersion relation, the pure viscous acoustic wave TTI wave equation represented by the memory variable is derived;

[0014] Step 2.3. Decouple the two-dimensional pure-viscous acoustic TTI wave equation containing the coupling amplitude attenuation and phase dispersion operators derived in Step 2.2, and derive a pure-viscous acoustic TTI wave equation containing the decoupled amplitude attenuation and phase dispersion operators.

[0015] Further, the Step 2.1 includes:

[0016] Use a single SLS element combined with the viscoacoustic wave equation to characterize the seismic wave properties in an attenuating medium; for the SLS model, the constitutive equation of the SLS model is:

[0017]

[0018] where p1, q0, q1 are coefficients obtained through the forward and inverse Laplace transforms, σ represents stress, which is the force per unit area inside the material, represents the rate of change of stress with respect to time, that is, the derivative of stress ε represents strain, which is the deformation of the material under the action of force, represents the rate of change of strain with respect to time, that is, the derivative of strain

[0019] Derive the viscoacoustic wave frequency-wavenumber domain equation through the constitutive relationship of the SLS model:

[0020]

[0021] where p is the pressure wave field, r is the memory variable, P and are the spatial Fourier transforms of p and r respectively, k x , k y , k z are the wave fields in x, y, and z, and are the stress relaxation time and strain relaxation time respectively, where Q is the quality factor, ω represents the angular frequency, ω0 is the reference angular frequency, and in the Fourier transform, the time derivative corresponds to iω in the frequency domain, and iω represents the derivative with respect to time.

[0022] Further, the Step 2.2 includes:

[0023] Characterize the seismic wave behavior in an anisotropic medium through the dispersion relations of P-waves and SV-waves, set the phase velocity of SV-waves to zero at all phase angles, add the viscous effect, and transform to the time-space domain to derive the pure-viscous acoustic TTI wave equation:

[0024]

[0025] In the formula, The variable v represents the change of the pressure field p with time. p0 is the P-wave velocity, respectively represent the second-order partial derivatives of the pressure with respect to the spatial coordinates x, y, and z. S is used to describe the wave propagation characteristics in the medium and is related to anisotropy and viscosity. The variable represents the evolution of the memory variable with time.

[0026] Furthermore, the said Step 2.3 includes:

[0027] Combining the coordinate transformation of VTI and TTI media with the pure viscous acoustic wave TTI wavefield propagation operator equation:

[0028] W vis = W ac + W phase + W amp

[0029] where W ac is the pure acoustic wave TTI wavefield propagation operator, W phase is the phase dispersion control operator, and W amp is the amplitude loss control operator;

[0030] By decoupling through the addition and subtraction operations of these three operators, the pure viscous acoustic wave TTI wave equation with decoupled amplitude attenuation and phase dispersion operators based on the SLS model is derived:

[0031]

[0032] where k1 and k2 are the coefficients of the control phase dispersion term and the control amplitude loss term respectively. When k1 = 0 and k2 = 1, the equation is the pure viscous acoustic wave TTI wave equation dominated by phase dispersion. When k1 = 1 and k2 = 0, the equation is the pure viscous acoustic wave TTI wave equation dominated by amplitude loss. When k1 = 1 and k2 = 1, the equation includes both phase dispersion and amplitude loss terms and can comprehensively describe the propagation characteristics of viscous acoustic waves in TTI media.

[0033] Furthermore, the said Step 3 includes:

[0034] Performing forward modeling on the new wave equation and integrating the perfectly matched layer, and using the high-order staggered grid finite difference method as the numerical solution method to solve the pure viscous acoustic wave TTI wave equation to obtain the difference iteration equation at any point in the simulation area; the specific steps include:

[0035] Step 3.1: According to the one-dimensional second-order viscoacoustic medium wave equation derived in Step 2, determine the staggered grid finite difference format. Expand the computational variable matrix by a certain number of layers around as the perfectly matched layer (PML) for absorbing seismic waves. The main region is the computational interval and also the target region for simulation. The simulation object is the tunnel isomer. The expanded range is the PML region of the perfectly matched layer.

[0036] Step 3.2: Use the high-order staggered grid finite difference formula derived in Step 3.1 as the numerical solution method to solve the pure viscoacoustic wave TTI wave equation to obtain the difference iteration equation at any point in the simulation region. Discretize the entire new wave equation derivation formula to meet the requirements of numerical simulation, and then perform iteration to obtain a seismic wave simulation in which the seismic wave energy can rapidly decay in the PML region.

[0037] Furthermore, Step 3.1 includes:

[0038] The PML region is divided into three sub-regions: the upper and lower boundary regions with the x-axis as the main direction, the left and right boundary regions with the z-axis as the main direction, and the corner intervals at the four corners; each sub-region will independently apply the attenuation factor α according to its own characteristics and requirements to achieve precise and efficient seismic wave absorption.

[0039]

[0040] Among them, α(Z) represents the attenuation factor in the PML region, which is used to describe the attenuation intensity of seismic waves in the PML region. Z represents the distance from the boundary of the PML region. β is a parameter that controls the shape of the attenuation factor and is used to adjust the smoothness of the attenuation curve. Make the function meet the basic requirements of the attenuation function. Fine-tune the function, K, γ, δ are adjustable coefficients, and K is a known number, n is the order, L is the thickness of the PML layer; R is the theoretical boundary reflection coefficient; ρ is the medium density, and u is the elastic modulus of the medium.

[0041] Assign values to the frequency shift factor η x and the scale factor β x in the PML region. In the main region, the frequency shift factor η x takes the value of 0, and the scale factor β x is 1. The calculation formula for the scale factor β x is as follows;

[0042]

[0043] Among them, β0 represents the reference value of the scale factor β x and η0 represents the frequency shift factor η xThe reference value, x is the distance measured from the boundary of the PML region, P η and P β is the power exponent that controls η x and β x 's rate of change, and are normalized power functions that respectively describe η x and β x 's variation law with position x, f is a function used to adjust the η x distribution;

[0044] For the second-order acoustic wave equation, the wave field function is set to Let Δx, Δz, Δt represent the spatial step sizes in the x and z directions and the time step size respectively. Let i, j, k be the grid numbers in space and time, then there is Using the Taylor formula to expand p at time t gives and Derive the coefficient matrix equation:

[0045]

[0046] Solve for the difference coefficient C n , select the absorption function and substitute the Ricker wavelet as the boundary condition; the expression formula of the absorption function is as follows:

[0047]

[0048] Among them, L is the total number of layers in the absorption region, k is the number of layers from the forward modeling region, d(k) is the absorption function used to describe the absorption intensity of the k-th layer in the PML region, N is the total number of layers in the absorption region, V p is the P-wave velocity, R = 0.0001;

[0049] The expression formula of the Ricker wavelet is:

[0050]

[0051] Among them, r(t) represents the amplitude value of the Ricker wavelet at time t, f m is the main frequency of the Ricker wavelet;

[0052] Solve for the difference coefficient of the first-order derivative with 2N-order accuracy. The high-precision operator of the first-order derivative is expressed as:

[0053]

[0054] In the formula, and respectively represent the offsets from the grid point i in the x direction by and The pressure value of a grid unit;

[0055] The difference formula is:

[0056]

[0057] where f (2) (x) represents the second-order derivative of f(x) in the x direction, a0 is the central point coefficient, f(x0) represents the value of f(x) at x0, and f is the velocity and stress component to be differentiated; Δx and Δz are the spatial step sizes in the x and z directions; t is the simulation time; Δt is the time step; a m is the coefficient for the high-order accurate difference approximation of the second-order derivative, m and M are positive integers from one to infinity, and x0 is the unknown in the process of partial derivative calculation, represents the change in the function value within the time step Δt.

[0058] Furthermore, the Step 3.2 includes:

[0059] The difference iteration equation for any point includes obtaining the iteration formula for the pressure field by combining the spatial derivative of the velocity field and the influence of the memory variable r:

[0060]

[0061] where is the value of the pressure field at time step n + 1, and are the spatial derivatives of the velocity field, and K is the bulk modulus;

[0062] Combining the partial derivative of the pressure field p with respect to x obtains the iteration formula for the velocity component V x in the x direction:

[0063]

[0064] where and represent the velocity components V x at the same grid point but different time steps;

[0065] Combining the partial derivative of the pressure field p with respect to z obtains the iteration formula for the velocity component V z in the z direction:

[0066]

[0067] where and represent the velocity components V z .

[0068] Furthermore, Step 4 includes:

[0069] Using the pure viscous acoustic wave TTI wave equation obtained in Step 2, reverse the amplitude loss term to achieve attenuation compensation, and use the regularization method to suppress the high-frequency components generated by the inversion, to obtain an attenuation compensation wavefield modeling operator with a regularization term:

[0070]

[0071] After determining the appropriate regularization parameter ζ, use the source-normalized cross-correlation imaging condition as the imaging condition. The expression of the source-normalized cross-correlation imaging condition is:

[0072]

[0073] where I c represents the result of the source-normalized cross-correlation imaging condition, α and L respectively represent the attenuation coefficient and the propagation length; the subscripts up and down respectively represent the up-going wave and the down-going wave; S(x,t) and R(x,t) respectively represent the source wavefield and the received wavefield in the non-attenuating anisotropic medium, (x,t) represents the spatial coordinate x and the time t, S A (x,t) represents the source wavefield in the attenuating medium, R A (x,t) represents the received wavefield in the attenuating medium, R C (x,t) represents the compensated received wavefield, +αt represents the variation of the attenuation compensation intensity with time, e +αt is an exponential attenuation compensation factor, used to compensate for the energy lost by the wavefield due to attenuation during propagation, -αL down describes the variation of the attenuation intensity with the propagation length, is used to describe the energy attenuation of the down-going wave during propagation, represents the energy lost by the up-going wave due to attenuation during propagation, represents the energy recovered for the up-going wave lost due to attenuation.

[0074] The present invention also provides a visco-acoustic medium reverse time migration tunnel advanced detection system with a high-order staggered grid, and the system includes: a module for executing the visco-acoustic medium reverse time migration tunnel advanced detection method with a high-order staggered grid described above.

[0075] The present invention also provides an electronic device, including a memory, a processor, and a computer program stored on the memory and executable on the processor, and when the processor executes the program, it implements the visco-acoustic medium reverse time migration tunnel advanced detection method with a high-order staggered grid described above.

[0076] The beneficial effects of the present invention are:

[0077] 1. The reverse time migration method is used for tunnel advanced detection imaging. Compared with the traditional Kirchhoff imaging technology, this technology has a more superior imaging effect on large dip interfaces, significantly reducing the occurrence of the "arc-drawing" phenomenon; the reverse time migration method can better handle the propagation and reflection of seismic waves under complex geological conditions, providing more accurate and reliable image information for tunnel detection;

[0078] 2. By combining means such as staggered grids, variable step-size grids, optimized PML absorption boundaries, and spatio-temporal domain differential coefficients, the time complexity and space complexity of finite difference calculations are effectively reduced; this improvement makes the algorithm more efficient and more suitable for lightweight embedded formula systems with limited resources; therefore, even in harsh working environments, the present invention can operate stably and efficiently, providing technical support for tunnel detection;

[0079] 3. The normalized cross-correlation imaging condition is adopted. This innovation helps to reduce the ray noise and dispersion noise that may occur during the imaging process. By optimizing the imaging condition, the present invention significantly improves the signal-to-noise ratio of the tunnel detection imaging result, making the imaging effect clearer and more accurate. This improvement not only enhances the imaging quality of tunnel detection but also provides a more reliable data basis for subsequent geological analysis and tunnel design. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] Figure 1 is the flow chart of the present invention;

[0081] Figure 2 is the staggered grid finite difference schematic diagram of the first-order velocity-stress elastic wave equation of the present invention;

[0082] Figure 3 is the division diagram of the perfectly matched layer PML absorption boundary region of the present invention;

[0083] Figure 4 is the schematic diagram of the reverse time migration imaging principle of the present invention;

[0084] Figure 5 is the simulation model diagram of the present invention;

[0085] Figure 6 is the wave field snapshot diagram sampled every 20 ms of the present invention;

[0086] Figure 7 is the general reverse time migration effect diagram;

[0087] Figure 8 is the reverse time migration effect diagram of the model of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0089] Example 1: As Figures 1 - 8As shown, a high-order staggered grid viscoacoustic medium reverse time migration tunnel advance detection method, the specific steps of the method are as follows:

[0090] Step 1: Construct a velocity model and set the spatial and temporal sampling intervals, then discretize the velocity model to prepare for subsequent numerical simulations;

[0091] like Figure 5 As shown, a groove velocity model is constructed, the upper and lower layers are set at different speeds, the sampling interval is 10 meters, and the time step is 1ms;

[0092] Step 2: Based on the SLS model (Standard Linear Solid), a viscoacoustic wave equation is established to characterize the propagation characteristics of seismic waves in attenuating media; the SLS theory is combined with the pure qP wave dispersion relation to derive the pure viscous acoustic wave TTI wave equation represented by memory variables;

[0093] In underground media, viscoelastic properties are often distributed in a linear relationship. Linear viscoelasticity represents the phenomenon that stress and strain are always interdependent during the entire deformation process. This stress-strain relationship can be established by a function that changes with time. The attenuation characteristics of underground media can be described by standard linear solids (SLS). The viscoelastic equation based on the attenuation mathematical model can be used for the propagation of seismic waves in attenuating anisotropic media, but its amplitude loss and phase dispersion are coupled, which is not conducive to the realization of RTM. Therefore, by combining SLS with the pure qP wave dispersion relationship, a pure viscous TTI wave equation represented by a memory variable is derived, which has a decoupling function and higher computational efficiency.

[0094] Based on the constitutive relation of the SLS model, the viscoacoustic wave equation is derived by considering the motion equation and constitutive relation of the medium, and then combined with the pure qP wave dispersion relation. Then, the memory variable is introduced to simplify the solution of the acoustic wave equation, and finally the pure viscous acoustic wave TTI wave equation represented by the memory variable is obtained.

[0095] Step 2 specifically includes:

[0096] Step 2.1. Based on the SLS model, establish the viscoacoustic wave equation to characterize the propagation characteristics of seismic waves in the attenuation medium;

[0097] SLS solid model usually refers to Standard Linear Solid Model, also known as Zener model. It is a classic viscoelastic material model, consisting of the following two parts: the spring represents the elastic properties of the material and follows Hooke's law. The spring-damper series represents the viscous properties of the material, which is used to describe the characteristics of the material that exhibits both elastic and viscous behavior when subjected to force;

[0098] The characteristics of seismic waves in an attenuating medium can be characterized by using a single SLS element in combination with the viscoacoustic wave equation; for the SLS model, its constitutive equation is:

[0099]

[0100] where p1, q0, and q1 are coefficients obtained through the forward and inverse Laplace transforms, σ represents stress, which is the force per unit area inside the material, represents the rate of change of stress with respect to time, i.e., the derivative of stress ε represents strain, which is the deformation of the material under the action of force, represents the rate of change of strain with respect to time, i.e., the derivative of strain

[0101] The equation in the viscoacoustic frequency-wavenumber domain is obtained through the SLS constitutive relation:

[0102]

[0103]

[0104] where p is the pressure wave field, r is the memory variable, P and are the spatial Fourier transforms of p and r respectively, k x , k y , k z are the wave fields in the x, y, and z directions, and are the stress relaxation time and strain relaxation time respectively, where Q is the quality factor, ω represents the angular frequency, ω0 is the reference angular frequency, and in the Fourier transform, the time derivative corresponds to iω in the frequency domain, and iω represents the derivative with respect to time.

[0105] Step 2.2: According to the viscoacoustic wave equation established in Step 2.1, by combining it with the dispersion relation of the pure qP wave, the pure viscous acoustic wave TTI wave equation expressed by the memory variable is derived.

[0106] Due to the anisotropy of the underground medium, the propagation speeds of seismic waves in different directions (P waves and S waves) are different. When the S wave passes through the anisotropic medium, shear wave splitting occurs, and the dispersion relation of seismic waves becomes more complex. The behavior of seismic waves in the anisotropic medium is characterized by the dispersion relations of P waves and SV waves; the expressions for the exact dispersion relations of P waves and SV waves are:

[0107]

[0108] where Equations (4) and (5) represent the dispersion relations of P waves and SV waves respectively, v p0and v s0 represent the P-wave and SV-wave velocities respectively, ω p and ω sv represent the angular frequencies of the P-wave and SV-wave respectively. E is an intermediate variable that describes the dispersion relationship between the P-wave and SV-wave. δ and ε are Thomsen anisotropy parameters. By setting the phase velocity v of the SV-wave sv to zero at all phase angles, a pure qP-wave dispersion formula containing only three parameters is obtained:

[0109]

[0110] In equations (2) and (3), is applicable to isotropic media. When sound waves pass through anisotropic media underground, dispersion occurs. (1 + 2ε reflects the anisotropic characteristics of the medium in the x and y directions. S k is a term related to viscosity, describing the energy loss of waves during propagation in viscous media. By replacing with (1 + 2ε the anisotropic effect is introduced. By introducing the term the viscous effect is introduced. When ε = 0 and S k = 0,

[0111] In equations (2) and (3), is replaced with Thus, the anisotropic effect and viscous effect are introduced into the dispersion relationship of the pure qP-wave, and a pure viscous acoustic wave equation is obtained in the frequency-wavenumber domain:

[0112]

[0113] By transforming equations (9) and (10) into the time-space domain, the pure viscous acoustic wave VTI wave equation is obtained:

[0114]

[0115] In the formula, represents the variable of the pressure field p changing with time, represent the second-order partial derivatives of the pressure with respect to the spatial coordinates x, y, and z respectively. S is used to describe the propagation characteristics of waves in the medium and is related to anisotropy and viscosity, represents the variable of the memory variable evolving with time. δ is one of the Thomsen anisotropy parameters and is used to describe the anisotropic characteristics of transversely isotropic (VTI) media.

[0116] Among them, equations (11) and (12) can be solved using the finite difference method. Equation (13) is the spatial domain representation of equation (8) and can only be solved by a spectral-based method. To solve equation (13) using the finite difference method subsequently, the equivalent form of equation (8) is:

[0117]

[0118] where n = (n x , n y , n z ) = k / |k| = (k x , k y , k z ) / |k| represents the propagation direction of seismic waves, and n x , n y , n z correspond to in the spatial domain respectively. The expression formula of formula (14) in the spatial domain is:

[0119]

[0120] where S n represents the propagation characteristics of waves in the spatial domain; substituting S in equations (11) and (12) with S n , a new time - space domain pure viscous acoustic VTI wave equation is obtained:

[0121]

[0122] Through the above transformation, equations (16) and (17) can be solved using the finite difference method. By rotating the grid to align the grid with the symmetry axis of the medium, the transformed equation is:

[0123]

[0124] where θ and represent the tilt angle and azimuth angle respectively. and are the new wave numbers in the rotated coordinate system. According to equation (18), the two - dimensional pure viscous acoustic TTI wave equation can be expressed as:

[0125]

[0126] In the formula, is the approximate form of the anisotropic viscous term, which describes the propagation characteristics of waves in the tilted coordinate system, T x , T z is the first - order differential operator, and T xx is the second - order differential operator, which is used to describe the propagation characteristics of waves in the tilted coordinate system.

[0127] Step 2.3. Based on the two-dimensional pure-viscous acoustic TTI wave equation containing the coupled amplitude attenuation and phase dispersion operators derived in Step 2.2, decouple it to derive a pure-viscous acoustic TTI wave equation containing the decoupled amplitude attenuation and phase dispersion operators. The steps are as follows:

[0128] Pure-viscous acoustic TTI wavefield propagation operator W vis Can generally be expressed as:

[0129] W vis = W ac + W phase + W amp (23)

[0130] Where W ac is the pure acoustic TTI wavefield propagation operator, W phase is the phase dispersion control operator, and W amp is the amplitude attenuation control operator. The phase dispersion-dominated pure-viscous acoustic TTI wavefield propagation operator W vpd is:

[0131] W vpd = W ac + W phase (24)

[0132] The amplitude attenuation-dominated pure-viscous acoustic TTI wavefield propagation operator W vamp is expressed as:

[0133] W vamp = W ac + W amp (25)

[0134] Subtracting the phase dispersion control operator W vis from the pure-viscous acoustic TTI propagation operator W phase , the amplitude attenuation-dominated pure-viscous acoustic TTI wavefield propagation operator W vamp can be indirectly obtained:

[0135] W vamp = W vis - W phase (26)

[0136] Based on the above analysis, by deriving the phase dispersion-dominated and amplitude attenuation-dominated pure-viscous acoustic TTI wavefield propagation operators, a standard linear solid (SLS) model can be used to further derive a pure-viscous acoustic TTI wave equation containing the decoupled amplitude attenuation and phase dispersion terms.

[0137] First, calculate the memory variable related to the relaxation time

[0138]

[0139] Substituting Equation (27) into Equation (9), a new pure-viscous acoustic wave equation in the frequency-wavenumber domain is obtained:

[0140]

[0141] where can be rewritten as and and correspond to the phase dispersion factor W phase and the amplitude attenuation factor W vamp respectively, where A is a parameter used to describe the phase dispersion and amplitude attenuation of the medium, Q is the quality factor, Converting Equation (28) to the time-space domain, the pure-viscous acoustic wave TTI wave equation in the time-space domain can be obtained:

[0142]

[0143] The pure-viscous acoustic wave TTI wave equation dominated by phase dispersion is expressed as:

[0144]

[0145] The pure-viscous acoustic wave TTI wave equation dominated by amplitude attenuation is expressed as:

[0146]

[0147] Based on Equations (30) and (31), a pure-viscous acoustic wave TTI wave equation with decoupled amplitude attenuation and phase dispersion based on the SLS model is finally obtained:

[0148]

[0149] where k1 and k2 are the coefficients for controlling the phase dispersion term and the amplitude loss term respectively. When k1 = 0 and k2 = 1, the equation is the pure-viscous acoustic wave TTI wave equation dominated by phase dispersion. When k1 = 1 and k2 = 0, the equation is the pure-viscous acoustic wave TTI wave equation dominated by amplitude loss. When k1 = 1 and k2 = 1, the equation contains both the phase dispersion and amplitude loss terms, and can comprehensively describe the propagation characteristics of viscous acoustic waves in TTI media.

[0150] Step3: Use the derived wave equation to simulate the propagation of seismic waves inside the velocity model. In this way, the wave field data at each grid point at each moment are obtained; through the wave field absorption boundary equation and the wave field absorption coefficient, the wave field values of the grid points on the boundary of the velocity model at each moment are calculated, so as to obtain the boundary wave field data;

[0151] In Step 3, forward modeling of the new wave equation is performed and the perfectly matched layer is integrated. The high-order staggered grid finite difference method is used as the numerical solution method to solve the pure visco-acoustic TTI wave equation, and the difference iteration equation of any point in the simulation area is obtained.

[0152] Step 3.1: According to the one-dimensional second-order visco-acoustic medium wave equation derived in Step 2, determine the staggered grid finite difference format. Expand the computational variable matrix by a certain number of layers around as the perfectly matched layer for absorbing seismic waves. The main area is the computational interval and also the target area of the simulation. The simulation object is the tunnel isomer, and the expanded range is the perfectly matched layer PML area;

[0153] The PML area is divided into three sub-areas: the upper and lower boundary areas with the x-axis as the main direction, the left and right boundary areas with the z-axis as the main direction, and the corner intervals at the four corners; Each sub-area will independently apply the attenuation factor α according to its own characteristics and requirements to achieve more accurate and efficient seismic wave absorption. The PML absorption schematic diagram is as Figure 3 shown;

[0154]

[0155] In the above formula, α(Z) represents the attenuation factor in the PML area, which is used to describe the attenuation intensity of seismic waves in the PML area. Z represents the distance from the PML area boundary, and β is a parameter that controls the shape of the attenuation factor and is used to adjust the smoothness of the attenuation curve. Make the function meet the basic requirements of the attenuation function. Fine-tune the function. In the formula, K, γ, δ are adjustable coefficients, and K is a known number, n is the order, L is the layer thickness of the PML; R is the theoretical boundary reflection coefficient, ρ is the medium density, and u is the elastic modulus of the medium.

[0156] Assign values to the frequency shift factor and the scale factor in the PML area. In the main area, the frequency shift factor takes the value of 0 and the scale factor is 1. The specific formula is as follows, β x is the scale factor, η x is the frequency shift factor;

[0157]

[0158] β0 represents the reference value of the scale factor β x and η0 represents the reference value of the frequency shift factor η x The reference value, x is the distance measured from the PML area boundary, P η and P β are power exponents that control the change rate of η x and β x respectively. and are normalized power functions, respectively describing the variations of η x and β x with position x. f is a function used to adjust the η x distribution.

[0159] For the second-order acoustic wave equation, the wave field function is set to Let Δx, Δz, and Δt represent the spatial step sizes in the x and z directions and the time step size respectively. Let i, j, and k be the spatial and time grid numbers. Then there is Expand p at time t using the Taylor formula to obtain and Derive the coefficient matrix equation:

[0160]

[0161] Solve for the difference coefficient C n , select the absorption function and substitute the Ricker wavelet as the boundary condition. The expression of the absorption function is as follows:

[0162]

[0163] In the formula, k is the number of layers from the forward modeling region, d(k) is the absorption function, used to describe the absorption intensity of the k-th layer in the PML region. N is the total number of layers in the absorption region, V p is the P-wave velocity, and R = 0.0001.

[0164] The expression of the Ricker wavelet is:

[0165]

[0166] In the formula, r(t) represents the amplitude value of the Ricker wavelet at time t, and f m is the main frequency of the Ricker wavelet.

[0167] Solve for the difference coefficient of the first-order derivative with 2N-order accuracy. The high-order accuracy operator of the first-order derivative can be expressed as:

[0168]

[0169] In the formula, and respectively represent the pressure values at the grid points offset by and grid units in the x direction from grid point i.

[0170] Use the high-order staggered grid finite difference method as the numerical solution method to solve the pure visco-acoustic TTI wave equation to obtain the difference iteration equation for any point in the simulation region. Among them, the schematic diagram of the high-order staggered grid is asFigure 2 As shown, the difference formula is:

[0171]

[0172] In the formula, f (2) (x) represents the second-order derivative of f(x) in the x direction, a0 is the central point coefficient, f(x0) represents the value of f(x) at x0, f is the velocity and stress component to be differentiated; Δx, Δz are the spatial step sizes in the x and z directions, t is the simulation time, Δt is the time step, a m is the coefficient for the high-order accurate difference approximation of the second-order derivative, m and M are positive integers from one to infinity, x0 is the unknown in the process of partial derivative, represents the change in the function value within the time step Δt.

[0173] Step3.2. Use the high-order staggered grid finite difference formula derived in Step3.1 as the numerical solution method to solve the pure visco-acoustic TTI wave equation to obtain the difference iteration equation for any point in the simulation region. Discretize the entire new wave equation derivation formula to meet the requirements of numerical simulation, and then perform iteration to obtain a seismic wave simulation in which the seismic wave energy can rapidly decay in the PML region. The difference iteration equation for any point includes combining the spatial derivative of the velocity field and the influence of the memory variable r to obtain the iteration formula for the pressure field:

[0174]

[0175] In the formula, is the value of the pressure field at time step n + 1, and the spatial derivative of the velocity field, K is the bulk modulus.

[0176] Combining the partial derivative of the pressure field p with respect to x to obtain the iteration formula for the x-direction velocity component V x :

[0177]

[0178] In the formula, and represent the velocity components V at the same grid point and different time steps x .

[0179] Combining the partial derivative of the pressure field p with respect to z to obtain the iteration formula for the z-direction velocity component V z :

[0180]

[0181] In the formula, and represents the velocity component V at the same grid point but different time steps z .

[0182] Step4: Combine the boundary wavefield data, internal wavefield data, and corner wavefield data to obtain the complete wavefield data, and output the final wavefield snapshot of the velocity model. The wavefield snapshot is as shown in Figure 6 and perform reverse time migration using the source-normalized correlation (SNC) imaging condition. The schematic diagram of reverse time migration imaging is as shown in Figure 4 .

[0183] Obtain the pure viscoacoustic TTI wave equation through Step2, reverse the amplitude loss term to achieve attenuation compensation, and use the regularization method to suppress the high-frequency components generated by the inversion to obtain the attenuation compensation wavefield modeling operator with a regularization term:

[0184]

[0185] After determining the appropriate regularization parameter ζ, use the source-normalized correlation (SNC) imaging condition as the imaging condition, and its expression is

[0186]

[0187] where I c represents the result of the source-normalized correlation imaging condition, and α and L represent the attenuation coefficient and propagation length respectively. The subscripts up and down represent the up-going wave and down-going wave respectively. S(x,t) and R(x,t) represent the source wavefield and received wavefield in the non-attenuating anisotropic medium respectively, (x,t) represents the spatial coordinate x and time t, and S A (x,t) represents the source wavefield in the attenuating medium, R A (x,t) represents the received wavefield in the attenuating medium, R C (x,t) represents the compensated received wavefield, +at represents the variation of the attenuation compensation intensity with time, and e +αt is an exponential attenuation compensation factor used to compensate for the energy lost by the wavefield due to attenuation during propagation, and -αL down describes the variation of the attenuation intensity with the propagation length, is used to describe the energy attenuation of the down-going wave during propagation, represents the energy lost by the up-going wave due to attenuation during propagation, represents the energy recovered for the up-going wave lost due to attenuation.

[0188] As shown in Figure 7 , the general reverse time migration effect diagram shows that the general migration imaging is insufficient in correcting the amplitude attenuation and phase distortion of seismic waves, resulting in amplitude distortion and reduced resolution of the imaging result, while Figure 8The shown imaging effect of the acoustic medium migration can more accurately reflect the true structure of the underground medium after introducing the regularization parameter, and has higher imaging accuracy and resolution.

[0189] The specific implementation formula of the present invention has been described in detail above in conjunction with the accompanying drawings. However, the present invention is not limited to the above implementation formula, and various changes can be made without departing from the gist of the present invention within the scope of knowledge possessed by those of ordinary skill in the art.

Claims

1. An inverse time migration tunneling advanced detection method for viscoacoustic media with high-order staggered grids, characterized in that: The method includes: Step 1: Construct a velocity model, set the spatial and temporal sampling intervals, and then discretize the velocity model to prepare for subsequent numerical simulations. Step 2: Based on the SLS model, establish the viscoacoustic wave equation to characterize the propagation characteristics of seismic waves in an attenuating medium; combine the SLS theory with the pure qP-wave dispersion relation to derive the pure viscous acoustic wave TTI wave equation represented by memory variables. Step 3: Use the derived wave equation to simulate the propagation of seismic waves inside the velocity model to obtain the wave field data at each grid point at each moment; calculate the wave field values at the grid points on the boundary of the velocity model at each moment through the wave field absorption boundary equation and the wave field absorption coefficient, so as to obtain the boundary wave field data. Step 4: Combine the boundary wave field data, the internal wave field data, and the corner wave field data to obtain the complete wave field data, output the final wave field snapshot of the velocity model, and perform reverse time migration using the source-normalized cross-correlation imaging condition.

2. The reverse time migration tunnel advanced detection method for viscoacoustic media with a high-order staggered grid according to claim 1, characterized in that: The said Step 2 includes: Based on the constitutive relation of the SLS model, consider the motion equation and constitutive relation of the medium to derive the viscoacoustic wave equation, and combine it with the pure qP-wave dispersion relation, then introduce memory variables to simplify the solution of the acoustic wave equation, and finally obtain the pure viscous acoustic wave TTI wave equation represented by memory variables; specifically including: Step 2.1: Use a single SLS element combined with the viscoacoustic wave equation to characterize the seismic wave characteristics in an attenuating medium; for the SLS model, the constitutive equation of the SLS model is: Step 2.2: Through the dispersion relations of P-waves and SV-waves to characterize the seismic wave behavior in an anisotropic medium, set the phase velocity of the SV-wave to zero at all phase angles, add the viscous effect, and transform it to the time-space domain to derive the pure viscous acoustic wave TTI wave equation. Step 2.3: Decouple the two-dimensional pure viscous acoustic wave TTI wave equation containing the coupled amplitude attenuation and phase dispersion operators derived in Step 2.2 to derive a pure viscous acoustic wave TTI wave equation containing the decoupled amplitude attenuation and phase dispersion operators.

3. A reverse time migration tunnel advanced detection method for viscoacoustic media with a high-order staggered grid according to claim 2, characterized in that: The said Step 2.1 includes: Use a single SLS element combined with the viscoacoustic wave equation to characterize the seismic wave characteristics in an attenuating medium; for the SLS model, the constitutive equation of the SLS model is: Among them, p1, q0, and q1 are coefficients obtained through the forward and inverse Laplace transforms. σ represents stress, which is the force per unit area inside the material. represents the rate of change of stress with respect to time, that is, the derivative of stress ε represents strain, which is the amount of deformation of the material under the action of force. represents the rate of change of strain with respect to time, that is, the derivative of strain Derive the equation of the viscoacoustic wave in the frequency-wavenumber domain through the constitutive relation of the SLS model: where p is the pressure wave field, r is the memory variable, P and are the spatial Fourier transforms of p and r respectively, k x , k y , k z are the wave fields in x, y, and z, and are the stress relaxation time and the strain relaxation time respectively, where Q is the quality factor, ω represents the angular frequency, ω0 is the reference angular frequency, and in the Fourier transform, the time derivative corresponds to iω in the frequency domain, and iω represents the derivative with respect to time.

4. A reverse time migration tunnel advanced detection method for viscoacoustic media with a high-order staggered grid according to claim 2, characterized in that: The said Step 2.2 includes: Characterize the seismic wave behavior in an anisotropic medium through the dispersion relations of P-waves and SV-waves, set the phase velocity of the SV-wave to zero at all phase angles, add the viscous effect, and transform it to the time-space domain to derive the pure viscous acoustic wave TTI wave equation. In the formula, represents the variable of the pressure field p changing with time, v p0 is the P-wave velocity, respectively represent the second-order partial derivatives of the pressure with respect to the spatial coordinates x, y, and z. S is used to describe the wave propagation characteristics in the medium and is related to anisotropy and viscosity, represents the variable of the memory variable evolving with time.

5. A method for advanced detection of tunnels in viscoacoustic media by reverse time migration of high-order staggered grids according to claim 1, characterized in that: The said Step 2.3 includes: Perform coordinate transformation of VTI and TTI media, combined with the pure viscous acoustic wave TTI wave field propagation operator equation: W vis = W ac + W phase + W amp Among them, W ac is the pure acoustic wave TTI wavefield propagation operator, W phase is the phase dispersion control operator, W amp is the amplitude loss control operator; Through the addition and subtraction operations of these three operators for decoupling, derive the pure viscous acoustic wave TTI wave equation with decoupled amplitude attenuation and phase dispersion operators based on the SLS model: where k1 and k2 are the coefficients for controlling the phase dispersion term and the amplitude loss term respectively. When k1 = 0 and k2 = 1, the equation is a pure visco-acoustic TTI wave equation dominated by phase dispersion. When k1 = 1 and k2 = 0, the equation is a pure visco-acoustic TTI wave equation dominated by amplitude loss. When k1 = 1 and k2 = 1, the equation includes both the phase dispersion and amplitude loss terms, and can comprehensively describe the propagation characteristics of visco-acoustic waves in TTI media.

6. The inverse time migration tunnel advanced detection method for viscoacoustic media with a high-order staggered grid according to claim 1, wherein: The said Step3 includes: Carry out forward modeling on the new wave equation and integrate the perfectly matched layer. Use the high-order staggered grid finite difference method as the numerical solution method to solve the pure visco-acoustic TTI wave equation to obtain the difference iteration equation at any point in the simulation area. The specific steps include: Step3.1: According to the one-dimensional second-order visco-acoustic medium wave equation derived in Step2, determine the staggered grid finite difference format. Expand the computational variable matrix by a certain number of layers around as the perfectly matched layer for absorbing seismic waves. The main area is the computational interval and also the target area of the simulation. The simulation object is the tunnel isomer. The expanded range is the perfectly matched layer PML area; Step3.2: Use the high-order staggered grid finite difference formula derived in Step3.1 as the numerical solution method to solve the pure visco-acoustic TTI wave equation to obtain the difference iteration equation at any point in the simulation area. Discretize the entire new wave equation acquisition formula to meet the requirements of numerical simulation, and then perform iteration to obtain a seismic wave simulation in which the seismic wave energy can rapidly decay in the PML area.

7. A reverse time migration tunnel advanced detection method for visco-acoustic media with a high-order staggered grid according to claim 6, characterized in that: The said Step3.1 includes: The PML area is divided into three sub-areas: the upper and lower boundary areas with the x-axis as the main direction, the left and right boundary areas with the z-axis as the main direction, and the corner intervals at the four corners; each sub-area will independently apply the attenuation factor α according to its own characteristics and requirements to achieve accurate and efficient seismic wave absorption; Among them, α(Z) represents the attenuation factor in the PML region, which is used to describe the attenuation intensity of seismic waves in the PML region. Z represents the distance from the boundary of the PML region, and β is a parameter that controls the shape of the attenuation factor and is used to adjust the smoothness of the attenuation curve. Make the function meet the basic requirements of the attenuation function. Fine-tune the function, K, γ, δ are adjustable coefficients, and K is a known number, n is the order, L is the layer thickness of the PML; R is the theoretical boundary reflection coefficient; ρ is the medium density, and u is the elastic modulus of the medium. Assign values to the frequency shift factor η x and the scale factor β x in the PML region. In the main region, the frequency shift factor η x is set to 0, and the scale factor β x is 1. The calculation formula for the scale factor β x is as follows; where, β0 represents the reference value of the scaling factor β x and η0 represents the reference value of the frequency shift factor η x ; x is the distance measured from the PML region boundary, P η and P β are the power exponents that control the variation rate of η x and β x ; and are the normalized power functions that respectively describe the variation laws of η x and β x with the position x, and f is the function used to adjust the distribution of η x . For the second-order acoustic wave equation, the wave field function is set as Let Δx, Δz, and Δt denote the spatial step sizes in the x and z directions and the time step size, respectively. Let i, j, and k be the grid numbers in space and time. Then we have Using the Taylor formula to expand p at time t, we get and Derive the coefficient matrix equation: Solving the difference coefficient C n , select the absorption function and substitute the Ricker wavelet as the boundary condition; the expression formula of the absorption function is as follows: Among them, L is the total number of layers in the absorption region, k is the number of layers from the forward modeling region, d(k) is the absorption function used to describe the absorption intensity of the k-th layer in the PML region, N is the total number of layers in the absorption region, V p is the P-wave velocity, and R = 0.0001; The expression formula of the Ricker wavelet is: where r(t) represents the amplitude value of the Ricker wavelet at time t, and f m is the dominant frequency of the Ricker wavelet; Solve the difference coefficients of the first-order derivative with 2N-order accuracy. The high-order accuracy operator of the first-order derivative is expressed as: In the formula, and respectively represent the pressure values that are offset by and grid units in the x direction from the grid point i; The difference formula is: where, f (2) (x) represents the second-order derivative of f(x) in the x direction, a0 is the center point coefficient, f(x0) represents the value of f(x) at x0, f is the velocity and stress component to be differentiated; Δx, Δz are the spatial step sizes in the x and z directions; t is the simulation time; Δt is the time step; a m is the coefficient for the high-order accuracy difference approximation of the second-order derivative, m and M are positive integers from one to infinity, x0 is the unknown in the process of partial derivative, represents the change in the function value within the time step Δt.

8. A reverse time migration tunnel advanced detection method for viscoacoustic media with a high-order staggered grid according to claim 6, characterized in that: The said Step3.2 includes: The difference iteration equation at any point includes obtaining the iteration formula of the pressure field by combining the spatial derivative of the velocity field and the influence of the memory variable r; where, is the value of the pressure field at time step n + 1, and the spatial derivative of the velocity field, and K is the bulk modulus; Combined with the partial derivative of the pressure field p with respect to x The velocity component V in the x-direction is obtained x The iterative formula for: Among them, and represent the velocity components V at the same grid point but different time steps x ; Combined with the partial derivative of the pressure field p with respect to z The z-direction velocity component V is obtained z The iterative formula for: Among them, and represent the velocity components V at the same grid point but different time steps z .

9. A reverse time migration tunnel advanced detection method for a viscoacoustic medium with a high-order staggered grid according to claim 1, characterized in that: The said Step4 includes: Through the pure viscous acoustic wave TTI wave equation obtained in Step2, reverse the amplitude loss term to achieve attenuation compensation, and use the regularization method to suppress the high-frequency components generated by the inversion to obtain an attenuation compensation wavefield modeling operator with a regularization term; After determining the appropriate regularization parameter ζ, use the source-normalized correlation imaging condition as the imaging condition. The expression formula of the source-normalized correlation imaging condition is: Among them, I c represents the result of the source-normalized correlation imaging condition, where α and L represent the attenuation coefficient and the propagation length respectively; the subscripts up and down represent the up-going wave and the down-going wave respectively; S(x,t) and R(x,t) represent the source wave field and the received wave field in the non-attenuating anisotropic medium respectively, and (x,t) represents the spatial coordinate x and the time t. S A (x,t) represents the source wave field in the attenuating medium, and R A (x,t) represents the received wave field in the attenuating medium, and R C (x,t) represents the compensated received wave field, +αt represents the variation of the attenuation compensation intensity with time, and e +αt is an exponential attenuation compensation factor used to compensate for the energy lost by the wave field due to attenuation during propagation, and -αL down describes the variation of the attenuation intensity with the propagation length, is used to describe the energy attenuation of the down-going wave during propagation, represents the energy lost by the up-going wave due to attenuation during propagation, represents the energy recovered for the up-going wave lost due to attenuation.

10. An inverse time migration tunnel advanced detection system for viscoacoustic media with a high-order staggered grid, characterized in that, The said system includes: a module for executing a visco-acoustic medium reverse time migration tunnel advance detection method of a high-order staggered grid as described in any one of claims 1 to 9.

Images