A method for azimuthal seismic inversion in biaxial stress-induced anisotropic media

Abstract

The present invention relates to the field of oil and gas exploration and development, and in particular to a method for azimuthal seismic inversion of biaxial stress-induced anisotropic media. The method comprises the following steps: step S1, parameterization of the stiffness tensor of biaxial stress-induced anisotropic media; step S2, derivation of the PP wave reflection coefficient equation of biaxial stress-induced anisotropic media; step S3, azimuthal seismic inversion of biaxial stress-induced anisotropic media. The present invention utilizes the amplitude difference of seismic reflection data in different azimuths to invert the model parameters of the anisotropic part, and then uses the inverted model parameters of the anisotropic part as input to further invert the model parameters of the isotropic part. The method effectively reduces the pathological nature of multi-parameter synchronous inversion and improves the stability and accuracy of the inversion.

Description

Technical Field

[0001] The invention relates to the field of oil and gas exploration and development, and in particular to a biaxial stress-induced anisotropic medium azimuthal seismic inversion method. Background Art

[0002] Seismic anisotropy reveals the directional dependence of seismic wave velocity and amplitude, which is of great significance for geophysical exploration and development. The physical mechanisms that cause anisotropy are diverse, including preferentially oriented crystals or minerals, fine stratification or arrangement of clay particles, accumulation of strata with different elastic properties, directional fractures, and geostress. Like intrinsic anisotropy and fracture-induced anisotropy, stress-induced anisotropy is also widely present in the underground, but it is often overlooked in seismic anisotropy studies.

[0003] At present, the PP wave reflection coefficient equation of stress-induced anisotropic media is usually derived based on the acoustoelastic theory, that is, the rock is regarded as a continuous medium, and the relationship between the inherent stress state and the elastic properties is explored from a macroscopic perspective. For example, the technical solution of the patent document "A method and system for azimuthal seismic inversion of stress-horizontally anisotropic media" (application publication number: CN115840250A) is: the first step is to establish a transversely isotropic medium (HTI medium) model with a horizontal symmetry axis under linear slip boundary conditions and apply a uniaxial small stress parallel to the symmetry axis to the HTI medium model; the second step is to decouple the stiffness tensor induced by cracks and stress in the horizontal uniaxial stress HTI medium; the third step is to establish the PP wave approximate equation for decoupling cracks and stress-induced anisotropic parameters; the fourth step is to use the difference in reflection coefficient amplitude between different azimuths for azimuthal seismic inversion. This method focuses on HTI media under horizontal stress and uses the solid acoustoelastic theory to study the influence of stress on rock properties.

[0004] However, the stress dependence of rock elastic properties can be explained as the response of microcracks in rocks to the applied stress field, that is, when isotropic rocks are subjected to anisotropic stress, microcracks perpendicular to the stress direction will be closed first, thereby inducing elastic anisotropy of the rock. However, the acoustoelastic theory cannot explain the mechanism of stress-induced elastic anisotropy from a microscopic perspective. In addition, the PP wave reflection coefficient of anisotropic media contains many unknown parameters, and the numerical contributions of different parameters to the reflection coefficient are significantly different.

[0005] Therefore, the existing multi-parameter synchronous inversion of anisotropic media often has problems such as low parameter inversion stability and poor noise resistance. Summary of the invention

[0006] The present invention takes into account the effect of stress on microcrack closure, and based on this derives a PP wave reflection coefficient equation for coupling stress and rock microstructure changes, and proposes an azimuthal seismic step-by-step inversion method based on azimuthal amplitude differences under a Bayesian framework to overcome the above-mentioned defects of the prior art.

[0007] In order to solve the above technical problems, the technical solution adopted by the present invention is: a method for azimuthal seismic inversion of biaxial stress-induced anisotropic media, which comprises the following steps:

[0008] Step S1, parameterization of the stiffness tensor of biaxial stress-induced anisotropic media;

