A seismic wave forward modeling and two-way wave illumination method for two-phase viscoelastic TTI media

Abstract

The present invention discloses a seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium, which relates to the field of petroleum geophysical exploration technology, including: introducing the Biot two-phase medium theory to describe the two-phase medium; deriving the constitutive equation of the TTI medium; obtaining the stress-strain relationship of the solid skeleton; deriving the wave equation applicable to the two-phase viscoelastic TTI medium; using the high-order staggered grid finite difference method to perform seismic wave forward modeling on the derived first-order velocity stress equation of the two-phase viscoelastic TTI medium; establishing an illumination operator, defining source and reverse illumination operators, and calculating the source illumination through forward modeling; based on the reciprocal principle, analyzing the reverse illumination results, and optimizing the excitation position and shot point density. The present invention combines TTI medium and Biot two-phase viscoelastic medium for the first time, realizing effective exploration illumination of complex two-phase media with viscoelasticity and TTI anisotropy such as oil and gas formations, and making up for the current lack of knowledge of the propagation characteristics of seismic waves in such media.

Description

Technical Field

[0001] The invention relates to the technical field of petroleum geophysical exploration, and in particular to a seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium. Background Art

[0002] In oil and gas exploration, oil and gas-bearing strata, due to their unique stratification and directional development of fractures, coexist with solid and fluid phases, forming a two-phase complex medium with viscoelasticity and inclined transverse anisotropy (TTI).

[0003] Therefore, for the needs of oil and gas exploration, it is particularly important to understand the propagation characteristics of seismic waves in this complex medium, and it is urgent to study a seismic wave forward simulation and two-way wave illumination method for two-phase viscoelastic TTI medium. Summary of the invention

[0004] In order to achieve forward simulation of complex oil and gas media and describe underground media more realistically and accurately, the present invention discloses a seismic wave forward simulation and two-way wave illumination method for a two-phase viscoelastic TTI medium.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

[0006] A seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium comprises the following steps:

[0007] Step S1, introducing the Biot two-phase medium theory to describe the two-phase medium, and using the coupling coefficient to represent the interaction between the volume change of the solid phase and the liquid phase;

[0008] Step S2, considering the anisotropy of the two-phase medium, taking the TTI medium as the solid skeleton, and deriving the constitutive equation of the TTI medium;

[0009] Step S3, introduce Carcione's viscoelastic anisotropic medium model and combine it with the Boltzmann superposition principle to obtain the stress-strain relationship of the solid skeleton:

[0010] Step S4, embedding the viscoelastic constitutive relation of the solid phase into Biot's two-phase medium theory, and deriving a wave equation applicable to the two-phase viscoelastic TTI medium;

[0011] Step S5, using a high-order staggered grid finite difference method to perform seismic wave forward simulation on the derived first-order velocity stress equation of the two-phase viscoelastic TTI medium;

[0012] Step S6, establishing an illumination operator, defining source direction and reverse illumination operators, calculating source direction illumination through forward simulation, and analyzing the illumination effects on the target area under different configurations;

[0013] Step S7: Based on the reciprocity principle, analyze the reverse illumination results, identify the surface locations that contribute more to the target area, and optimize the excitation positions and shot point density to achieve optimal coverage.

[0014] Optionally, in step S2, the constitutive equation of the TTI medium is:

[0015] T TTI =C TTI S TTI ;

[0016] Where, T TTI is the TTI medium stress vector, S TTI is the TTI medium strain vector.

[0017] Optionally, in step S3, the stress-strain relationship of the solid skeleton is:

[0018]

[0019] In the formula, T=(σ xx ,σ yy ,σ zz ,σ yz ,σ zx ,σ xy ) is the stress vector, is the first-order time derivative of the strain vector, ψ represents the relaxation matrix;

[0020] Separate the integral and non-integral terms in the equation and substitute them into the equation T TTI =C TTI S TTI ,get:

[0021]

[0022] In the formula, H TTI The elastic matrix representing the non-integral term, H TTI v represents the viscoelastic matrix of the integral term, represents the memory variable, v = 1, 2 represents the expansion and contraction motion and shear motion of the particles in the medium, L v Represents the number of units that make up the generalized standard linear device.

[0023] Optionally, in step S4, the wave equation applicable to the two-phase viscoelastic TTI medium is:

[0024]

[0025] In the formula, U=(u x ,u y ,u z ,U x ,Uy ,U z ) is the displacement vector; represents the second-order derivative of U with respect to time; F = (f x ,f y ,f z ,F x ,F y ,F z ) is the body force vector, where f x 、f y 、f z represents the solid phase fraction, F x 、F y 、F z represents the liquid phase fraction, Compared to A 7×7 matrix with a row vector and a column vector with 0 elements is added respectively; ρ is the density matrix; G is the partial differential operator; Β is the dissipation coefficient matrix; J is the sum and difference operator.

[0026] Optionally, in step S5, the first-order time derivatives of the stress and strain of the wave equation of the two-phase viscoelastic TTI medium are obtained to obtain the first-order velocity stress equation of the two-phase viscoelastic TTI medium; in the two-dimensional case, the velocity stress equation is:

[0027]

[0028]

[0029] In the formula, ρ 11 is the solid effective mass of the medium volume element; ρ 12 is the mass coupling coefficient between the solid phase and the fluid phase, also called coupling mass; ρ 22 is the effective mass of the liquid phase per volume of the medium; v x and v z are the velocities of the solid phase along the x and z directions respectively; b ij ,i,j=1,2,3 represents the dissipation coefficient; V x and V z are the velocities of the liquid phase along the x and z directions, respectively;

