CN116794716B - A dispersive AVO simulation method for mesoscopic fracture rock physics model - Google Patents

Abstract

The invention discloses a dispersion AVO simulation method of a mesoscopic fracture rock physical model, which belongs to the technical field of oil and gas exploration. It includes a frequency-varying stiffness matrix acquisition stage of the mesoscopic fracture rock and a model dispersion AVO simulation stage; the mesoscopic In the acquisition stage of the frequency-varying stiffness matrix of fractured rock, the corresponding additional flexibility matrix can be constructed for different fracture arrangements. Combined with the frequency-dependent fracture normal flexibility, the corresponding dynamic flexibility matrix of dry rock can be obtained. In the fluid characteristics After adding, the frequency-varying stiffness matrix of the rock is obtained; in the dispersion AVO simulation stage of the model, the frequency-dependent anisotropy coefficient is obtained based on the frequency-varying stiffness matrix of the rock, and combined with the frequency-varying interface reflection coefficient and the source wavelet, the calculation is Multi-interface seismic dispersion AVO simulation recording. The invention solves the problem of difficulty in rock model construction and dispersion AVO analysis of randomly distributed mesoscopic fracture dynamic models.

Description

A Dispersion AVO Simulation Method for Mesoscopic Fracture Rock Physics Model

Technical Field

The invention belongs to the technical field of oil and gas exploration, and in particular relates to a dispersion AVO simulation method of a mesoscopic fracture rock physics model.

Background Art

When waves pass through fluid-containing rocks underground, velocity dispersion and energy attenuation (abbreviated as dispersion) usually occur. Currently, the widely accepted theory about the mechanism of this dispersion effect is wave-induced fluid flow (WIFF), and the dispersion attenuation mechanism occurring in the seismic exploration frequency band is mainly mesoscopic flow.

The mesoscale mainly refers to the scale that is much smaller than the wavelength of artificial seismic waves and much larger than the size of rock pores, usually at the centimeter level. There are many reasons for the mesoscale heterogeneity of rocks, and cracks are one of the main factors. In addition, cracks may have a certain directionality, which will cause the dispersion and attenuation of waves to be anisotropic. At present, the description of fractured rocks, especially at the mesoscale, mostly only considers directional arrangements, and lacks modeling schemes for complex arrangements or even randomly arranged mesoscopic fracture models. Therefore, a method that can construct a dynamic model of randomly distributed and arranged mesoscopic fractures and perform dispersion AVO simulation on the constructed model is urgently needed.

Summary of the invention

In view of the above-mentioned deficiencies in the prior art, the present invention provides a dispersion AVO simulation method for a mesoscopic fracture rock physics model, which solves the problem that it is difficult to construct a rock model and perform dispersion AVO analysis on a randomly distributed mesoscopic fracture dynamic model.

In order to achieve the above-mentioned object of the invention, the technical solution adopted by the present invention is:

The present invention provides a dispersion AVO simulation method for a mesoscopic fracture rock physics model, comprising a frequency-variable stiffness matrix acquisition stage of the mesoscopic fracture rock, and a model dispersion AVO simulation stage;

The frequency-dependent stiffness matrix acquisition stage of the mesoscopic fractured rock includes the following steps:

A1. Obtain the rock crack density, crack number, crack radius and shear modulus of the background rock, and calculate the normal flexibility and tangential flexibility of the crack;

A2. Based on the normal flexibility of the fracture and the P-wave modulus of the saturated fractured rock and background rock, the frequency-dependent normal flexibility of the fracture is calculated.

A3. According to the arrangement of the cracks, the dynamic flexibility matrix of the dry rock is calculated based on the tangential flexibility of the cracks, the frequency-dependent normal flexibility of the cracks and the dry rock background flexibility matrix;

A4. Based on the dry rock dynamic flexibility matrix and the anisotropic Gassmann equation, the complete rock frequency-dependent flexibility matrix is calculated;

A5. Invert the complete rock frequency-variant flexibility matrix to obtain the rock frequency-variant stiffness matrix;

The model dispersion AVO simulation stage includes the following steps:

B1. Based on the frequency-dependent stiffness matrix of rock, the frequency-dependent stiffness coefficient is obtained;

B2. Based on the real part of the frequency-dependent stiffness coefficient, the frequency-dependent anisotropy coefficient is calculated;

B3. Based on the anisotropy coefficient and reflection coefficient equation, the frequency-dependent interface reflection coefficient is obtained;

B4. Based on the frequency-varying interface reflection coefficient and source wavelet, the multi-interface seismic dispersion AVO simulation record is calculated.

The beneficial effects of the present invention are as follows: the present invention provides a dispersion AVO simulation method for a mesoscopic fracture rock physics model, which constructs a frequency-dependent fracture flexibility parameter through the frequency-dependent stiffness matrix acquisition stage of the mesoscopic fracture rock, and a complex dispersion model can be constructed based on the frequency-dependent fracture flexibility parameter; the model constructed by this scheme is not limited to describing HTI media, and can realize the construction of a dynamic flexibility matrix of dry rock under random distribution and any fracture arrangement, and the calculation method for constructing the model is more convenient and faster than the existing construction method; by introducing dispersion characteristics, the construction of a dynamic model containing fluid can be realized, which greatly enriches the diversity of mesoscopic fracture rock physics models; by analyzing the AOV frequency variation and frequency-related synthetic seismic records of the dispersion rock physics model in the model dispersion AVO simulation stage, more parameters or phenomena can be provided to describe the fluid and fracture conditions of underground rocks.

Furthermore, the calculation expressions of the normal flexibility and the tangential flexibility in A1 are as follows:

in, and They represent the normal flexibility and tangential flexibility of the dry crack, λ represents the Lame coefficient, μ represents the shear modulus of the background rock, e represents the crack density, N represents the number of cracks, a represents the crack radius, and V represents the rock volume.

The beneficial effect of adopting the above further scheme is: providing a calculation method for normal flexibility and tangential flexibility, obtaining the normal flexibility and tangential flexibility of cracks in dry rock based on relevant characteristics such as crack density in rock and shear modulus of background rock, providing a basis for obtaining frequency-dependent crack normal flexibility.

Furthermore, the A2 comprises the following steps:

A21. Based on the P-wave modulus of saturated fractured rock and background rock, the time scale parameter and shape parameter are calculated;