[0009] Step S2, deriving a PP wave reflection coefficient equation for a biaxial stress-induced anisotropic medium;

[0010] Step S3: azimuthal seismic inversion of biaxial stress-induced anisotropic media.

[0011] Preferably, the step S1 comprises:

[0012] Considering that the unstressed rock is elastically isotropic and full of randomly distributed and randomly oriented microcracks, the elastic stiffness tensor of the unstressed rock is expressed as:

[0013] (1),

[0014] In formula (1), and represent the longitudinal wave modulus and shear modulus of unstressed rock respectively;

[0015] The equivalent flexibility tensor caused by the closure of microcracks due to horizontal biaxial stress is expressed as:

[0016] (2),

[0017] In formula (2), represents the equivalent flexibility matrix of the rock without microcracks, and They represent the normal flexibility and tangential flexibility of microcracks in unstressed rock, respectively. The normalized stress is expressed as and , The tensor represents the typical closure pressure of the fracture. and represents the compliance contribution of the initial microcracks, the tensor and Represents stress The flexibility contribution that causes crack closure, tensor and Represents stress The tensor expressions of the flexibility contribution that causes crack closure are:

[0018] (3),

[0019] (4),

[0020] (5),

[0021] Rearranging formula (2) yields:

[0022] (6),

[0023] In formula (6), represents the flexibility tensor of unstressed rock, and Represents horizontal stress , Compliance contribution that causes microcrack closure;

[0024] Assuming biaxial stress-induced anisotropic media, the equivalent stiffness tensor is derived as:

[0025] (7),

[0026] Combining formulas (1)-(7), the stiffness tensor of biaxial stress-induced anisotropic media is derived as follows:

[0027] (8),

[0028] The biaxial stress-induced anisotropic medium is approximately equivalent to the orthotropic anisotropic medium, and the stress-induced anisotropic parameter is defined as:

[0029] (9),

[0030] (10),

[0031] The components of the elastic stiffness tensor in formula (8) are expressed as:

[0032] (11),

[0033] In formula (11), , , .

[0034] Preferably, the step S2 comprises:

[0035] The PP wave reflection coefficient equation is expressed as:

[0036] (12),

[0037] In formula (12), represents the density of the background medium, is the angle of incidence, represents the perturbation of the stiffness tensor components, The expression is as follows:

[0038] (13),

[0039] In formula (13), represents the observation azimuth, represents the longitudinal wave velocity of the background medium;

[0040] Consider the interface of two biaxial stress-induced anisotropic media with different properties, and derive the perturbation terms of each component of the stiffness tensor by combining formula (11), ignoring the , , , , , , , , , , , Proportional terms, we get:

[0041] (14),

[0042] In formula (14), △ represents the difference in the properties of the media above and below the reflection interface;

[0043] Combining formulas (12)-(14), the PP wave reflection coefficient equation of biaxial stress-induced anisotropic media is derived as:

[0044] (15),

[0045] In formula (15),

[0046] (16),

[0047] (17),

[0048] (18),

[0049] In formulas (16)-(18), , , ,

[0050] ,

[0051] ,

[0052] ,

[0053] .

[0054] Preferably, step S3 comprises:

[0055] Step S31: According to formula (15), the azimuth difference reflection coefficient is expressed as:

[0056] (19),

[0057] In formula (19), and Represent two different observation azimuths, and,

[0058] (20),

[0059] Considering the case of M reflection interfaces and L incident angles, based on the convolution model, the azimuth difference of seismic data is expressed in matrix form as follows:

[0060] (twenty one),

[0061] In formula (21),

[0062] ,

[0063] The incident angle is And the observation azimuth is The seismic traces;

[0064] ,

[0065] ,

[0066] ,

[0067] , , ,

[0068] , ,

[0069] Where W represents the wavelet matrix, represents the Mth reflection interface, and Respectively represent matrix transposition and generate an M×M diagonal matrix;