[0030] The first-order time derivative equation of the memory variable is:

[0031]

[0032] In the formula, the relaxation time τ can be calculated by the quality factor, and the calculation formula is:

[0033]

[0034] Where v = 1, 2, Q1 and Q2 are the longitudinal wave quality factor and the shear wave quality factor, and ω is the reciprocal of the main frequency of the source wavelet.

[0035] Optionally, in step S6, the energy and illumination are calculated based on the two-way wave forward modeling, the wave field value is recorded as U, and the shadow area pseudo-reverse excitation point in the observation system source illumination (also called forward illumination) result is marked as X. o (x o ,y o ,z o ), the area where the surface observation system is located is recorded as D, and the position coordinates within D are generally expressed as X g (x g ,y g ,0), the source excitation point of the source illumination is recorded as X s (x s ,y s ,0), then the illumination of the observation system to the reverse excitation point is:

[0036]

[0037] The illumination function of the reversely excited wave field on the surface observation system is:

[0038]

[0039] According to the coordinate point X corresponding to the function extreme value of the above formula m (x m ,y m ,0), by function I g (X g ) indicates that the coordinate point X m (x m ,y m ,0) is the maximum value in the surface observation system.

[0040] Optionally, in step S7, the encrypted surface observation system source excitation point is recorded as X′ s (x′ s ,y′ s ,0), the coordinates of any point in the underground medium model area Ω are marked as X(x,y,z), then the illumination function of the optimized observation system for the medium model is:

[0041]

[0042] The above formula is the optimized lighting function for lighting enhancement in shadow areas after two-way wave lighting analysis, i.e., the lighting operator.

[0043] According to the derivation logic of the equilibrium equation of motion of two-phase medium, it is not difficult to see that the energy of two-phase medium is composed of the energy of the solid phase and liquid phase, as well as the energy lost by the solid-liquid coupling effect. If the x component and z component of the solid phase wave field of the two-phase medium are denoted as U x and U z , the x-component and z-component of the liquid phase wave field are V x and V z , then the solid phase energy is:

[0044]

[0045] The liquid phase energy is:

[0046]

[0047] The energy of the coupling effect is:

[0048] E d =ρ 12 (U x V x +U z V z );

[0049] So the total energy of the two-phase medium is:

[0050]

[0051] In two-dimensional illumination, if the coordinates of each position in the surface observation range are generally expressed as X g (x g ,0), the excitation positions of the observation system are X s (x s ,0), the coordinates of any point D in the two-phase viscoelastic medium model area are X(x,z), then according to the principle of statistical illumination method, the source-directed illumination function of the observation system for the two-phase medium is:

[0052]

[0053] I(X) can be used as the source-direction illumination operator of the dual-phase viscoelastic TTI medium of the present invention;

[0054] The target point of the source illumination in the shadow area is X o (x o ,z o ), the shadow area is recorded as D′, and several excitation points are arranged in the shadow area for reverse illumination. The statistically obtained illumination function within the observation range of the surface is:

[0055]

[0056] I g (X g) can be used as the reverse illumination operator of the two-phase viscoelastic TTI medium of the present invention;

[0057] Remember I g (X g ) has an extreme point coordinate of X m (x m ,0),According to the lighting optimization scheme based on the reciprocal principle,X m (x m ,0) contributes the most to the illumination of the shadow area, so the original observation system's shot point excitation position is adjusted to m (x m ,0) to encrypt the shot points and optimize the observation system. The encrypted surface observation system is recorded as X′ s (x′ s ,0). The illumination function of the optimized and encrypted observation system for the two-phase medium is:

[0058]

[0059] The beneficial effect of the present invention is that the present invention is based on the Biot two-phase medium theory, introduces the Carcione system and combines the constitutive relationship of the TTI medium to describe the viscoelastic anisotropy of the solid skeleton. By embedding the constitutive relationship of the solid phase into the Biot two-phase system, the wave equation of the two-phase viscoelastic TTI medium is derived; then, the first-order velocity-stress equation of the two-phase viscoelastic TTI medium is simulated by seismic wave forward modeling using a high-order staggered grid finite difference method, and a two-phase illumination operator is established using the derived wave equation, thereby performing two-way wave illumination on the two-phase viscoelastic TTI medium; finally, the wave field characteristics of the two-phase viscoelastic TTI medium are analyzed through numerical simulation results. The present invention realizes effective exploration illumination of complex two-phase media with viscoelasticity and TTI anisotropy, such as oil and gas formations, and makes up for the current lack of understanding of the propagation characteristics of seismic waves in such media. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1 It is a flow chart of a seismic wave forward modeling and two-way wave illumination method of a two-phase viscoelastic TTI medium according to the present invention;

[0061] Figure 2 It is the wave field degradation diagram of the two-phase viscoelastic TTI medium;

[0062] Figure 3 The wave field degradation diagram of the two-phase viscoelastic TTI medium after adjusting the order;

[0063] Figure 4 Establish the two-phase viscoelastic TTI medium model parameters for the depression model construction of a certain work area;

[0064] Figure 5AThis is the result diagram of the solid phase operator source illumination of the depression model in the work area;

[0065] Figure 5B This is the result diagram of the source direction illumination of the liquid phase operator of the depression model in this work area;

[0066] Figure 5C This is the source-direction illumination result diagram of the two-phase operator of the depression model in this work area;