The calculation expressions of the time scale parameter and shape parameter are as follows:

Among them, τ and denote the time scale parameter and shape parameter respectively, _Cb denotes the P-wave modulus of saturated fractured rock and background rock under high-frequency limit conditions, _C0 denotes the P-wave modulus of saturated fractured rock and background rock under low-frequency limit conditions, G denotes the low-frequency scale parameter, and T denotes the high-frequency scale parameter;

A22, based on the time scale parameter and the shape scale parameter, the flexibility frequency relationship is calculated;

The calculation expression of the flexibility frequency relationship is as follows:

Where f _fra (ω) represents the flexibility frequency relationship, i represents the imaginary part of the complex number, and ω represents the frequency;

A23, based on the flexibility frequency relationship and the normal flexibility of the crack, the frequency-dependent normal flexibility of the crack is calculated;

The calculation expression of the frequency-dependent crack normal flexibility is as follows:

where Z _N (ω) represents the frequency-dependent crack normal compliance.

The beneficial effects of adopting the above further scheme are: based on the time scale parameters, shape parameters, P-wave moduli of saturated fractured rock and background rock under low-frequency and high-frequency extreme conditions respectively, and the characteristics of fluid, fracture and rock mineral particles, the flexibility frequency relationship is obtained, and the normal flexibility is assigned frequency characteristics, providing a basis for constructing a rock dynamic flexibility matrix for different arrangements of fractures in the vertical direction.

Furthermore, the A3 comprises the following steps:

A31. Obtain the arrangement of fractures in the mesoscopic fracture rock physics model;

A32, determine whether the arrangement of the cracks is vertical random arrangement, if so, proceed to A33, otherwise proceed to A34;

A33, based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, calculate the additional flexibility matrix of the vertical random crack, and use the additional flexibility matrix of the vertical random crack as the frequency-dependent additional flexibility matrix, and enter step A37;

The calculation expression of the additional flexibility matrix of the vertical random crack is as follows:

Where ΔS(ω) _s represents the additional flexibility matrix of vertical random cracks, Z _N1 (ω) represents the frequency-dependent normal flexibility of cracks when the cracks are arranged vertically randomly, It represents the tangential flexibility of dry cracks when cracks are arranged vertically and randomly;

A34, judging whether the vertical arrangement of the cracks is vertical directional arrangement, if so, proceeding to A35, otherwise proceeding to A36;

A35, based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, calculate the additional flexibility matrix of the vertical oriented crack, and use the additional flexibility matrix of the vertical oriented crack as the frequency-dependent additional flexibility matrix, and enter step A37;

The calculation expression of the additional flexibility matrix of the vertical oriented crack is as follows:

where ΔS(ω) _d represents the dry rock additional flexibility matrix for fracture orientation, Z _N (ω) ₂ represents the frequency-dependent fracture normal flexibility when the fractures are arranged vertically, It represents the tangential flexibility of dry cracks when the cracks are arranged vertically;

A36, judging whether the vertical arrangement of the cracks is vertically oriented and randomly arranged, if so, proceeding to A37, otherwise proceeding to A38;

A37. Based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, a directional and random additional flexibility matrix of the vertical crack is calculated, and the directional and random additional flexibility matrix of the vertical crack is used as the frequency-dependent additional flexibility matrix;

The calculation expression of the additional flexibility matrix of the vertical crack orientation and randomness is as follows:

Where, ΔS(ω) _sd represents the additional flexibility matrix of the vertical crack which is oriented and random;

A38. Based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, a completely random additional flexibility matrix of the crack is calculated, and the completely random additional flexibility matrix of the crack is used as the frequency-dependent additional flexibility matrix;

The calculation expression of the additional flexibility matrix of the completely random crack is as follows:

Where ΔS(ω) ₃ represents the additional flexibility matrix when the cracks are completely random, Z _N (ω) ₃ represents the frequency-dependent normal flexibility of the cracks when the cracks are completely randomly arranged, It represents the tangential flexibility of dry cracks when the cracks are arranged completely randomly;

A39, based on the dry rock background flexibility matrix and the frequency-dependent additional flexibility matrix, the dry rock dynamic flexibility matrix is calculated;

The calculation expression of the dry rock dynamic flexibility matrix is as follows:

Where S ⁰ (ω) represents the dynamic flexibility matrix of dry rock, represents the dry rock background flexibility matrix, and ΔS(ω) represents the frequency-dependent flexibility matrix.

The beneficial effects of adopting the above further scheme are: a calculation method of the additional flexibility matrix is provided for the cases of random arrangement direction and random distribution in the vertical direction of the mesoscopic fracture rock physics model, and the additional flexibility matrix is combined with the dry rock background flexibility matrix to obtain the dry rock dynamic flexibility matrix, which provides a basis for obtaining a complete rock frequency-variable flexibility matrix.

Furthermore, the calculation expression of the complete rock frequency-dependent flexibility matrix in A4 is as follows:

in, and denote the flexibility tensors of saturated rock, dry rock and mineral particles, respectively, represents the sum of the flexibility of the i′th row in the matrix corresponding to the dry rock flexibility tensor, represents the sum of the flexibility of the i′th row in the matrix corresponding to the flexibility tensor of the mineral particles, represents the sum of the flexibility of the j′th column in the matrix corresponding to the dry rock flexibility tensor, represents the sum of the flexibility of the j′th column in the matrix corresponding to the flexibility tensor of the mineral particles, β _dry represents the generalized compressibility coefficient of dry rock with cracks, β _g represents the compressibility coefficient of mineral particles, β _f represents the compressibility coefficient of fluid, Represents the total porosity, where i′=1,2,3,4,5,6 and j′=1,2,3,4,5,6.

The beneficial effect of adopting the above further scheme is: the dynamic flexibility matrix of dry rock is substituted into the anisotropic Gassmann equation, and the residual fluid effect is added. On the basis of dry rock, the complete rock frequency-variant flexibility matrix can be obtained by combining the fluid characteristics, which provides a basis for obtaining the frequency-variant stiffness matrix of rock and performing dispersion AVO simulation on the constructed rock model.

Furthermore, the calculation expression of the anisotropy coefficient in B2 is as follows:

Wherein, α(ω) represents the first frequency-dependent anisotropy coefficient, c ₃₃ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 3rd row and 3rd column in the frequency-dependent stiffness matrix, ρ represents the medium density, β(ω) represents the second frequency-dependent anisotropy coefficient, c ₄₄ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 4th row and 4th column in the frequency-dependent stiffness matrix, ε(ω) represents the third frequency-dependent anisotropy coefficient, c ₁₁ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 1st column in the frequency-dependent stiffness matrix, γ(ω) represents the fourth frequency-dependent anisotropy coefficient, c ₆₆ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 6th row and 6th column in the frequency-dependent stiffness matrix, δ(ω) represents the fifth frequency-dependent anisotropy coefficient, c ₁₃ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 3rd column in the frequency-dependent stiffness matrix, c ₅₅ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 3rd column in the frequency-dependent stiffness matrix, (ω) represents the real part of the frequency-dependent stiffness coefficient at the 5th row and 5th column in the frequency-dependent stiffness matrix.

The beneficial effect of adopting the above further scheme is: by calculating the coefficients in the frequency-dependent stiffness matrix of the rock corresponding to the mesoscopic fracture rock physics model with respect to the medium density, the frequency-dependent anisotropic coefficient is obtained, which provides a basis for calculating the interface reflection coefficient.

Furthermore, the calculation expression of the frequency-dependent interface reflection coefficient in B3 is as follows:

Δγ(ω)=γ ₂ (ω)-γ ₁ (ω)

Δα(ω)＝α ₂ (ω)-α ₁ (ω)Δε(ω)＝ε ₂ (ω)-ε ₁ (ω)Δδ(ω)＝δ ₂ (ω)-δ ₁ (ω)

ΔZ(ω)=Z ₂ (ω)-Z ₁ (ω)Z(ω)=ρα(ω)

Δμ(ω)=μ ₂ (ω)-μ ₁ (ω)μ(ω)=ρβ(ω) ²

Where R _pp (ω,ξ,θ) represents the frequency-dependent interface reflection coefficient, ΔZ(ω) represents the frequency-dependent vertical P-wave impedance difference between the lower and upper layers, represents the average value of the vertical P-wave impedance of the lower and upper layers related to frequency, Δα(ω) represents the difference of the first anisotropy coefficient of the lower and upper layers dependent on frequency, represents the average value of the first anisotropy coefficient of the lower and upper layers in terms of frequency dependence, represents the average value of the second anisotropy coefficient of the lower and upper layers in terms of frequency, Δμ(ω) represents the difference in the vertical shear wave modulus of the lower and upper layers in terms of frequency, represents the average value of the vertical shear wave modulus of the lower and upper media related to frequency, Δδ(ω) represents the difference of the fifth anisotropy coefficient of the lower and upper media depending on frequency, θ represents the P-wave incident angle, Δε(ω) represents the difference of the third anisotropy coefficient of the lower and upper media depending on frequency, ξ represents the P-wave incident azimuth, α ₁ (ω) and α ₂ (ω) represent the first anisotropy coefficient of the upper and lower media depending on frequency, β ₁ (ω) and β ₂ (ω) represent the second anisotropy coefficient of the upper and lower media depending on frequency, ε ₁ (ω) and ε ₂ (ω) represent the third anisotropy coefficient of the upper and lower media depending on frequency, δ ₁ (ω) and δ ₂ (ω) represent the fifth anisotropy coefficient of the upper and lower media depending on frequency, γ ₁ (ω) and γ ₂ (ω) (ω) represents the fourth anisotropy coefficient of the upper medium and the lower medium in dependence of frequency, Δγ(ω) represents the difference of the fourth anisotropy coefficient of the lower medium and the upper medium in dependence of frequency, Z ₁ (ω) and Z ₂ (ω) represent the vertical P-wave impedance of the lower medium and the upper medium in dependence of frequency, Z(ω) represents the vertical P-wave impedance in dependence of frequency, μ ₁ and μ ₂ represent the vertical shear wave modulus of the lower medium and the upper medium in dependence of frequency, respectively.

The beneficial effect of adopting the above further scheme is: substituting the frequency-dependent anisotropy coefficient into the HIT medium P-wave reflection coefficient formula to obtain the frequency-dependent interface reflection coefficient, which provides a basis for seismic recording combined with the AVO simulation results of the source wavelet.

Further, the B4 comprises the following steps:

B41. Based on the frequency-dependent interface reflection coefficient, the frequency-dependent reflection coefficient of the interface is obtained;

B42, based on the interface-frequency-dependent reflection coefficient and source wavelet, the multi-interface seismic dispersion AVO simulation record is calculated;

The calculation expression of the multi-interface seismic dispersion AVO simulation record is as follows:

where x _t (θ) is the earthquake record at time depth t, represents the unit reflection coefficient of the nth interface, k represents the total number of interface layers, Δt represents the sampling interval, W(ω) represents the source wavelet, R ⁿ (θ,ω) represents the frequency-dependent reflection coefficient of the nth interface, * represents the convolution operation, and FFT ^-1 represents the inverse Fourier transform.

The beneficial effect of adopting the above further scheme is: in the frequency domain, the frequency-variant reflection coefficient is multiplied by the source wavelet, and an inverse Fourier transform is performed to obtain a synthetic seismic record. For multi-layer interfaces, the target reservoir is divided according to different time depths to obtain interfaces of different time depths, so that synthetic seismic records of multiple interfaces can be obtained, thereby realizing multi-interface seismic dispersion AVO simulation.

Compared with the prior art, the present invention has the following advantages:

The present invention provides a dispersion AVO simulation method for a mesoscopic fracture rock physics model, which can realize dynamic rock simulation of complex fracture conditions. When rocks are affected by different geological conditions at different times, multiple groups of fracture conditions are generated. When the rock fracture arrangement has vertical directional arrangement, vertical random arrangement and completely random arrangement at the same time, a calculation method for constructing a frequency-dependent additional flexibility matrix can still be provided. The construction method provided by the present invention is based on frequency characteristics, which is simple and efficient but not accurate.

Other advantages of the present invention will be analyzed in more detail in subsequent embodiments.

BRIEF DESCRIPTION OF THE DRAWINGS

In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments are briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without creative work.

FIG1 is a flow chart showing the steps of a dispersion AVO simulation method for a mesoscopic fracture rock physics model according to an embodiment of the present invention.

FIG. 2 is a schematic diagram of a horizontal layer in which cracks are arranged in a vertically oriented manner in an embodiment of the present invention.