[0070] In the Bayesian theory framework, the posterior probability distribution of the model parameters is Proportional to the likelihood function of the seismic data and the prior probability distribution of the model parameters :

[0071] (twenty two),

[0072] In formula (22), Assuming it obeys Cauchy distribution, it can be expressed as:

[0073] (twenty three),

[0074] In formula (23), represents the variance of the model parameters;

[0075] Assuming that the likelihood function follows a Gaussian distribution, it can be expressed as:

[0076] (twenty four),

[0077] In formula (24), represents the variance of noise in seismic data;

[0078] Combining formulas (22)-(24), the target functional of step S31 is constructed as:

[0079] (25),

[0080] In formula (25), represents the regularization coefficient of the model parameters, Represents the low-frequency components of the model parameters;

[0081] The iterative weighted least squares algorithm is used to optimize and solve formula (25). After multiple iterations, the unknown model parameters are obtained;

[0082] Step S32: The model parameters of the anisotropic part obtained by inversion in step S31 are , as input, invert the model parameters of the isotropic part;

[0083] The seismic forward modeling equation in step S32 is expressed as:

[0084] (26),

[0085] In formula (26),

[0086] ,

[0087] ,

[0088] ,

[0089] ,

[0090] In step S32, the forward equation (21) in step S31 is replaced by equation (26), and the target functional of step S32 is constructed using the same inversion method as step S31:

[0091] (27),

[0092] In formula (27), and Seismic data The variance of the noise in and the variance of the isotropic part of the model parameters, represents the regularization coefficient of the isotropic part of the model parameters, The low-frequency components of the model parameters of the isotropic part are represented. Finally, the iterative weighted least squares algorithm is used to solve the target functional of step S32 to obtain the model parameters M, μ, and ρ of the isotropic part.

[0093] Beneficial technical effects brought by the present invention:

[0094] The present invention utilizes the amplitude difference of seismic reflection data in different azimuths to invert the model parameters of the anisotropic part, and then uses the inverted model parameters of the anisotropic part as input to further invert the model parameters of the isotropic part. This method effectively reduces the pathological nature of multi-parameter synchronous inversion and improves the stability and accuracy of the inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0095] Figure 1 This is a flow chart of an embodiment of the present invention.

[0096] Figure 2 Schematic diagram of the biaxial stress-induced anisotropic medium model.

[0097] Figure 3 Schematic diagram of the reflection coefficient curve changing with the observation azimuth.

[0098] Figure 4 Schematic diagram of the reflection coefficient changing with the incident angle.

[0099] Figure 5 is the model parameter Schematic diagram of the inversion results.

[0100] Figure 6 is the model parameter Schematic diagram of the inversion results.

[0101] Figure 7 is the model parameter Schematic diagram of the inversion results.

[0102] Figure 8 is the model parameter Schematic diagram of the inversion results.

[0103] Fig. 9 is the model parameter Schematic diagram of the inversion results.

[0104] Fig.10 is the model parameter Schematic diagram of the inversion results.

[0105] Fig.11 is the model parameter Schematic diagram of the inversion results. DETAILED DESCRIPTION

[0106] The detailed description and technical contents of the present invention are described as follows with reference to the accompanying drawings. However, the accompanying drawings are only provided for reference and description and are not intended to limit the present invention.

[0107] like Figure 1 As shown, an embodiment of the present invention provides a method for azimuthal seismic inversion of biaxial stress-induced anisotropic media, comprising the following steps:

[0108] Step S1, parameterization of the stiffness tensor of biaxial stress-induced anisotropic media;

[0109] Step S2, deriving a PP wave reflection coefficient equation for a biaxial stress-induced anisotropic medium;

[0110] Step S3: azimuthal seismic inversion of biaxial stress-induced anisotropic media.

[0111] Step S1, the specific steps of parameterizing the stiffness tensor of biaxial stress-induced anisotropic media are as follows.