[0067] Fig. 6A This is the source-direction illumination diagram of the depression model in the work area;

[0068] Figure 6B This is the reverse lighting diagram of the depression model in the work area;

[0069] Figure 7 This is the normalized contour map of the energy below 400m depth from the source direction of the depression model in the work area;

[0070] Figure 8 Back-illuminate the surface energy curve for the depression model of the work area;

[0071] Fig. 9 The lighting results after optimizing and encrypting the sunken model of the work area;

[0072] Fig.10 Normalize the energy below 400m depth before and after the sag model of the work area is densified;

[0073] Fig.11 The distribution of the parameter values ​​of the Marmours model of two-phase viscoelastic TTI properties;

[0074] Fig. 12A The source-directed illumination result of the Marmours model of the two-phase viscoelastic TTI property obtained by the method of the present invention;

[0075] Fig. 12B This is the energy distribution diagram of the Marmours model below 1200m depth;

[0076] Fig.13A The reverse illumination result of the Marmours model of the two-phase viscoelastic TTI property obtained by the method of the present invention;

[0077] Fig. 13B is the normalized curve of reverse illumination surface energy;

[0078] Fig.14A This is the source illumination result of the Marmours model after shot point optimization and encryption;

[0079] Fig. 14B Energy normalization for the encrypted portion below 1200m for the Marmours model;

[0080] Fig.15 Comparison of normalized energy before and after optimization of the Marmous model below 1200m depth. DETAILED DESCRIPTION

[0081] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the invention claimed for protection, but merely represents the selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0082] The present invention proposes a seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium. Figure 1 As shown, the specific steps include:

[0083] Step S1, introducing the Biot two-phase medium theory to describe the two-phase medium; in the Biot theory, the geometric equation of the two-phase medium (i.e., the Cauchy equation) is:

[0084]

[0085] In the formula, u=(u x ,u y ,u z ) represents the solid phase displacement of the medium, and U = (U x ,U y ,U z ) represents the liquid phase displacement of the medium, e represents the solid phase strain, ε represents the liquid phase strain,

[0086] The general form of the constitutive equation for a two-phase anisotropic medium is:

[0087]

[0088] In the formula, C represents the solid elastic matrix, and the matrix element is c ij ,i,j=1,2,…6;T=(σ xx ,σ yy ,σ zz ,σ yz ,σ xz ,σxy ) T and T′=(s) T are the solid and liquid phase force vectors, S = (e xx ,e yy ,e zz ,e yz ,e xz ,e xy ) T and S′=(ε) T are the solid and liquid phase deformation vectors respectively; P = (P1, P2, P3, P4, P5, P6) T represents the coupling coefficient vector, R = (R) T represents the liquid phase pressure coefficient vector; P i ,i=1,2,…,6 is the coupling coefficient, which represents the interaction between the volume changes of the solid phase and the liquid phase; R is the liquid phase pressure coefficient, which represents a measure of the force required on the fluid to force a certain volume of fluid to flow into the volume element of the two-phase medium while keeping the volume constant.

[0089] Step S2: Using the TTI medium as a solid skeleton, deriving the constitutive equation of the TTI medium.

[0090] Based on the anisotropy of the two-phase medium, taking the TTI medium as the solid skeleton, the constitutive equation of the TTI medium is first derived. The derivation process of the constitutive equation of the TTI medium is as follows:

[0091] Assume that the constitutive equation of VTI medium is:

[0092] T VTI =C VTI S VTI (3);

[0093] Where, T VTI is the VTI medium stress vector, S VTI is the strain vector of the VTI medium, C VTI is the elastic matrix of the VTI medium; according to the Bond transformation, the elastic matrix is ​​transformed from the constitutive coordinate system to the observation coordinate system to obtain the elastic matrix of the TTI medium:

[0094]

[0095] Where, θ is the polarization angle of the TI medium; is the azimuth angle of the TI medium; M θ is the polarization anisotropy matrix, is the azimuthal anisotropy matrix;

[0096] Substituting equation (4) into equation (3), the constitutive equation of TTI medium is obtained as follows:

[0097] T TTI =C TTI S TTI (5);

[0098] Where, T TTI is the TTI medium stress vector, S TTI is the TTI medium strain vector.

[0099] Step S3: Carcione's viscoelastic anisotropic medium model is introduced, and combined with the Boltzmann superposition principle, the stress-strain relationship of the solid skeleton is obtained.

[0100] Carcione introduced the superposition principle of the cumulative integral method and the anisotropy of the medium to obtain a viscoelastic anisotropic medium model based on a generalized standard linear body. According to Carcione's viscoelastic anisotropy theory and Boltzmann's superposition principle, the stress-strain relationship of the solid skeleton is:

[0101]

[0102] In the formula, T=(σ xx ,σ yy ,σ zz ,σ yz ,σ zx ,σ xy ) is the stress vector, is the first-order time derivative of the strain vector, and ψ represents the relaxation matrix.

[0103] Step S4: embed the viscoelastic constitutive relation of the solid phase into Biot's two-phase medium theory to derive the wave equation applicable to the two-phase viscoelastic TTI medium, which specifically includes:

[0104] Separating the integral term and the non-integral term in equation (6) and substituting them into equation (5), we get:

[0105]

[0106] In the formula, H TTI The elastic matrix representing the non-integral term, H TTI v represents the viscoelastic matrix of the integral term, represents the memory variable, v = 1, 2 represents the expansion and contraction motion and shear motion of the particles in the medium, L v Indicates the number of units that make up the generalized standard linear device;