FIG. 3( a ) is a graph showing how the reflection coefficient of an interface numbered H varies with the incident angle in an embodiment of the present invention.

FIG. 3( b ) is a graph showing how the reflection coefficient of the interface numbered H varies with frequency in an embodiment of the present invention.

FIG3( c ) is a schematic diagram of an AVO simulation record of an interface seismic dispersion numbered H in an embodiment of the present invention.

FIG. 3( d ) is a waveform diagram of superimposed seismic records of the interface numbered H in an embodiment of the present invention that is frequency-dependent and frequency-independent.

FIG. 4( a ) is a graph showing how the reflection coefficient of an interface numbered M varies with the incident angle in an embodiment of the present invention.

FIG. 4( b ) is a graph showing how the reflection coefficient of an interface numbered M varies with frequency in an embodiment of the present invention.

FIG4( c ) is a schematic diagram of an AVO simulation record of seismic dispersion at an interface numbered M in an embodiment of the present invention.

FIG. 4( d ) is a waveform diagram of superimposed seismic records of the interface numbered M in an embodiment of the present invention that is frequency-dependent and frequency-independent.

FIG. 5( a ) is a graph showing how the reflection coefficient of an interface numbered S varies with the incident angle in an embodiment of the present invention.

FIG5( b ) is a schematic diagram of an AVO simulation record of seismic dispersion at an interface numbered S in an embodiment of the present invention.

FIG. 5( c ) is a waveform diagram of superimposed seismic records of the interface numbered S in an embodiment of the present invention that is frequency-dependent and frequency-independent.

DETAILED DESCRIPTION

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 only a part of the embodiments of the present invention, rather than all the embodiments. The components of the embodiments of the present invention generally described and shown in the drawings here can be arranged and designed in various different configurations. 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 claimed invention, but merely represents the selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making creative work belong to the scope of protection of the present invention.

Example 1

As shown in FIG1 , in one embodiment of the present invention, the present invention provides a dispersion AVO simulation method for a mesoscopic fracture rock physics model, including a frequency-variant stiffness matrix acquisition stage of the mesoscopic fracture rock, and a model dispersion AVO simulation stage;

The frequency-dependent stiffness matrix acquisition stage of the mesoscopic fractured rock includes the following steps:

A1. Obtain the rock crack density, crack number, crack radius and shear modulus of the background rock, and calculate the normal flexibility and tangential flexibility of the crack;

The calculation expressions of the normal flexibility and tangential flexibility in A1 are as follows:

in, and They represent the normal flexibility and tangential flexibility of the dry crack, λ represents the Lame coefficient, μ represents the shear modulus of the background rock, e represents the crack density, N represents the number of cracks, a represents the crack radius, and V represents the rock volume.

A2. Based on the normal flexibility of the fracture and the P-wave modulus of the saturated fractured rock and background rock, the frequency-dependent normal flexibility of the fracture is calculated.

A2 comprises the following steps:

A21. Based on the P-wave modulus of saturated fractured rock and background rock, the time scale parameter and shape parameter are calculated;

The calculation expressions of the time scale parameter and shape parameter are as follows:

Among them, τ and denote the time scale parameter and shape parameter respectively, _Cb denotes the P-wave modulus of saturated fractured rock and background rock under high-frequency limit conditions, _C0 denotes the P-wave modulus of saturated fractured rock and background rock under low-frequency limit conditions, G denotes the low-frequency scale parameter, and T denotes the high-frequency scale parameter;

At present, the most common crack is a coin-shaped crack. This scheme provides a calculation method for the high-frequency scale parameters and low-frequency scale parameters of a coin-shaped crack. The calculation expressions for the high-frequency scale parameters and low-frequency scale parameters of a coin-shaped crack are as follows:

Wherein, α _b represents the Biot-Willis coefficient of the background rock, M _b represents the fluid-solid coupling modulus of the background rock, g _b represents the ratio of the shear modulus to the P-wave modulus in the background rock, η represents the fluid viscosity coefficient, μ _b represents the shear modulus of the background rock, _LB represents the P-wave modulus of the background rock, κ _b represents the permeability of the background medium, K _b , K _g , and K _f represent the bulk moduli of the dry background rock, mineral particles, and fluid, respectively. Represents the background porosity;

A22, based on the time scale parameter and the shape scale parameter, the flexibility frequency relationship is calculated;

The calculation expression of the flexibility frequency relationship is as follows:

Where f _fra (ω) represents the flexibility frequency relationship, i represents the imaginary part of the complex number, and ω represents the frequency;

A23, based on the flexibility frequency relationship and the normal flexibility of the crack, the frequency-dependent normal flexibility of the crack is calculated;

The calculation expression of the frequency-dependent crack normal flexibility is as follows:

where Z _N (ω) represents the frequency-dependent crack normal compliance.

A3. According to the arrangement of the cracks, the dynamic flexibility matrix of the dry rock is calculated based on the tangential flexibility of the cracks, the frequency-dependent normal flexibility of the cracks and the dry rock background flexibility matrix;

A3 comprises the following steps:

A31. Obtain the arrangement of fractures in the mesoscopic fracture rock physics model;

A32, determine whether the arrangement of the cracks is vertical random arrangement, if so, proceed to A33, otherwise proceed to A34;

A33, based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, calculate the additional flexibility matrix of the vertical random crack, and use the additional flexibility matrix of the vertical random crack as the frequency-dependent additional flexibility matrix, and enter step A37;

The calculation expression of the additional flexibility matrix of the vertical random crack is as follows:

Where ΔS(ω) _s represents the additional flexibility matrix of vertical random cracks, Z _N1 (ω) represents the frequency-dependent normal flexibility of cracks when the cracks are arranged vertically randomly, It represents the tangential flexibility of dry cracks when cracks are arranged vertically and randomly;

A34, judging whether the vertical arrangement of the cracks is vertical directional arrangement, if so, proceeding to A35, otherwise proceeding to A36;

A35, based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, calculate the additional flexibility matrix of the vertical oriented crack, and use the additional flexibility matrix of the vertical oriented crack as the frequency-dependent additional flexibility matrix, and enter step A37;

The calculation expression of the additional flexibility matrix of the vertical oriented crack is as follows:

where ΔS(ω) _d represents the dry rock additional flexibility matrix for fracture orientation, Z _N (ω) ₂ represents the frequency-dependent fracture normal flexibility when the fractures are arranged vertically, It represents the tangential flexibility of dry cracks when the cracks are arranged vertically;