[0112] Consider that the unstressed rock is elastically isotropic and full of randomly distributed and randomly oriented microcracks. According to the linear Hooke's law, the elastic stiffness tensor of the unstressed rock can be expressed as:

[0113] (1),

[0114] In formula (1), and represent the longitudinal wave modulus and shear modulus of unstressed rock respectively.

[0115] like Figure 2 As shown, when isotropic rock is subjected to horizontal biaxial stress and When the stress acts, the microcracks perpendicular to the stress direction will close first, causing elastic anisotropy of the rock. Among them, ψ represents the angle between the projection of the seismic wave propagation direction on the horizontal plane and the X1 axis, which is defined as the observation azimuth; θ represents the angle between the incident seismic wave propagation direction and the X3 axis, which is defined as the seismic wave incident angle.

[0116] According to Collet et al. (Geophysics, 2014, 79(5): C135-C145, doi: 10.1190 / geo2013-0311.1.), the equivalent flexibility tensor caused by microcrack closure due to horizontal biaxial stress can be expressed as:

[0117] (2),

[0118] In formula (2), Represents the equivalent flexibility matrix of rock without microcracks. and They represent the normal and tangential flexibility of microcracks in unstressed rock. The normalized stress can be expressed as and ,in, represents the typical closure pressure of the fracture. and represents the compliance contribution of the initial microcrack, the tensor and Represents stress The flexibility contribution that causes crack closure, tensor and Represents stress The flexibility contribution that causes crack closure is expressed as follows:

[0119] (3),

[0120] (4),

[0121] (5),

[0122] In order to highlight the flexibility contribution of different parts, formula (2) is rearranged to obtain:

[0123] (6),

[0124] In formula (6), represents the flexibility tensor of unstressed rock, and Represents horizontal stress , Compliance contribution leading to microcrack closure.

[0125] Furthermore, assuming that the biaxial horizontal stress is small and the anisotropy induced by the stress is weaker than that of the unstressed rock, the equivalent stiffness tensor of the biaxial stress-induced anisotropic medium can be derived as:

[0126] (7),

[0127] Combining formulas (1)-(7), the stiffness tensor of biaxial stress-induced anisotropic media can be derived as:

[0128] (8),

[0129] From the symmetry of the stiffness tensor, it can be seen that the biaxial stress-induced anisotropic medium is approximately equivalent to the orthogonal anisotropic medium. In order to simplify the stiffness tensor, the stress-induced anisotropic parameters are defined as follows:

[0130] (9),

[0131] (10),

[0132] Finally, the components of the elastic stiffness tensor in formula (8) can be expressed as:

[0133] (11),

[0134] In formula (11), , , .

[0135] Step S2, the specific steps of deriving the PP wave reflection coefficient equation of the biaxial stress-induced anisotropic medium are as follows.

[0136] According to Born scattering theory, the PP wave reflection coefficient equation can be expressed as:

[0137] (12),

[0138] In formula (12), represents the density of the background medium, is the angle of incidence, represents the perturbation of the stiffness tensor components, The specific expression is as follows:

[0139] (13),

[0140] In formula (13), represents the observation azimuth, Represents the longitudinal wave velocity of the background medium.

[0141] Consider the interface of two biaxial stress-induced anisotropic media with different properties. Combined with formula (11), the perturbation terms of each component of the stiffness tensor are derived. At the same time, assuming that the difference in the properties of the media on both sides of the reflection interface is small, the influence of , , , , , , , , , , , Proportional terms, we can get:

[0142] (14),

[0143] In formula (14), △ represents the difference in the properties of the media above and below the reflection interface.

[0144] Combining formulas (12)-(14), the PP wave reflection coefficient equation of biaxial stress-induced anisotropic media can be derived as follows:

[0145] (15),

[0146] In formula (15),

[0147] (16),

[0148] (17),

[0149] (18),

[0150] In formulas (16)-(18), , , ,