[0107] By inserting the viscoelastic TTI solid skeleton into the mechanism of the two-phase medium, ignoring the viscosity of the fluid, and replacing the solid term in equation (2) with equation (7), the constitutive equation of the two-phase viscoelastic TTI medium based on the standard linear viscoelastic body is obtained:

[0108]

[0109] Where v = 1, 2, l = 1, ..., L v Without considering the viscosity of the fluid, P = (P1, P2, P3, 0, 0, 0) T ;

[0110] Using the principles of Newtonian mechanics to analyze the dynamic characteristics of a two-phase medium, take any volume element of the medium, and the kinetic energy of the unit cube is:

[0111]

[0112] In the formula, ρ 11 is the solid effective mass of the medium volume element; ρ 12 is the mass coupling coefficient between the solid phase and the fluid phase, also called coupling mass; ρ 22 is the effective mass of the liquid phase per volume element of the medium;

[0113] Due to the interaction between the solid phase and the liquid phase, there is energy loss per unit volume of the two-phase medium relative to the single-phase medium, so the dissipation function is defined as:

[0114]

[0115] Where B is the dissipation coefficient, satisfying η is the fluid viscosity coefficient, k is the permeability;

[0116] After considering the energy dissipation of the solid and liquid phases, the energy of the volume element is:

[0117] E k ′=E k -D(11);

[0118] The momentum of the solid phase p k =(p x ,p y ,p z )for:

[0119]

[0120] Momentum of the liquid phase P k =(P x ,P y ,P z )for:

[0121]

[0122] The force on the solid phase in the x direction is:

[0123]

[0124] The force on the liquid phase in the x direction is:

[0125]

[0126] According to the relationship between kinetic energy and momentum, the relationship between momentum and force, and equation (11), the force in the solid phase in the x direction is obtained as:

[0127]

[0128] The force on the liquid phase in the x direction is:

[0129]

[0130] Considering the anisotropy of the medium, the motion of the solid and liquid phases satisfies the generalized Darcy's law, and the dissipation coefficient is a symmetric matrix B, expressed as:

[0131]

[0132] The stress-displacement relationship does not consider viscoelasticity. The simultaneous equations of stress and momentum are used to express the motion equilibrium equation of the two-phase viscoelastic TTI medium. The stress-displacement equation of the solid phase is:

[0133]

[0134] The stress-displacement equation of the liquid phase is:

[0135]

[0136] Substituting equation (1) and constitutive equation (8) into the kinematic equilibrium equations (19) and (20), the wave equation of the two-phase viscoelastic TTI medium expressed in displacement form is obtained:

[0137]

[0138] In the formula, U=(u x ,u y ,u z ,U x ,U y ,U z ) is the displacement vector; represents the second-order derivative of U with respect to time; F = (f x ,f y ,f z ,F x ,F y ,F z ) is the body force vector, where f x 、f y 、f z represents the solid phase fraction, F x 、Fy 、F z represents the liquid phase fraction, Compared to A 7×7 matrix with a row vector and a column vector with 0 elements is added respectively; ρ is the density matrix; G is the partial differential operator; Β is the dissipation coefficient matrix; J is the sum and difference operator;

[0139] In formula (21), the correlation coefficient, operator matrix and calculation formula are:

[0140]

[0141]

[0142]

[0143] In the formula,

[0144] M represents the number of viscoelastic units in the solid phase, and N represents the number of viscoelastic units in the liquid phase. The values ​​of M and N determine the relaxation function and viscoelastic properties in the solid phase and liquid phase, respectively. ij ,i,j=1,…,6 represents the elastic coefficient, χ v represents the relaxation function, v = 1 and v = 2 correspond to the expansion strain and shear strain, that is, the relaxation function of the longitudinal wave and the transverse wave, and They represent the strain relaxation time and stress relaxation time of the first device assembly in the standard linear body respectively; τ0 is the reciprocal of the center frequency, Q v is the longitudinal or transverse wave quality factor corresponding to the center frequency, L v Corresponding to the number of device assembly units in the viscoelastic model.

[0145] Step S5, using a high-order staggered grid finite difference method to perform seismic wave forward simulation on the derived first-order velocity stress equation of the two-phase viscoelastic TTI medium; taking the first-order time derivatives of the stress and strain of the wave equation of the two-phase viscoelastic TTI medium to obtain the first-order velocity stress equation of the two-phase viscoelastic TTI medium.

[0146] In two dimensions, the velocity stress equation is:

[0147]

[0148] In formula (22), ρ 11 is the solid effective mass of the medium volume element; ρ 12 is the mass coupling coefficient between the solid phase and the fluid phase, also called coupling mass; ρ 22 is the effective mass of the liquid phase per volume of the medium; vx and v z are the velocities of the solid phase along the x and z directions respectively; b ij ,i,j=1,2,3 represents the dissipation coefficient; V x and V z are the velocities of the liquid phase along the x direction and the z direction respectively; in formula (23),

[0149] The first-order time derivative equation of the memory variable is:

[0150]

[0151] In the formula, the relaxation time τ is calculated by the quality factor, and the calculation formula is:

[0152]

[0153] Where v = 1, 2, Q1 and Q2 are the longitudinal wave quality factor and the shear wave quality factor, and ω is the reciprocal of the main frequency of the source wavelet;