A36, judging whether the vertical arrangement of the cracks is vertically oriented and randomly arranged, if so, proceeding to A37, otherwise proceeding to A38;

A37. Based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, a directional and random additional flexibility matrix of the vertical crack is calculated, and the directional and random additional flexibility matrix of the vertical crack is used as the frequency-dependent additional flexibility matrix;

The calculation expression of the additional flexibility matrix of the vertical crack orientation and randomness is as follows:

Where, ΔS(ω) _sd represents the additional flexibility matrix of the vertical crack which is oriented and random;

A38. Based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, a completely random additional flexibility matrix of the crack is calculated, and the completely random additional flexibility matrix of the crack is used as the frequency-dependent additional flexibility matrix;

The calculation expression of the additional flexibility matrix of the completely random crack is as follows:

Where ΔS(ω) ₃ represents the additional flexibility matrix when the cracks are completely random, Z _N (ω) ₃ represents the frequency-dependent normal flexibility of the cracks when the cracks are completely randomly arranged, It represents the tangential flexibility of dry cracks when the cracks are arranged completely randomly;

A39, based on the dry rock background flexibility matrix and the frequency-dependent additional flexibility matrix, the dry rock dynamic flexibility matrix is calculated;

The calculation expression of the dry rock dynamic flexibility matrix is as follows:

Where S ⁰ (ω) represents the dynamic flexibility matrix of dry rock, represents the dry rock background flexibility matrix, and ΔS(ω) represents the frequency-dependent flexibility matrix.

As a preferred solution, in this embodiment, the background flexibility matrix of the rock is obtained from the background stiffness matrix; the parameters in the background stiffness matrix can be obtained by calculating the mineral particle parameters and background porosity of the rock, or by measuring the high-pressure experiment.

A4. Based on the dry rock dynamic flexibility matrix and the anisotropic Gassmann equation, the complete rock frequency-dependent flexibility matrix is calculated;

The calculation expression of the complete rock frequency-dependent flexibility matrix in A4 is as follows:

in, and denote the flexibility tensors of saturated rock, dry rock and mineral particles, respectively, represents the sum of the flexibility of the i′th row in the matrix corresponding to the dry rock flexibility tensor, represents the sum of the flexibility of the i′th row in the matrix corresponding to the flexibility tensor of the mineral particles, represents the sum of the flexibility of the j′th column in the matrix corresponding to the dry rock flexibility tensor, represents the sum of the flexibility of the j′th column in the matrix corresponding to the flexibility tensor of the mineral particles, β _dry represents the generalized compressibility coefficient of dry rock with cracks, β _g represents the compressibility coefficient of mineral particles, β _f represents the compressibility coefficient of fluid, Represents the total porosity, where i′=1,2,3,4,5,6 and j′=1,2,3,4,5,6.

A5. Invert the complete rock frequency-variant flexibility matrix to obtain the rock frequency-variant stiffness matrix;

The model dispersion AVO simulation stage includes the following steps:

B1. Based on the frequency-dependent stiffness matrix of rock, the frequency-dependent stiffness coefficient is obtained;

B2. Based on the real part of the frequency-dependent stiffness coefficient, the frequency-dependent anisotropy coefficient is calculated;

The calculation expression of the anisotropy coefficient in B2 is as follows:

Wherein, α(ω) represents the first frequency-dependent anisotropy coefficient, c ₃₃ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 3rd row and 3rd column in the frequency-dependent stiffness matrix, ρ represents the medium density, β(ω) represents the second frequency-dependent anisotropy coefficient, c ₄₄ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 4th row and 4th column in the frequency-dependent stiffness matrix, ε(ω) represents the third frequency-dependent anisotropy coefficient, c ₁₁ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 1st column in the frequency-dependent stiffness matrix, γ(ω) represents the fourth frequency-dependent anisotropy coefficient, c ₆₆ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 6th row and 6th column in the frequency-dependent stiffness matrix, δ(ω) represents the fifth frequency-dependent anisotropy coefficient, c ₁₃ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 3rd column in the frequency-dependent stiffness matrix, c ₅₅ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 3rd column in the frequency-dependent stiffness matrix, (ω) represents the real part of the frequency-dependent stiffness coefficient at the 5th row and 5th column in the frequency-dependent stiffness matrix.

B3. Based on the anisotropy coefficient and reflection coefficient equation, the frequency-dependent interface reflection coefficient is obtained;

The calculation expression of the frequency-dependent interface reflection coefficient in B3 is as follows:

Δγ(ω)=γ ₂ (ω)-γ ₁ (ω)

Δα(ω)＝α ₂ (ω)-α ₁ (ω)Δε(ω)＝ε ₂ (ω)-ε ₁ (ω)Δδ(ω)＝δ ₂ (ω)-δ ₁ (ω)

ΔZ(ω)=Z ₂ (ω)-Z ₁ (ω)Z(ω)=ρα(ω)

Δμ(ω)=μ ₂ (ω)-μ ₁ (ω)μ(ω)=ρβ(ω) ²

Where R _pp (ω,ξ,θ) represents the frequency-dependent interface reflection coefficient, ΔZ(ω) represents the frequency-dependent vertical P-wave impedance difference between the lower and upper layers, represents the average value of the vertical P-wave impedance of the lower and upper layers related to frequency, Δα(ω) represents the difference of the first anisotropy coefficient of the lower and upper layers dependent on frequency, represents the average value of the first anisotropy coefficient of the lower and upper layers in terms of frequency dependence, represents the average value of the second anisotropy coefficient of the lower and upper layers in terms of frequency, Δμ(ω) represents the difference in the vertical shear wave modulus of the lower and upper layers in terms of frequency, represents the average value of the vertical shear wave modulus of the lower and upper media related to frequency, Δδ(ω) represents the difference of the fifth anisotropy coefficient of the lower and upper media depending on frequency, θ represents the P-wave incident angle, Δε(ω) represents the difference of the third anisotropy coefficient of the lower and upper media depending on frequency, ξ represents the P-wave incident azimuth, α ₁ (ω) and α ₂ (ω) represent the first anisotropy coefficient of the upper and lower media depending on frequency, β ₁ (ω) and β ₂ (ω) represent the second anisotropy coefficient of the upper and lower media depending on frequency, ε ₁ (ω) and ε ₂ (ω) represent the third anisotropy coefficient of the upper and lower media depending on frequency, δ ₁ (ω) and δ ₂ (ω) represent the fifth anisotropy coefficient of the upper and lower media depending on frequency, γ ₁ (ω) and γ ₂ (ω) (ω) represents the fourth anisotropy coefficient of the upper medium and the lower medium in dependence of frequency, Δγ(ω) represents the difference of the fourth anisotropy coefficient of the lower medium and the upper medium in dependence of frequency, Z ₁ (ω) and Z ₂ (ω) represent the vertical P-wave impedance of the lower medium and the upper medium in dependence of frequency, Z(ω) represents the vertical P-wave impedance in dependence of frequency, μ ₁ and μ ₂ represent the vertical shear wave modulus of the lower medium and the upper medium in dependence of frequency, respectively.