[0151] ,

[0152] ,

[0153] ,

[0154] .

[0155] Formula (15) is compared with the reflection coefficient equation for arbitrary weakly anisotropic media proposed by Pšenčík and Martins to verify the accuracy of the equation, as follows: Figure 3 and Figure 4 shown. Figure 3The curves showing the variation of the reflection coefficient with the observation azimuth angle are shown. The dotted line represents the reflection coefficient curve calculated using the Pšenčík and Martins equations, and the solid line represents the reflection coefficient curve calculated using formula (15). Each curve in the figure corresponds to the incident angle θ = 0°, 10°, and 30° from top to bottom. Figure 4 The curves showing the change of reflection coefficient with incident angle are shown. The dotted line represents the reflection coefficient curve calculated using the Pšenčík and Martins equations, and the solid line represents the reflection coefficient curve calculated using formula (15). Each curve in the figure corresponds to the observation azimuth angle ψ = 90°, 45°, and 0° from top to bottom. Figure 3 and Figure 4 It can be seen that the PP wave reflection coefficient equation (15) of the biaxial stress-induced anisotropic medium proposed in the present invention has a high degree of consistency with the Pšenčík and Martins equation, which proves the rationality of the derived equation.

[0156] Step S3, the specific steps of azimuthal seismic inversion of biaxial stress-induced anisotropic media are as follows.

[0157] Based on the azimuthal variation characteristics of the reflection coefficient of stress-induced anisotropic media, the model parameters of the anisotropic part are first inverted using the azimuthal amplitude difference method based on the Bayesian framework, considering the Cauchy prior information and low-frequency regularization constraints. Then, the inversion results are used as input to invert the model parameters of the isotropic part.

[0158] Step S31: According to formula (15), the azimuth difference reflection coefficient can be expressed as:

[0159] (19),

[0160] In formula (19), and Represent two different observation azimuths, and,

[0161] (20),

[0162] Considering the case of M reflection interfaces and L incident angles, based on the convolution model, the azimuth difference of seismic data can be expressed in matrix form as follows:

[0163] (twenty one),

[0164] In formula (21),

[0165] ,

[0166] The incident angle is And the observation azimuth is The seismic traces;

[0167] ,

[0168] ,

[0169] ,

[0170] , , ,

[0171] , ,

[0172] Where W represents the wavelet matrix, represents the Mth reflection interface. and They represent matrix transposition and generation of an M×M diagonal matrix respectively.

[0173] In the Bayesian theory framework, the posterior probability distribution of the model parameters is Proportional to the likelihood function of the seismic data and the prior probability distribution of the model parameters ,Right now

[0174] (twenty two),

[0175] In formula (22), It is usually assumed to obey the Cauchy distribution, expressed as

[0176] (twenty three),

[0177] In formula (23), Represents the variance of the model parameters.

[0178] Assuming that the likelihood function follows a Gaussian distribution, it can be expressed as:

[0179] (twenty four),

[0180] In formula (24), represents the variance of the noise in the seismic data. The background constraints of the model parameters can effectively improve the stability of the inversion. Therefore, combined with formulas (22)-(24), the target functional can be constructed as:

[0181] (25),

[0182] In formula (25), represents the regularization coefficient of the model parameters, Represents the low-frequency components of the model parameters.

[0183] The iterative weighted least squares algorithm is used to optimize and solve formula (25). After multiple iterations, the unknown model parameters can be obtained.

[0184] Step S32: The model parameters of the anisotropic part obtained by inversion in step S31 (assuming that it is ) as input to invert the model parameters of the isotropic part.

[0185] The seismic forward modeling equation of step S32 can be expressed as:

[0186] (26),

[0187] In formula (26),

[0188] ,

[0189] ,

[0190] ,

[0191] ,

[0192] In the inversion of step S32, the forward equation (21) in the inversion of step S31 is replaced by equation (26), and the inversion target functional of step S32 can be constructed using the same inversion method as step S31, and the following is obtained:

[0193] (27),

[0194] in, and Seismic data The variance of the noise in and the variance of the isotropic part of the model parameters, represents the regularization coefficient of the isotropic part of the model parameters, The low-frequency component of the isotropic model parameters is represented. Finally, the iterative weighted least squares algorithm is used to solve the target functional to obtain the isotropic model parameters M, μ, and ρ.

[0195] Case study testing the inversion method:

[0196] An oil and gas reservoir in southwest China is taken as an example to verify the inversion method proposed in the present invention. The oil and gas reservoir is a tight shale reservoir with low porosity and low permeability, and the structural fracture system is an important channel for oil and gas seepage. In the target work area, it is assumed that the oil-bearing reservoir can be regarded as an orthogonal medium induced by horizontal biaxial stress. Before being used in the present invention, the seismic traces collected on-site are preprocessed according to the data standard of azimuthal seismic inversion, including denoising, amplitude compensation, rotation angle domain, etc., and azimuthal seismic data are generated. Then, azimuthal seismic inversion is carried out to estimate stress-induced anisotropy parameters and elastic parameters. The results are as follows. Figure 5-Figure 11 As shown. It can be found that in the oil-bearing target reservoir (the location circled by the black ellipse), the inverted P-wave modulus, shear modulus and density all show abnormally low values ​​( Figure 5-Figure 7 ). This indicates that fluids reduce the elastic modulus of rocks, which is consistent with the prediction of Gassmann's fluid substitution formula, and indirectly reflects the rationality of the inversion results of the isotropic model parameters. The inversion results of stress-induced anisotropic parameters are obviously higher ( Figure 8-Figure 11 ), indicating that there is obvious stress-induced anisotropy in the target formation. In addition, the correlation between the seven model parameters confirms the reliability of the inversion results.

[0197] The above embodiments are only preferred embodiments for fully illustrating the present invention, and the protection scope of the present invention is not limited thereto. Any equivalent substitution or change made by a person skilled in the art based on the present invention is within the protection scope of the present invention.

Claims

1. A method for azimuthal seismic inversion of biaxial stress-induced anisotropic media, characterized in that: It includes the following steps: Step S1, parameterization of the stiffness tensor of biaxial stress-induced anisotropic media; Step S2, deriving a PP wave reflection coefficient equation for a biaxial stress-induced anisotropic medium; PP wave reflection coefficient equation for biaxial stress-induced anisotropic media for: （15）， In formula (15), is the angle of incidence, represents the observation azimuth, （16）， （17）， （18）， In formulas (16)-(18), and They represent the longitudinal wave modulus and shear modulus of unstressed rock, respectively. ΔM represents the difference in longitudinal wave modulus of unstressed rock above and below the reflection interface. Δμ represents the difference in shear modulus of unstressed rock above and below the reflection interface. represents the density of the background medium, represents the difference in density of the media above and below the reflection interface, , , , ， ， ， , in, , represents the ratio of the unstressed rock shear modulus to the unstressed rock longitudinal wave modulus; Step S3: azimuthal seismic inversion of biaxial stress-induced anisotropic media.