[0154] The staggered grid finite difference method is used to differentiate equations (22), (23) and (24) to obtain iterative formulas for calculating velocity, stress and memory variables. The perfectly matched layer (PML) method is used in the forward modeling to decompose the velocity and stress into two components parallel and perpendicular to the artificial boundary, and substituted them into the difference equation. By introducing the attenuation factor, the seismic wave is strongly attenuated in the artificially established attenuation layer to eliminate boundary reflections.

[0155] In order to make the difference stable, relevant stability conditions need to be met. Taking into account the viscoelasticity and TI anisotropy of the medium, the necessary stability conditions for the forward simulation of the first-order velocity stress equation of the two-phase viscoelastic TTI medium using the high-order staggered grid finite difference method are:

[0156]

[0157] Where Δt is the time sampling interval, Δx and Δz are the spatial steps in the x-direction and z-direction, i.e., the grid spacing. is the 2N-order spatial staggered grid finite difference coefficient;

[0158] Using the Routh criterion, the stability condition that needs to be satisfied by the staggered grid finite difference forward modeling of viscoelastic media is derived as follows:

[0159]

[0160] By analyzing the stability of the staggered grid finite difference method for viscoelastic TTI media, the stability condition is obtained as follows:

[0161]

[0162] In the formula,

[0163] Step S6: Establish an illumination operator, define source direction and reverse illumination operators, calculate source direction illumination by forward simulation, and analyze the illumination effect on the target area under different configurations; use the introduced two-phase viscoelastic TTI medium wave equation to perform illumination work, the process is as follows:

[0164] Assume that seismic waves propagate in the underground space, and any particle in a certain area Ω underground is denoted as X(x, y, z), and the wave field of the particle is denoted as U. The observation system has m excitation points, and each source is denoted by S1, S2, ..., S m According to the relationship between seismic wave energy and wave field value, the particle X0 (x0, y0, z0) at t i The energy at that moment is recorded as Assume that the total number of time sampling points in the whole process of seismic wave propagation in the forward simulation is n, then the single source S j The illumination of the particle is:

[0165]

[0166] The illumination of the particle by the entire observation system is:

[0167]

[0168] The illumination of all particles in the model area D by the observation system is statistically analyzed to obtain the illumination result of the medium model by the observation system:

[0169]

[0170] Formula (31) is the expression of the illumination function, also known as the illumination operator of the observation system on the medium model. The illumination process of the observation system excited from the surface to the underground medium is called source-directed illumination.

[0171] Based on the two-way wave forward modeling, the energy and illumination are calculated. The wave field value is denoted as U, and the shadow area of ​​the pseudo-reverse excitation point in the observation system source illumination (also called forward illumination) result is marked as X. o (x o ,y o ,z o ), the area where the surface observation system is located is denoted as D, and the position coordinates within D are denoted as X g (x g ,y g ,0), the source excitation point of the source illumination is recorded as X s (x s ,y s ,0), then the illumination of the observation system to the reverse excitation point is:

[0172]

[0173] The illumination function of the reversely excited wave field on the surface observation system is:

[0174]

[0175] According to the function extreme value of formula (33), the coordinate point X m (x m ,y m ,0), by function I g (X g ) indicates that the coordinate point X m (x m ,y m ,0) is the maximum value in the surface observation system.

[0176] Step S7: Based on the reciprocity principle, analyze the reverse illumination results, identify the surface positions that contribute more to the target area, and optimize the excitation positions and shot point density to achieve the best coverage. Based on the reciprocity principle followed by seismic wave propagation, reverse illumination is performed with the target area as the source. The energy results on the surface can indicate the most favorable excitation position for improving the wave field energy in the target area, thereby guiding the adjustment of the excitation position and the number of shot points of the observation system, thereby optimizing the observation system.

[0177] According to the inverse principle, the seismic wave excited by the coordinate point with the extreme value in equation (33) affects the target point X o (x o ,y o ,z o ) has the largest lighting energy, that is, if m (x m ,y m ,0) can enhance the illumination contribution to the shadow area; the encrypted surface observation system source excitation point is denoted as X′ s (x′ s ,y′ s ,0), the coordinates of any point in the underground medium model area Ω are marked as X(x,y,z), then the illumination function of the optimized observation system for the medium model is:

[0178]

[0179] Formula (34) is the optimized lighting function for lighting enhancement in shadow areas after two-way wave lighting analysis, i.e., the lighting operator;

[0180] The energy of the two-phase medium is composed of the energy of the solid phase and the liquid phase, as well as the energy lost due to the solid-liquid coupling effect. The x component and z component of the solid phase wave field of the two-phase medium are denoted as U x and U z , the x-component and z-component of the liquid phase wave field are V x and V z , the solid phase energy is:

[0181]

[0182] The liquid phase energy is:

[0183]

[0184] The energy of the coupling effect is:

[0185] E d =ρ 12 (U x V x +U z V z )(37);

[0186] The total energy of the two-phase medium is:

[0187]

[0188] In two-dimensional lighting, the coordinates of each position in the surface observation range are recorded as X g (x g ,0), the excitation positions of the observation system are X s (x s ,0), the coordinates of any point D in the two-phase viscoelastic medium model area are X(x,z), then according to the principle of statistical illumination method, the source-directed illumination function of the observation system for the two-phase medium is:

[0189]

[0190] In formula (39), I(X) is used as the source-direction illumination operator of the two-phase viscoelastic TTI medium of the present invention;

[0191] The target point of the source illumination in the shadow area is X o (x o ,z o ), the shadow area is recorded as D′, and several excitation points are arranged in the shadow area for reverse illumination. The statistically obtained illumination function within the observation range of the surface is:

[0192]

[0193] In formula (40), I g (X g) as the reverse illumination operator of the two-phase viscoelastic TTI medium of the present invention;

[0194] Remember I g (X g ) has an extreme point coordinate of X m (x m ,0),According to the lighting optimization scheme based on the reciprocal principle,X m (x m ,0) contributes the most to the illumination of the shadow area, so the original observation system's shot point excitation position is adjusted to m (x m ,0) to encrypt the shot points and optimize the observation system. The encrypted surface observation system is recorded as X′ s (x′ s ,0); the illumination function of the optimized and encrypted observation system for the two-phase medium is:

[0195]

[0196] Adjust the excitation position and shot point density in the observation system, especially in the shadow area where the lighting effect is originally weak. Through the illumination analysis of reverse excitation, select the energy maximum point for encryption, optimize the design of seismic exploration, and ensure the best coverage of seismic wave energy in the target area.

[0197] Application Examples

[0198] (a) Wave field numerical simulation analysis based on parameter changes

[0199] A uniform medium model with a grid size of 301×301 was established, and a Ricker wavelet longitudinal wave source was used to excite at the center of the model (151,151). The source frequency was 25Hz, the time sampling interval was 0.5ms, and the spatial sampling interval was Δx=Δz=10m. The verification work started with the forward modeling of the two-phase viscoelastic medium model. After that, each subsequent stage of the forward modeling took the relevant parameters as zero based on the previous step. The wave equation derived in this application example was verified and analyzed based on the forward modeling results of each step. The initial parameters of the two-phase viscoelastic TTI medium model are shown in Table 1.

[0200] Table 1 Parameters of the uniform model of two-phase viscoelastic TTI medium

[0201]

[0202] Note: P0 、v S0 The unit is m / s; ρ 11 , 12 , 22 The unit is kg / m 3 ; The unit of P1 and R is 10 9kg·m -1 ·s -2 ; b 11 The unit is kg·m -3 ·s -1 ; The remaining parameters are dimensionless.

[0203] like Figure 2 As shown, the wave field snapshots at each step in the process of the two-phase viscoelastic TTI medium first degenerating into various types and simple two-phase media, and finally degenerating into a single-phase homogeneous isotropic perfectly elastic medium. The specific steps are as follows: from left to right, first, let the polarization angle θ = 0, the energy of the wave field of the two-phase medium does not change much, the wavefront morphology of the longitudinal and transverse waves changes, the symmetry axis rotates to the direction parallel to the z-axis, and has the anisotropy and transverse wave splitting characteristics of the fast longitudinal wave. At the same time, there are also slow longitudinal waves. Its wavefield characteristics are consistent with the wavefield characteristics of the two-phase viscoelastic VTI medium, indicating that when the polarization angle is zero, the equation degenerates into the wave equation of the two-phase viscoelastic VTI medium; secondly, on the basis of θ = 0, let the anisotropy parameter ε = γ = δ = 0, the total energy of the wavefield still does not change much, but the wavefront morphology has changed greatly, the fast longitudinal wave wavefront changes from the original distance excitation point position to the direction axis symmetric anisotropy, and evolves into a regular circle, while the transverse wave no longer appears, and the slow longitudinal wave characteristics remain basically unchanged, which is consistent with the wavefield characteristics of the two-phase isotropic medium. Combined with the characteristics of the energy change before and after, it can be judged that on this basis, the anisotropy parameter is taken to zero, and the medium degenerates into a two-phase viscoelastic isotropic medium. Again, let the longitudinal and transverse wave quality factors Q P and Q S is infinite, the wavefront shape in the wave field is consistent, but the energy is enhanced, the longitudinal and transverse wave quality factors are infinite, and the seismic wave energy attenuation phenomenon caused by the viscosity of the medium no longer occurs, which conforms to the wave field characteristics of a two-phase isotropic perfect elastic medium, and the medium degenerates into a two-phase isotropic perfect elastic medium. Finally, let the relevant parameter ρ that describes the two-phase medium be 12 =ρ 22 =P i =R = b ij =0, the wave field energy is further enhanced, and the slow longitudinal wave no longer exists. This is because the parameter describing the interaction between the solid and liquid phases of the two-phase medium is zero and no longer plays a role in the wave equation. The forward modeling result of the medium no longer has the characteristics of a two-phase medium, but evolves into a uniform single-phase isotropic perfectly elastic medium.