B4. Calculate the multi-interface seismic dispersion AVO simulation record based on the frequency-varying interface reflection coefficient and source wavelet;

The B4 comprises the following steps:

B41. Based on the frequency-dependent interface reflection coefficient, the frequency-dependent reflection coefficient of the interface is obtained;

B42, based on the interface-frequency-dependent reflection coefficient and source wavelet, the multi-interface seismic dispersion AVO simulation record is calculated;

The calculation expression of the multi-interface seismic dispersion AVO simulation record is as follows:

where x _t (θ) is the earthquake record at time depth t, represents the unit reflection coefficient of the nth interface, k represents the total number of interface layers, Δt represents the sampling interval, W(ω) represents the source wavelet, R ⁿ (θ,ω) represents the frequency-dependent reflection coefficient of the nth interface, * represents the convolution operation, and FFT ^-1 represents the inverse Fourier transform.

The present invention adds a relaxation function to the crack parameters to construct a frequency-dependent crack flexibility parameter. This improvement facilitates the construction of complex dispersion models. Based on this method, the construction of a randomly distributed mesoscopic crack dynamic model can be realized, and the corresponding formula form is simpler than the previous model.

The rock physics model of the present invention is built and the simulation of related AVO seismic records is realized, which can provide a rock physics basis for the forward modeling and interpretation of seismic and logging data.

Based on the relevant dispersive rock physics model, analyzing its AVO frequency variation and frequency-related synthetic seismic records can provide more parameters or phenomena to describe the fluid and fracture conditions of underground rocks, so it is of great significance.

Example 2

As shown in FIG. 2 , in a practical example of the present invention, a horizontal layer is designed, the upper layer is dense mudstone, which is isotropic; the lower layer is sandstone, which contains fluid and fractures, wherein the fractures are arranged vertically (HTI type arrangement), and the mudstone parameters and the background rock parameters of the sandstone are shown in Table 1:

Table 1

In this embodiment, the background rock parameters are constant;

The fracture fluid parameters related to the horizontal layer are shown in Table 2:

Table 2

In this embodiment, according to the dispersion AVO simulation method of a mesoscopic fracture rock physics model provided in Example 1, a 25 Hz Ricker wavelet and the characteristic frequencies in Table 2 were used to perform dispersion AVO simulation on the designed horizontal layer, and multi-interface seismic dispersion AVO simulation records of a series of models in Table 2 were obtained.

When the sandstone parameters are model H, the characteristic frequency is about 158Hz, which is greater than the main frequency of the Ricker wavelet, 25Hz. As shown in Figure 3(a), there is no obvious reflection coefficient dispersion from 1Hz to 60Hz; as shown in Figure 3(b), the interface reflection coefficient is negative at vertical incidence, and as the incident angle increases, the reflection coefficient becomes positive and continues to increase; as shown in Figure 3(c), 0.3s is the frequency-dependent dynamic model simulation result, 0.7s is the frequency-independent low-frequency limit result simulated by the Gassmann equation, and the rightmost one is the result after the angle gather is superimposed; as shown in Figure 3(d), the dynamic model and Gassmann low-frequency limit results are almost consistent in the waveform diagram of the synthetic seismic record.

When the sandstone parameters are model M, the characteristic frequency is about 40Hz, which is close to the main frequency of the Ricker wavelet. In Figure 4(a), it is observed that with the change of the incident angle, there is an obvious reflection coefficient dispersion from 1Hz to 60Hz; as shown in Figure 4(b), in the synthetic seismic record, there is a difference between the low-frequency limit results simulated by the dynamic model of this scheme at 0.3s and the Gassmann equation at 0.7s; as shown in Figure 4(c), the waveform of the dynamic model angle gather stack is distorted and the waveform is no longer symmetrical; as shown in Figure 4(d), the dispersion difference shown at an incident angle of 40° is greater than that at a small angle, and the same phenomenon can be observed in subsequent simulations.

The sandstone parameters are model S, and the characteristic frequency is about 2Hz. As shown in Figure 5(a), the dispersion of the reflection coefficient is no longer obvious after 20Hz; as shown in Figure 5(b), in the synthetic seismic record, the difference between the dynamic model and the low-frequency extreme static model is further increased; as shown in Figure 5(c), the amplitude value of the angle trace of the dynamic model is small, and the amplitude polarity of the near trace and the far trace of the dynamic model is opposite, which cancel each other out after superposition, and the amplitude of the superimposed trace is weak.

The same method is used to calculate the multi-interface seismic dispersion AVO simulation records of sandstone M and mudstone at azimuth angles of 45 and 90 degrees. Comparing the interface of model M at 0, 45 and 90 degrees azimuth, from the change of reflection coefficient, the proportion of dispersion difference at large azimuth (parallel to the fracture direction) is increasing, but from the numerical difference, it is the smallest. This is because the dispersion is strongest when incident perpendicular to the fracture surface, and the dispersion is weakest when parallel to the fracture. For the phenomenon of waveform distortion, it can also be observed that it decreases with the increase of azimuth angle.

Combined with the analysis of synthetic seismic records at different characteristic frequencies, when the characteristic frequency is near the main frequency of the seismic wavelet, the dispersion effect is obvious on the seismic profile, which is mainly reflected in the waveform amplitude difference and waveform distortion. As the characteristic frequency becomes smaller, the difference between the result after considering the dispersion and the result of the static model using the traditional Gassmann equation becomes greater. It should be noted that when the characteristic frequency is approximately one order of magnitude greater than the main frequency of the seismic exploration wavelet, the characteristics of the dispersion on the seismic profile are not obvious, and the dispersion effect can be ignored. However, in terms of well-seismic combination, the role of dispersion has to be considered. In the directional HTI fracture medium, the dispersion phenomenon on the seismic profile gradually weakens with the increase of the azimuth angle. In other words, if you want to observe obvious waveform distortion in the fracture medium, the phenomenon is most obvious in the profile perpendicular to the fracture strike.