2. A method for azimuthal seismic inversion of biaxial stress-induced anisotropic media according to claim 1, characterized in that: The step S1 comprises: Considering that the unstressed rock is elastically isotropic and filled with randomly distributed and oriented microcracks, the elastic stiffness tensor of the unstressed rock is It is expressed as: （1）， In formula (1), and represent the longitudinal wave modulus and shear modulus of unstressed rock respectively; Equivalent flexibility tensor caused by microcrack closure due to horizontal biaxial stress It is expressed as: （2）， In formula (2), represents the equivalent flexibility matrix of the rock without microcracks, and They represent the normal flexibility and tangential flexibility of microcracks in unstressed rock, respectively. The normalized stress is expressed as and , The tensor represents the typical closure pressure of the fracture, and They represent the tangential compliance contribution rate and compliance difference of the initial microcracks respectively. Contribution rate, tensor and Represents horizontal stress Tangential compliance contribution and compliance difference caused by crack closure Contribution rate, tensor and Represents horizontal stress Tangential compliance contribution and compliance difference caused by crack closure Contribution rate, the tensor expressions are: （3）， （4）， （5）， Rearranging formula (2) yields: （6）， In formula (6), represents the flexibility tensor of unstressed rock, and Represents horizontal stress , The compliance contribution that causes microcrack closure; Assuming biaxial stress-induced anisotropic media, the equivalent stiffness tensor It is derived as follows: （7）， Combining formulas (1)-(7), the equivalent stiffness tensor of biaxial stress-induced anisotropic media is derived as follows: （8）， The biaxial stress-induced anisotropic medium is approximately equivalent to the orthotropic anisotropic medium, and the stress-induced anisotropic parameter is defined as: （9）， （10）， The components of the elastic stiffness tensor in formula (8) are expressed as: （11）， In formula (11), , represents the first Lamé constant, , which represents the ratio of the first Lame constant to the longitudinal wave modulus of unstressed rock.

3. A method for azimuthal seismic inversion of biaxial stress-induced anisotropic media as claimed in claim 2, characterized in that: The step S2 comprises: The PP wave reflection coefficient equation It is expressed as: （12）， In formula (12), represents the density of the background medium, represents the difference in density of the media above and below the reflection interface, is the angle of incidence, Represents the stiffness matrix of the media above and below the reflection interface The difference is used to represent the perturbation of the stiffness tensor component, i and j are the matrix orders, The expression is as follows: （13）， In formula (13), represents the observation azimuth, represents the longitudinal wave velocity of the background medium; Consider the interface of two biaxial stress-induced anisotropic media with different properties, and derive the perturbation terms of each component of the stiffness tensor by combining formula (11), ignoring the , , , , , , , , , , , Proportional terms, we get: （14）， In formula (14), △ represents the difference in the properties of the media above and below the reflection interface; Combining formulas (12)-(14), the PP wave reflection coefficient equation for biaxial stress-induced anisotropic media is derived: .

4. A method for azimuthal seismic inversion of biaxial stress-induced anisotropic media as claimed in claim 3, characterized in that: The step S3 comprises: Step S31: According to formula (15), the azimuth difference reflection coefficient It is expressed as: （19）， In formula (19), and Represent two different observation azimuths, and, （20）， Considering the case of B reflection interfaces and L incident angles, based on the convolution model, the azimuth difference of seismic data is expressed in matrix form as follows: （21）， In formula (21), ， The incident angle is And the observation azimuth is The earthquake channel, The incident angle is And the observation azimuth is The earthquake channel, ， ， ， ， ， ， ， ， Where W represents the wavelet matrix, and Respectively represent matrix transposition and generate a B×B diagonal matrix; In the Bayesian theory framework, the posterior probability distribution of the model parameters is Proportional to the likelihood function of the seismic data and the prior probability distribution of the model parameters : （22）， In formula (22), Assuming it obeys Cauchy distribution, it can be expressed as: （23）， In formula (23), represents the standard deviation of the model parameters, represents the variance of the model parameters; Assuming that the likelihood function follows a Gaussian distribution, it can be expressed as: （24）， In formula (24), represents the standard deviation of the noise in the seismic data, represents the variance of noise in seismic data; Combining formulas (22)-(24), the target functional of step S31 is constructed: for: （25）， In formula (25), represents the regularization coefficient of the model parameters, Represents the low-frequency components of the model parameters; The iterative weighted least squares algorithm is used to optimize and solve formula (25). After multiple iterations, the unknown model parameters are obtained; Step S32: using the model parameters of the anisotropic part obtained by inversion in step S31 as input, inverting the model parameters of the isotropic part.

Images