[0204] Figure 3 is relative Figure 2 Numerical simulations of sequence-modified degradation were performed to validate the forward modeling results, and Figure 2The forward modeling series uses the same set of parameter tables. The idea of ​​this forward modeling is to first degenerate the two-phase viscoelastic TTI medium into a single-phase viscoelastic TTI medium, and then further eliminate the anisotropy and viscoelasticity of the medium, and finally degenerate into a single-phase homogeneous isotropic perfect elastic medium. The specific steps are as follows: from left to right, first assign the parameters describing the two-phase properties of the medium to ρ 12 =ρ 22 =P i =R = b ij =0, the wave field snapshot has no slow P-wave component, while other components are retained. The wave field morphology of fast P-wave and S-wave has not changed, and their symmetry is consistent, but the energy has been significantly enhanced. This is because there is no effect of the two-phase parameter, the medium does not show biphasicity, and there is no energy attenuation caused by the interaction and coupling effect between solid and fluid. Therefore, after this treatment, the medium degenerates from a two-phase viscoelastic TTI medium to a single-phase viscoelastic TTI medium; secondly, the polarization angle θ is set to 0, the splitting of fast P-wave and S-wave exists, and the overall energy of the two components does not change significantly, but the wave front morphology rotates, and the symmetry axis is tilted from 45° to 45°. It is changed to coincide with the central vertical line. Such a wave field conforms to the wave field characteristics of a single-phase VTI medium. Combined with the invariance of the medium energy, it can be judged that the medium has evolved from a single-phase viscoelastic TTI medium to a single-phase viscoelastic VTI medium. Again, based on the previous forward modeling, the anisotropy parameter ε = γ = δ = 0, and the seismic wave field has only one type of wave, namely, the longitudinal wave (because the forward simulation of this application example uses a longitudinal wave source, without separating the longitudinal and transverse wave fields, no transverse wave will be generated in a uniform isotropic medium), and the wave field energy is not much different. It can be seen that the medium degenerates into a single-phase viscoelastic medium. Finally, the longitudinal and transverse wave quality factors Q P and Q S When the value is infinite, the medium degenerates into a uniform isotropic perfectly elastic medium, the wave field morphology does not change, but the wave field energy is significantly enhanced, which conforms to the wave field characteristics of an isotropic perfectly elastic medium. Figure 2 and Figure 3 From the results, it can be seen that as long as the values ​​of relevant parameters are adjusted, the forward modeling results of various corresponding types of media can be obtained based on the finite difference forward simulation of the present invention, and finally the wave field results of isotropic perfectly elastic media can be degenerated, thereby verifying the accuracy and comprehensiveness of the wave equation.

[0205] (b) Optimization of two-way wave illumination and observation system for depression work area model

[0206] A two-phase viscoelastic TTI medium model is established based on a depression model of a certain work area. The model grid size is 288×151, the grid spacing is Δx=Δz=10m, and multiple shots are excited on the surface. The number of shots is 20, and the shot spacing is 140m. The source is uniformly a Ricker wavelet longitudinal wave source with a main frequency of 25Hz. The time sampling interval of the forward simulation is 0.5ms. The distribution of the parameter values ​​of the medium model is shown in Figure 4 , the first row from left to right is v P 、v S , P1, the second row from left to right is ρ 12 , 22 , 11 , the third row from left to right is δ, γ, ε. Figure 5A and Figure 5B As shown in the figure, the illumination results of the solid and liquid phases show abnormal values, among which the illumination energy of the liquid phase shows a negative number, which is obviously inconsistent with reality. Figure 5C The source illumination result calculated by the two-phase illumination operator is normal in shape. This result shows the correctness of the two-phase illumination operator, and the subsequent reverse illumination and optimized encryption results can further prove this point. Fig. 6A is the source lighting result, Figure 6B is the reverse illumination result of the model. Fig. 6A The illumination result of the source direction is intercepted, and the illumination result of the area below 400m is intercepted. The normalized isoline of the illumination energy is drawn. The shadow area is located according to the isoline, and the excitation area and excitation position of the reverse illumination are determined. The isoline is shown in Figure 7 . Figure 8 The normalized curve of the surface energy of reverse illumination is shown. From the figure, it can be determined that the energy maximum range is roughly around several wave peaks between 200 and 700 meters. The energy in the horizontal distance between 1500 and 2800 meters is very low. For energy optimization of the target area, it is possible to consider reducing the density of shot points in this range. According to the reverse illumination positioning encryption interval, the shot points are encrypted in the encrypted interval, and the number of shots in the low-energy interval indicated by the reverse illumination surface energy curve is correspondingly reduced. The encrypted source illumination results are shown in Fig. 9 . Fig.10 The comparison results of the normalized illumination energy of the shadow area of ​​the medium model before and after encryption are shown. From the reflection of this result, it can be seen that the illumination energy of the shadow area of ​​the observation system after the densification of the shot points has been significantly enhanced. This result conforms to the reciprocal principle. Therefore, the observation system is optimized, and the correctness of the two-phase illumination operator and the two-phase viscoelastic TTI medium wave equation are verified.

[0207] (c) Optimization of two-way wave illumination and observation system of Marmours model

[0208] The Marmours model of two-phase viscoelastic TTI properties is established. The model grid size is 298×170, the grid spacing is Δx=Δz=10m, and multiple shots are excited on the surface. The number of shots is 20, and the shot spacing is 150m. The source is uniformly a Ricker wavelet longitudinal wave source with a main frequency of 25Hz. The time sampling interval of the forward simulation is 0.5ms. The distribution of the parameter values ​​of the medium model is shown in Fig.11 , the first row from left to right is v P 、v S , P1, the second row from left to right is ρ 12 , 22 , 11 , the third row from left to right is δ, γ, ε. Fig. 12A It is the source illumination result, intercepting the area below 1200m in depth, and locating the shadow area according to the regional energy distribution. The energy value distribution is shown in Fig. 12B . Fig.13A is the reverse lighting result of the Marmours model, Fig. 13B The normalized curve of reverse illumination surface energy is obtained. The high energy interval determined from the figure is roughly around several wave peaks between 1250 and 2100 m, and the low energy interval is within the horizontal distance of 0 to 1100 m. The number of shot points in the high energy interval of reverse illumination surface energy is increased, and the number of shots in the low energy interval is reduced. The source illumination result after the shot point encryption is shown in Fig.14A and 14B . Fig.15 The comparison results of the normalized illumination energy of the shadow area of ​​the Marmours model before and after encryption are shown. The illumination energy of the shadow area has been significantly enhanced after encryption. This result also conforms to the reciprocity principle, which shows that the observation system has been optimized. At the same time, the correctness and applicability of the illumination operator and wave equation have been verified.