The above description is only a specific implementation mode of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can easily think of changes or substitutions within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims (7)

1. A dispersion AVO simulation method for a mesoscopic fracture rock physics model, characterized by comprising a frequency-variable stiffness matrix acquisition stage of the mesoscopic fracture rock and a model dispersion AVO simulation stage; The frequency-dependent stiffness matrix acquisition stage of the mesoscopic fractured rock includes the following steps: A1. Obtain the rock crack density, crack number, crack radius and shear modulus of the background rock, and calculate the normal flexibility and tangential flexibility of the crack; A2. Based on the normal flexibility of the fracture and the P-wave modulus of the saturated fractured rock and background rock, the frequency-dependent normal flexibility of the fracture is calculated. A3. According to the arrangement of the cracks, the dynamic flexibility matrix of the dry rock is calculated based on the tangential flexibility of the cracks, the frequency-dependent normal flexibility of the cracks and the dry rock background flexibility matrix; A3 comprises the following steps: A31. Obtain the arrangement of fractures in the mesoscopic fracture rock physics model; A32, determine whether the arrangement of the cracks is vertical random arrangement, if so, proceed to A33, otherwise proceed to A34; A33, based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, calculate the additional flexibility matrix of the vertical random crack, and use the additional flexibility matrix of the vertical random crack as the frequency-dependent additional flexibility matrix, and enter step A37; The calculation expression of the additional flexibility matrix of the vertical random crack is as follows: Where ΔS(ω) _s represents the additional flexibility matrix of vertical random cracks, Z _N1 (ω) represents the frequency-dependent normal flexibility of cracks when the cracks are arranged vertically randomly, It represents the tangential flexibility of dry cracks when cracks are arranged vertically and randomly; A34, judging whether the vertical arrangement of the cracks is vertical directional arrangement, if so, proceeding to A35, otherwise proceeding to A36; A35, based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, calculate the additional flexibility matrix of the vertical oriented crack, and use the additional flexibility matrix of the vertical oriented crack as the frequency-dependent additional flexibility matrix, and enter step A37; The calculation expression of the additional flexibility matrix of the vertical oriented crack is as follows: Where ΔS(ω) _d represents the dry rock additional flexibility matrix for fracture orientation, _ZN2 (ω) represents the frequency-dependent fracture normal flexibility when the fracture is vertically oriented, It represents the tangential flexibility of dry cracks when the cracks are arranged vertically; A36, judging whether the vertical arrangement of the cracks is vertically oriented and randomly arranged, if so, proceeding to A37, otherwise proceeding to A38; A37. Based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, a directional and random additional flexibility matrix of the vertical crack is calculated, and the directional and random additional flexibility matrix of the vertical crack is used as the frequency-dependent additional flexibility matrix; The calculation expression of the additional flexibility matrix of the vertical crack orientation and randomness is as follows: Where, ΔS(ω) _sd represents the additional flexibility matrix of the vertical crack which is oriented and random; A38. Based on the tangential flexibility of the crack and the frequency-dependent normal flexibility of the crack, a completely random additional flexibility matrix of the crack is calculated, and the completely random additional flexibility matrix of the crack is used as the frequency-dependent additional flexibility matrix; The calculation expression of the additional flexibility matrix of the completely random crack is as follows: Where ΔS(ω) ₃ represents the additional flexibility matrix when the cracks are completely random, Z _N (ω) ₃ represents the frequency-dependent normal flexibility of the cracks when the cracks are completely randomly arranged, It represents the tangential flexibility of dry cracks when the cracks are arranged completely randomly; A39, based on the dry rock background flexibility matrix and the frequency-dependent additional flexibility matrix, the dry rock dynamic flexibility matrix is calculated; The calculation expression of the dry rock dynamic flexibility matrix is as follows: Where S ⁰ (ω) represents the dynamic flexibility matrix of dry rock, represents the dry rock background flexibility matrix, ΔS(ω) represents the frequency-dependent flexibility matrix; A4. Based on the dry rock dynamic flexibility matrix and the anisotropic Gassmann equation, the complete rock frequency-dependent flexibility matrix is calculated; A5. Invert the complete rock frequency-variant flexibility matrix to obtain the rock frequency-variant stiffness matrix; The model dispersion AVO simulation stage includes the following steps: B1. Based on the frequency-dependent stiffness matrix of rock, the frequency-dependent stiffness coefficient is obtained; B2. Based on the real part of the frequency-dependent stiffness coefficient, the frequency-dependent anisotropy coefficient is calculated; B3. Based on the anisotropy coefficient and reflection coefficient equation, the frequency-dependent interface reflection coefficient is obtained; B4. Based on the frequency-varying interface reflection coefficient and source wavelet, the multi-interface seismic dispersion AVO simulation record is calculated. 2. The dispersion AVO simulation method of the mesoscopic fracture rock physics model according to claim 1, characterized in that the calculation expressions of the normal flexibility and the tangential flexibility in A1 are as follows: in, and They represent the normal flexibility and tangential flexibility of the dry crack, λ represents the Lame coefficient, μ represents the shear modulus of the background rock, e represents the crack density, N represents the number of cracks, a represents the crack radius, and V represents the rock volume. 3. The dispersion AVO simulation method of the mesoscopic fracture rock physics model according to claim 2, characterized in that said A2 comprises the following steps: A21. Based on the P-wave modulus of saturated fractured rock and background rock, the time scale parameter and shape parameter are calculated; The calculation expressions of the time scale parameter and shape parameter are as follows: Among them, τ and denote the time scale parameter and shape parameter respectively, _Cb denotes the P-wave modulus of saturated fractured rock and background rock under high-frequency limit conditions, _C0 denotes the P-wave modulus of saturated fractured rock and background rock under low-frequency limit conditions, G denotes the low-frequency scale parameter, and T denotes the high-frequency scale parameter; A22, based on the time scale parameter and shape parameter, the flexibility frequency relationship is calculated; The calculation expression of the flexibility frequency relationship is as follows: Where f _fra (ω) represents the flexibility frequency relationship, i represents the imaginary part of the complex number, and ω represents the frequency; A23, based on the flexibility frequency relationship and the normal flexibility of the crack, the frequency-dependent normal flexibility of the crack is calculated; The calculation expression of the frequency-dependent crack normal flexibility is as follows: where Z _N (ω) represents the frequency-dependent crack normal compliance. 4. The dispersion AVO simulation method of the mesoscopic fracture rock physics model according to claim 1, characterized in that the calculation expression of the complete rock frequency-dependent flexibility matrix in A4 is as follows: in, and denote the flexibility tensors of saturated rock, dry rock and mineral particles, respectively, represents the sum of the flexibility of the i′th row in the matrix corresponding to the dry rock flexibility tensor, represents the sum of the flexibility of the i′th row in the matrix corresponding to the flexibility tensor of the mineral particles, represents the sum of the flexibility of the j′th column in the matrix corresponding to the dry rock flexibility tensor, represents the sum of the flexibility of the j′th column in the matrix corresponding to the flexibility tensor of the mineral particles, β _dry represents the generalized compressibility coefficient of dry rock with cracks, β _g represents the compressibility coefficient of mineral particles, β _f represents the compressibility coefficient of fluid, Represents the total porosity, where i′=1,2,3,4,5,6 and j′=1,2,3,4,5,6. 5. The dispersion AVO simulation method of the mesoscopic fracture rock physics model according to claim 4, characterized in that the calculation expression of the anisotropy coefficient in B2 is as follows: Wherein, α(ω) represents the first frequency-dependent anisotropy coefficient, c ₃₃ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 3rd row and 3rd column in the frequency-dependent stiffness matrix, ρ represents the medium density, β(ω) represents the second frequency-dependent anisotropy coefficient, c ₄₄ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 4th row and 4th column in the frequency-dependent stiffness matrix, ε(ω) represents the third frequency-dependent anisotropy coefficient, c ₁₁ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 1st column in the frequency-dependent stiffness matrix, γ(ω) represents the fourth frequency-dependent anisotropy coefficient, c ₆₆ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 6th row and 6th column in the frequency-dependent stiffness matrix, δ(ω) represents the fifth frequency-dependent anisotropy coefficient, c ₁₃ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 3rd column in the frequency-dependent stiffness matrix, c ₅₅ (ω) represents the real part of the frequency-dependent stiffness coefficient at the 1st row and 3rd column in the frequency-dependent stiffness matrix, (ω) represents the real part of the frequency-dependent stiffness coefficient at the 5th row and 5th column in the frequency-dependent stiffness matrix. 6. The dispersion AVO simulation method of the mesoscopic fracture rock physics model according to claim 5, characterized in that the calculation expression of the frequency-dependent interface reflection coefficient in B3 is as follows: Δγ(ω)=γ ₂ (ω)-γ ₁ (ω) Δα(ω)＝α ₂ (ω)-α ₁ (ω) Δε(ω)＝ε ₂ (ω)-ε ₁ (ω) Δδ(ω)＝δ ₂ (ω)-δ ₁ (ω) ΔZ(ω)＝Z ₂ (ω)-Z ₁ (ω) Z(ω)＝ρα(ω) Δμ(ω)=μ ₂ (ω)-μ ₁ (ω) μ(ω)=ρβ(ω) ² Where R _pp (ω,ξ,θ) represents the frequency-dependent interface reflection coefficient, ΔZ(ω) represents the frequency-dependent vertical P-wave impedance difference between the lower and upper layers, represents the average value of the vertical P-wave impedance of the lower and upper layers related to frequency, Δα(ω) represents the difference of the first anisotropy coefficient of the lower and upper layers dependent on frequency, represents the average value of the first anisotropy coefficient of the lower and upper layers in terms of frequency dependence, represents the average value of the second anisotropy coefficient of the lower and upper layers in terms of frequency, Δμ(ω) represents the difference in the vertical shear wave modulus of the lower and upper layers in terms of frequency, represents the average value of the vertical shear wave modulus of the lower and upper media related to frequency, Δδ(ω) represents the difference of the fifth anisotropy coefficient of the lower and upper media depending on frequency, θ represents the P-wave incident angle, Δε(ω) represents the difference of the third anisotropy coefficient of the lower and upper media depending on frequency, ξ represents the P-wave incident azimuth, α ₁ (ω) and α ₂ (ω) represent the first anisotropy coefficient of the upper and lower media depending on frequency, β ₁ (ω) and β ₂ (ω) represent the second anisotropy coefficient of the upper and lower media depending on frequency, ε ₁ (ω) and ε ₂ (ω) represent the third anisotropy coefficient of the upper and lower media depending on frequency, δ ₁ (ω) and δ ₂ (ω) represent the fifth anisotropy coefficient of the upper and lower media depending on frequency, γ ₁ (ω) and γ ₂ (ω) (ω) represents the fourth anisotropy coefficient of the upper medium and the lower medium in dependence of frequency, Δγ(ω) represents the difference of the fourth anisotropy coefficient of the lower medium and the upper medium in dependence of frequency, Z ₁ (ω) and Z ₂ (ω) represent the vertical P-wave impedance of the lower medium and the upper medium in dependence of frequency, Z(ω) represents the vertical P-wave impedance in dependence of frequency, μ ₁ and μ ₂ represent the vertical shear wave modulus of the lower medium and the upper medium in dependence of frequency, respectively. 7. The dispersion AVO simulation method of the mesoscopic fracture rock physics model according to claim 6, characterized in that said B4 comprises the following steps: B41. Based on the frequency-dependent interface reflection coefficient, the frequency-dependent reflection coefficient of the interface is obtained; B42, based on the interface-frequency-dependent reflection coefficient and source wavelet, the multi-interface seismic dispersion AVO simulation record is calculated; The calculation expression of the multi-interface seismic dispersion AVO simulation record is as follows: where x _t (θ) is the earthquake record at time depth t, represents the unit reflection coefficient of the nth interface, k represents the total number of interface layers, Δt represents the sampling interval, W(ω) represents the source wavelet, R ⁿ (θ,ω) represents the frequency-dependent reflection coefficient of the nth interface, * represents the convolution operation, and FFT ^-1 represents the inverse Fourier transform.