[0209] The present invention combines the anisotropy theory of TI media with the Carcione viscoelastic anisotropy theory to describe the constitutive relationship of solid media, and applies its mechanism to the Biot two-phase medium system, and successfully derives the wave equation of the two-phase viscoelastic TTI medium. The wave equation is numerically simulated using the high-order staggered grid finite difference method, and the results clearly reflect the various wave field information of the two-phase viscoelastic TTI medium. The numerical simulation results include fast longitudinal waves, shear waves and slow longitudinal waves, and its wave field characteristics have the characteristics of TTI media, viscoelastic media and two-phase media at the same time, which is consistent with the relevant theoretical predictions. By comparing the forward modeling of seismic waves of various types of media, the wave field characteristics of the two-phase viscoelastic TTI medium can be clearly demonstrated, verifying the validity of the derived wave equation. The seismic record curve generated by the forward modeling is used to qualitatively study the influence of relevant parameters on the wave field of the two-phase viscoelastic TTI medium, thereby verifying the relevant theory to a certain extent. In addition, the two-way wave illumination and observation system optimization of the two-phase viscoelastic TTI medium are performed using the newly established two-phase illumination operator. The results show that the calculation method of the two-phase illumination operator can well verify the principle of seismic wave reciprocity and the observation system optimization principle of the statistical illumination method. The illumination analysis process demonstrates the application of the two-phase viscoelastic TTI medium wave equation in complex models, and further verifies the correctness and applicability of the wave equation.

[0210] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium, characterized in that: The steps include: Step S1, introducing Biot two-phase medium theory to describe the two-phase medium; Step S2, using the TTI medium as a solid skeleton, deriving the constitutive equation of the TTI medium; Step S3, introducing Carcione's viscoelastic anisotropic medium model and combining it with the Boltzmann superposition principle to obtain the stress-strain relationship of the solid skeleton; Step S4, embedding the viscoelastic constitutive relation of the solid phase into Biot's two-phase medium theory, and deriving a wave equation applicable to the two-phase viscoelastic TTI medium; Step S5, obtaining the first-order time derivative of the stress and strain of the wave equation of the two-phase viscoelastic TTI medium to obtain the first-order velocity stress equation of the two-phase viscoelastic TTI medium, and using a high-order staggered grid finite difference method to perform seismic wave forward simulation on the derived first-order velocity stress equation of the two-phase viscoelastic TTI medium; Step S6, establishing an illumination operator, defining source direction and reverse illumination operators, calculating source direction illumination through forward simulation, and analyzing the illumination effects on the target area under different configurations; Step S7: Based on the reciprocity principle, analyze the reverse illumination results, identify the surface positions that contribute more to the target area, and optimize the excitation positions and shot point density to achieve optimal coverage; In step S2, the constitutive equation of the TTI medium is: T TTI =C TTI S TTI ; Where, T TTI is the TTI medium stress vector, S TTI is the TTI medium strain vector; In step S3, according to Carcione viscoelastic anisotropy theory and Boltzmann superposition principle, the stress-strain relationship of the solid skeleton is: In the formula, T=(σ xx ,σ yy ,σ zz ,σ yz ,σ zx ,σ xy ) is the stress vector, is the first-order time derivative of the strain vector, ψ represents the relaxation matrix; In step S4, the wave equation applicable to the two-phase viscoelastic TTI medium is: In the formula, U=(u x ,u y ,u z ,U x ,U y ,U z ) is the displacement vector; represents the second-order derivative of U with respect to time; F = (f x ,f y ,f z ,F x ,F y ,F z ) is the body force vector, where f x 、f y 、f z represents the solid phase fraction, F x 、F y 、F z represents the liquid phase fraction, Compared to A 7×7 matrix with a row vector and a column vector with 0 elements is added respectively; ρ is the density matrix; G is the partial differential operator; Β is the dissipation coefficient matrix; J is the sum and difference operator.

2. The seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium according to claim 1, characterized in that: In step S5, the first-order time derivative equation of the memory variable is: In the formula, the relaxation time τ is calculated by the quality factor, and the calculation formula is: Where v = 1, 2, Q1 and Q2 are the longitudinal wave quality factor and the shear wave quality factor, and ω is the reciprocal of the main frequency of the source wavelet.

3. The seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium according to claim 2, characterized in that: Taking into account the viscoelasticity and TI anisotropy of the medium, the necessary stability condition for the forward simulation of the first-order velocity stress equation of the two-phase viscoelastic TTI medium using the high-order staggered grid finite difference method is: In the formula, 4. The seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium according to claim 3, characterized in that: In step S6, the illumination operator expression is: Source lighting operator: Where, X s (x s ,0) are the excitation positions of the observation system, X(x,z) are the coordinates of any point D in the two-phase viscoelastic medium model area; The reverse illumination operator is: In the formula, X g (x g ,0) are the coordinates of each position in the observation range, X o (x o ,z o ) is the target point of the source lighting in the shadow area, and D′ is the shadow area.

5. The seismic wave forward modeling and two-way wave illumination method for a two-phase viscoelastic TTI medium according to claim 4, characterized in that: In step S7, the illumination function of the optimized and encrypted observation system for the two-phase medium is: Where, X s ′(x s ′,0) is the encrypted surface observation system.

Images