Anisotropic medium multi-parameter full waveform inversion method including uncertainty analysis

Abstract

This invention presents a multi-parameter full-waveform inversion method for anisotropic media incorporating uncertainty analysis. The method includes: establishing an initial model for each independent elastic modulus in a transversely anisotropic medium with a tilted axis of symmetry; performing forward numerical simulation using the initial model of elastic modulus to obtain a synthetic seismic dataset; calculating the data residuals between observed and simulated seismic data; deriving the maximum a posteriori solution based on a Bayesian probabilistic system; and iteratively updating the model parameters and posterior covariance matrix for the next frequency from low to high frequencies using an iterative extended Kalman filter, with the inversion result of the previous frequency serving as the initial model for the next frequency, until the calculated frequency value reaches a set maximum frequency value and the final data residuals meet the accuracy requirements. This invention solves the existing problems of multi-parameter inversion and quantitative estimation of uncertainties in elastic anisotropic media.

Description

Technical Field

[0001] This invention belongs to the field of exploration geophysics, specifically involving a multi-parameter full waveform inversion method for anisotropic media that incorporates uncertainty analysis. Background Technology

[0002] With advancements in geophysics and continuous improvements in exploration capabilities, current seismic reservoir prediction focuses on solving more complex geophysical problems. In common sedimentary strata and complex rock structures, the propagation velocity of seismic waves changes with the direction of propagation. Therefore, by fully considering the anisotropic characteristics of the subsurface medium, more accurate inversion results can be obtained. However, as the number of parameters increases, the difficulty of parameter inversion increases, and crosstalk between different parameters can significantly affect the inversion results. According to Tarantola, A., 2005, Inverse problem theory and methods for model parameter estimation: SIAM, 89, the process of determining geophysical parameters from seismic reflection data can be understood as an inverse scattering problem, and scattering theory is essentially an analysis of disturbed media. Seismic scattering theory has significant advantages in dealing with localized inhomogeneities with shorter seismic wavelengths, such as small-scale media containing pores and fractures, as well as strong scatterers. Applying seismic scattering theory to full-waveform inversion can more accurately characterize the physical parameter information of complex anisotropic media.

[0003] The full-waveform inversion method fully utilizes the kinematic and dynamic characteristics of seismic waves to obtain subsurface model parameter information, offering advantages such as high imaging accuracy of complex structures and good inversion results of physical parameters (Lailly, P., 1983, Theseismic inverse problem as a sequence of before stack migrations: Conference on Inverse Scattering: Theory and Application, Society for Industrial and Applied Mathematics, Expanded Abstracts, 206–220.). However, in actual industrial production, the inaccuracy of the initial model, the lack of low-frequency information, the limitations of large-offset observation data, insufficient illumination, and multi-parameter crosstalk effects can introduce significant uncertainties into the full-waveform inversion results. According to Eikrem, KS, G. Naevdal, and M. Jakobsen, 2019, Iterated Extended Kalman Filter Method for Time-lapse Seismic Full Waveform Inversion: Geophysical Prospecting, 67, no. 2, 379-394, doi:10.1111 / 1365-2478.12730, compared with traditional inversion methods, the Bayesian inversion method has the advantage of being able to quantify the uncertainty of the inversion results, thus enabling the effective use of prior information. The scattering integral method presents the forward modeling results in the form of the Lippmann-Schwinger equation, replacing the traditional finite difference forward modeling method. This avoids the grid difference errors that the traditional method may cause, and is more conducive to achieving fine structural imaging of subsurface media (Jakobsen, M., E. Ivan, I. Psencik, and B. Ursin, 2020, Transition operator approach to seismic full-waveform inversion in arbitrary anisotropic elastic media: Communications in Computational Physics, 27, no. 1, 1–31, doi:10.4208 / cicp.OA-2018-0197.). The Bayesian full-waveform inversion method based on scattering theory can more reliably identify subsurface media parameters and provide uncertainty information. It is suitable for solving multi-parameter inversion problems in anisotropic media and is of great significance for interpreting geological structures and updating reservoir models. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide a multi-parameter full waveform inversion method for anisotropic media that includes uncertainty analysis. This method can obtain high-resolution multi-parameter information of subsurface anisotropic media when seismic data lacks low-frequency information and has a low signal-to-noise ratio, and quantitatively estimate the uncertainty of the full waveform inversion results.

[0005] This invention is implemented as follows:

[0006] A multi-parameter full-waveform inversion method for anisotropic media incorporating uncertainty analysis, the method comprising:

[0007] Step 1: Input the observed seismic data, extract the travel time information of the seismic data, perform tomographic imaging, and obtain the tomographic velocity model;

[0008] Step 2: Using Fourier transform, the time-domain seismic data is converted to the frequency domain, its spectrum is plotted, and the frequency band with the richest effective information is determined to obtain the frequency-divided seismic data d. obs ;

[0009] Step 3: Use the Gardner formula to obtain the density model, combine it with the tomographic velocity model to establish the independent elastic modulus models in the anisotropic VTI medium, and use the Bond transformation to obtain the independent elastic modulus of the anisotropic TTI medium as the initial model for inversion.

[0010] Step 4: Based on scattering theory, solve the frequency domain elastic dynamics wave equation of the anisotropic TTI medium for forward numerical simulation, and construct the displacement-strain coupled integral equation in the form of the Lippmann-Swinger scattering integral equation using third-order and fourth-order Green's functions.

[0011] Step 5: Using the preprocessing generalized minimum residual method including discrete wavelet transform preprocessing operators, solve the strain scattering integral equation in the Krylov subspace, and then solve for the particle displacement field u to obtain the simulated seismic dataset d. cal ;

[0012] Step 6: Calculate the data residuals between observed seismic data and simulated seismic data. Based on the mapping relationship between the model and the data, obtain the sensitivity kernel according to the data residuals.

[0013] Step 7: Based on the frequency-division seismic data d obs Based on the Bayesian probability system, the maximum a posteriori solution and a posteriori covariance of the full waveform inversion of anisotropic TTI media are constructed.

[0014] Step 8: Use the iterative extended Kalman filter method to iteratively obtain the maximum a posteriori solution of the Bayesian full waveform inversion model parameters and the update amount of the a posteriori covariance. If the iteration stopping condition is met, output the model parameters and a posteriori covariance of this frequency as the input for the next frequency.

[0015] Step 9: Repeat steps 4-8 to iteratively update the next frequency model parameters and the posterior covariance matrix until the calculated frequency value reaches the set maximum frequency value and the final data residual meets the accuracy requirements. The model parameters obtained at this time are the final inversion results of the anisotropic TTI medium.

[0016] Further, in step 3, the density model is obtained using the Gardner formula, and an independent elastic modulus model for each anisotropic VTI medium is established by combining it with the tomographic velocity model. The independent elastic modulus of the anisotropic TTI medium is obtained using the Bond transformation, serving as the initial model for the inversion. Specifically, this includes:

[0017] Using the P-wave velocity v in the initial velocity model P and S-wave velocity v S The density ρ and the Thomsen parameters ε and δ of the Thomsen anisotropic parameter model characterize the elastic modulus of the anisotropic VTI medium, forming the VTI elastic modulus matrix:

[0018]

[0019] Assuming the angle between the z-axis and the axis of symmetry of the anisotropic medium is θ, the elastic modulus matrix of the anisotropic VTI medium is subjected to a Bond transformation relative to the axis of symmetry. A new set of stiffness tensors is obtained, which is the elastic modulus matrix of the anisotropic TTI medium represented by formula (2), and serves as the initial model for the full-waveform inversion of the independent elastic moduli in the anisotropic TTI medium:

[0020]

[0021] Further, in step 4, the frequency domain elastodynamic wave equation of the anisotropic TTI medium is solved for forward numerical simulation. A displacement-strain coupled integral equation in the form of the Lippmann-Swinger scattering integral equation is constructed using third- and fourth-order Green's functions. Specifically, this includes:

[0022] The particle displacement field is obtained by forward modeling using the frequency domain elastic dynamics wave equation of formula (3):

[0023]

[0024] In the formula, ρ and ω represent density and angular frequency, respectively, u(x) is the particle displacement, and S(x) is the source. Let I represent the strain field at point x, and let I denote the identity matrix; decompose the elastic constant C(x) into the background medium C. (0) Two components, δC(x) and the perturbation medium δC(x), are introduced, along with the background Green's function. Solving equation (3), and then utilizing the local integral symmetry of the elastic modulus and the displacement-strain relationship of the elastodynamic particle, we obtain the integral equations of the displacement-strain field as equations (4) and (5):

[0025]

[0026]

[0027] The third and fourth Green's function tensors can be written as follows:

[0028] Furthermore, the displacement-strain equation is written in the form of a pair of integral coupled operators according to equations (4) and (5):

[0029]

[0030]

[0031] Among them, among them, These are third-order and fourth-order Green's function tensors, u (0) ε (0) The background medium displacement and strain are given by V(x1,x2)=δC(x1)(x1-x2), which is the scattering potential operator, and ε is the strain field. x1 and x2 are two points in the model region that are close but do not coincide. δC(x1) represents the difference in elastic modulus at point x1 before and after the iteration update. The simulated seismic dataset d is obtained by solving formula (6). cal .

[0032] Furthermore,

[0033] The strain field of formula (7) is solved by an iterative solution method. Formula (7) is written as:

[0034]

[0035] make for And matrix ε (0) Divide into n column matrix blocks, and...

[0036] Solving Equation 8 is transformed into solving for the matching ε. (0) and The approximate solution x that minimizes the relationship between the two (m) Its mathematical expression is:

[0037]

[0038] Constructing a matrix using the Arnoldi loop orthogonal Krylov subspace

[0039]

[0040] For any Let x = x (0) +V m The formulas for calculating y and residual r are:

[0041]

[0042] Among them, V m ={v1,v2,...,v m}yes An orthonormal basis that satisfies AV m =V m+1 H m+1,m The least squares problem in formula (9) is simplified to:

[0043]

[0044] Where e1 refers to the unit vector and β refers to the quadratic norm of the residual.

[0045] Furthermore, when the iteration number k is reached, let x (0) =x (k) A new round of iterations will begin again until convergence, and then... As a preprocessing operator, it makes:

[0046]

[0047] Furthermore, step 6 specifically includes the following steps:

[0048] Step 6.1: Solve for the data residual of the i-th iteration. This is used to measure whether the current accuracy of the model meets the requirements for stopping the iteration, where ΔC (i+1) (x) represents the difference in elastic modulus between the current and previous iterations:

[0049]

[0050] Step 6.2: Introduce the calculation formula for elastic modulus disturbance. Equation (14) can be rewritten as:

[0051]

[0052] Where, m (p) =(C(x)-C (0) (x)) / C (0) (x) represents the model parameters, and p represents the number of elastic moduli; B (p) (x) is the structure matrix related to the model parameters, and p represents one of the 21 elastic tensors;

[0053] Step 6.3: Based on the relationship between the model and the data δu (i) =J (i) ·δm (i) The sensitivity kernel is calculated as follows:

[0054]

[0055] Step 6.4: Write the sensitivity kernel and perturbation update in the form of a block matrix:

[0056] J (i) =[J (i,1) J (i,2) ,...,J (i,21) ] T (17)

[0057] δm (i) =[δm (i,1) ,δm (i,2) ,...,δm (i,21) ] T (18).

[0058] Furthermore, step 7 specifically includes:

[0059] The maximum a posteriori solution and the a posteriori covariance of the model parameters based on the Bayesian probability system are calculated using formulas (19) and (20):

[0060]

[0061]

[0062] Among them, C D C M These are the data covariance and the model covariance, respectively, m pri As a priori model, d k Represents the frequency-division seismic dataset d obs The kth observation dataset, d cal For calculating data, m pri These are the prior model parameters.

[0063] Furthermore, step 8 specifically includes:

[0064] The iterative solution for the maximum a posteriori solution and the a posteriori covariance of the Bayesian probability system, obtained using the iterative extended Kalman filter method, is as follows:

[0065]

[0066]

[0067] Where λ is the dynamic factor and I is the identity matrix.

[0068] Compared with the prior art, the beneficial effects of this invention are as follows:

[0069] This invention is used to process conventional frequency domain seismic data. It introduces a Bayesian probabilistic framework and establishes a maximum a posteriori solution based on a probabilistic system, replacing the traditional full-waveform inversion method that uses least squares of observed and simulated data as the objective function. Simultaneously, it inverts the independent elastic parameters of anisotropic TTI media, obtaining uncertainty information along with the inversion results. The model and covariance information update strategy adopts a multi-scale inversion strategy of "low frequency to high frequency, a posteriori is prior". Based on scattering theory, it fully considers the multi-scattering properties of the subsurface medium to achieve fine imaging of subsurface structures. The scattering integral method is used for forward numerical simulation of the model. By using the Green's function to calculate the multi-scattering sensitivity kernel, approximate Hessian information is obtained, and then the covariance matrix is ​​constructed. When the background model complexity is low, the scattering integral equation automatically satisfies the Sommerfeld radiation condition, eliminating the need to set absorption boundaries and reducing the uncertainty of the inversion results caused by grid discretization errors. The scattering integral method allows the inversion to discretize only the anomaly region. Based on this, combined with a linear equation iterative solution operator, the computational cost of large-scale model inversion results is reduced to some extent. This invention improves the quality of inversion results while increasing computational efficiency and enhancing its effectiveness for industrial applications. Even in the absence of accurate initial models and low-frequency seismic records, the inversion results exhibit good robustness, making it suitable for research and industrial production applications in various fields, including terrestrial and marine environments. Attached Figure Description

[0070] Figure 1 This is a flowchart of a traditional full-waveform inversion algorithm;

[0071] Figure 2 This is a flowchart of a multi-parameter full waveform inversion method for anisotropic media that incorporates uncertainty analysis, as proposed in this invention.

[0072] Figure 3 The model provided in this invention is a real parameter model for simulation testing, where (a)-(f) show the P-wave and S-wave velocities, densities, Thomsen parameters, and model anisotropic tilt angles, respectively.

[0073] Figure 4 This is the true elastic modulus model for simulation testing provided by the present invention, where (a)-(e) show different elastic moduli C. 11 C 33 C 55 C 66 C 13 ;

[0074] Figure 5 This is the initial elastic modulus model for simulation testing provided by the present invention, where (a)-(e) show different elastic moduli C. 11 C 33 C 55 C 66 C 13 ;

[0075] Figure 6 This is the final inversion result of the elastic modulus model in the simulation test provided by the present invention, where (a)-(e) show different elastic moduli C. 11 C 33 C 55 C 66 C 13 ;

[0076] Figure 7 The error between the elastic modulus model inversion results and the true model in the simulation test provided by this invention is shown in (a)-(e), where different elastic moduli C are represented. 11 C 33 C 55 C 66 C 13 ;

[0077] Figure 8 This is the posterior standard deviation of the elastic modulus model inversion results provided in this invention, where (a)-(e) show different elastic moduli C. 11 C 33 C 55 C 66 C 13 ;

[0078] Figure 9 The prior covariance matrix of any elastic modulus in the simulation test provided by this invention;

[0079] Figure 10 The prior covariance matrix of all elastic moduli in the simulation test provided by this invention is represented by the block matrices on the upper left to lower right diagonal, respectively: C 11 C 33 C 55 C66 C 13 ;

[0080] Figure 11 The posterior covariance matrix of all elastic moduli in the simulation test provided by this invention is represented by the block matrices on the top left to bottom right diagonal, respectively: C 11 C 33 C 55 C 66 C 13 . Detailed Implementation

[0081] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0082] For traditional full waveform inversion methods, see Figure 1 As shown, the flowchart of the multi-parameter full waveform inversion method for anisotropic media incorporating uncertainty analysis proposed in this invention is as follows: Figure 2 As shown, the implementation process, using a theoretical test example, includes the following steps:

[0083] Step 1: Input the observed seismic data, extract the travel time information of the seismic data, perform tomographic imaging, and obtain the tomographic velocity model;

[0084] Step 2: Using Fourier transform, the time-domain seismic data is converted to the frequency domain, its spectrum is plotted, and the frequency band with the richest effective information is determined to obtain the frequency-divided seismic data d. obs ;

[0085] Step 3: Use the Gardner formula to obtain the density model, combine it with the tomographic velocity model to establish the independent elastic modulus models in the anisotropic VTI medium, and use the Bond transformation to obtain the independent elastic modulus of the anisotropic TTI medium as the initial model for inversion.

[0086] Step 4: Based on scattering theory, solve the frequency domain elastic dynamics wave equation of the anisotropic TTI medium for forward numerical simulation, and construct the displacement-strain coupled integral equation in the form of the Lippmann-Swinger scattering integral equation using third-order and fourth-order Green's functions.

[0087] Step 5: Using the preprocessing generalized minimum residual method including discrete wavelet transform preprocessing operators, solve the strain scattering integral equation in the Krylov subspace, and then solve for the particle displacement field u to obtain the simulated seismic dataset d. cal ;

[0088] Step 6: Calculate the data residuals between observed seismic data and simulated seismic data. Based on the mapping relationship between the model and the data, obtain the sensitivity kernel according to the data residuals.

[0089] Step 7: Based on the frequency-division seismic data d obs Based on the Bayesian probability system, the maximum a posteriori solution and a posteriori covariance of the full waveform inversion of anisotropic TTI media are constructed.

[0090] Step 8: Use the iterative extended Kalman filter method to iteratively obtain the maximum a posteriori solution of the Bayesian full waveform inversion model parameters and the update amount of the a posteriori covariance. If the iteration stopping condition is met, output the model parameters and a posteriori covariance of this frequency as the input for the next frequency.

[0091] Step 9: Repeat steps 4-8 to iteratively update the next frequency model parameters and the posterior covariance matrix until the calculated frequency value reaches the set maximum frequency value and the final data residual meets the accuracy requirements. The model parameters obtained at this time are the final anisotropic TTI medium inversion results.

[0092] Step 3 specifically includes the following steps:

[0093] Step 1.1: Use P-wave velocity v P S-wave velocity v S Density ρ, and the Thomsen parameters ε and δ characterize the elastic modulus of anisotropic VTI media:

[0094]

[0095] Step 1.2: The elastic modulus matrix of the anisotropic VTI medium is constructed as follows:

[0096]

[0097] Step 1.3: Assuming the angle between the z-axis and the axis of symmetry of the anisotropic medium is θ, perform a Bond transformation on the elastic modulus matrix of the anisotropic VTI medium relative to the axis of symmetry to obtain a new set of stiffness tensors, which is the elastic modulus matrix of the anisotropic TTI medium:

[0098]

[0099] The formula for calculating the tilt angle matrix is:

[0100]

[0101] Step 1.5: Obtain the independent elastic moduli in the stiffness matrix of the anisotropic TTI medium based on the elastic modulus matrix of the anisotropic TTI medium, and use them as parameters for the anisotropic model in the full waveform inversion:

[0102]

[0103] In step 4, a forward numerical simulation is performed on the initial model of the anisotropic TTI medium based on scattering theory. A displacement-strain coupled integral equation in the form of the Lippmann-Swinger scattering integral equation is constructed. Then, the generalized minimum residual method based on Krylov subspace is used to accelerate the solution of the strain of the anisotropic TTI medium, thereby solving for the particle displacement field and obtaining the simulated seismic dataset. This includes the following steps:

[0104] Introducing the background Green's function Solving the frequency domain elastic dynamics wave equation, we obtain the particle displacement field expression:

[0105]

[0106] In the formula, x′ is a point near x, and u (0) (x) represents the wave field of the background medium;

[0107] Taking advantage of the local integral symmetry of the elastic modulus, equation (4) can be rewritten as:

[0108]

[0109] in It is a third-order Green's function tensor;

[0110] From the particle displacement-strain relationship in elastic dynamics, the integral equation of the strain field is obtained:

[0111]

[0112] in It is a fourth-order Green's function tensor.

[0113] Based on formulas (4) and (5), the displacement-strain equation can be written in the form of a pair of integral coupled operators:

[0114]

[0115]

[0116] Where V(x1,x2)=δC(x1)(x1-x2) is the scattering potential operator;

[0117] Introducing iterative solution methods to solve strain fields in linear systems:

[0118]

[0119] make for And matrix ε(0) Dividing the equation into n column matrix blocks, the solution for the dependent variable is rewritten as:

[0120]

[0121] Solving for the matching ε (0) and The approximate solution x that minimizes the relationship between the two (m) for:

[0122]

[0123] Constructing a matrix using the Arnoldi loop orthogonal Krylov subspace κ m :

[0124]

[0125] For any x = x (0) +κ m Let x = x (0) +V m The formula for calculating y and residual r is formula (11):

[0126]

[0127] Among them, V m ={v1,v2,...,v m}yes An orthonormal basis that satisfies AV m =V m+1 H m+1,m ;

[0128] Step 2.11: The least squares problem is simplified to:

[0129]

[0130] Here, e1 refers to the unit vector, and β refers to the quadratic norm of the residual.

[0131] Make As a preprocessing operator, when the restart iteration step k is reached, let x (0) =x (k) Restart a new round of iterations until convergence:

[0132]

[0133] Step 6 specifically includes the following steps:

[0134] Step 6.1: Solve for the data residual of the i-th iteration. This is used to measure whether the current accuracy of the model meets the requirements for stopping the iteration, where ΔC(i+1) (x) represents the difference in elastic modulus between the current and previous iterations:

[0135]

[0136] Step 6.2: Introduce the calculation formula for elastic modulus disturbance. Equation (14) can be rewritten as:

[0137]

[0138] Where, m (p) =C (0) (x)(I+ΔC(x)) are model parameters, B (p) (x) is the structure matrix related to the model parameters;

[0139] Step 6.3: Based on the relationship between the model and the data, the sensitivity kernel is obtained as follows:

[0140]

[0141] Step 6.4: Write the sensitivity kernel and perturbation update in the form of a block matrix:

[0142]

[0143] The data residuals are expressed as the product of the model perturbation and the sensitivity kernel:

[0144] δu (i) =J (i) ·δm (i)

[0145] Step 7 specifically includes the following steps:

[0146] Step 7.1: Based on the Bayesian probability system, the expressions for the maximum a posteriori solution and the a posteriori covariance of the full waveform inversion model parameters are obtained as follows:

[0147]

[0148]

[0149] Among them, C D C M These are the data covariance and the model covariance, respectively, m pri As a priori model, d k Represents the frequency-division seismic dataset d obs The kth observation dataset, d cal For calculating data, m pri These are the prior model parameters.

[0150] The model parameters are iteratively solved using the iterative extended Kalman filter method, yielding the maximum a posteriori solution and the iterative solution of the a posteriori covariance:

[0151]

[0152]

[0153] Example

[0154] Establish an initial parameter model for the inversion of the elastic modulus of anisotropic TTI media.

[0155] according to Figure 3 Based on the actual parameter information shown, and setting the tilt angle θ, a true elastic modulus model of the anisotropic TTI medium is constructed, as follows: Figure 4 As shown, the calculation formula is as follows:

[0156] in

[0157]

[0158] Based on existing seismic profile results, well logging data, geological and lithological data, an initial parameter model for the elastic modulus of the tested anisotropic TTI medium is established. In this example, the initial model is similar to a smooth gradient model, such as... Figure 5 As shown.

[0159] Set up an appropriate seismic observation system based on the actual parameter model size and add noise interference to make the experiment more closely resemble the actual situation.

[0160] The model mesh was set to 40×20 with a grid spacing of 25 meters, and the model size was 1.0×0.5km. A 7.5Hz Ricker wavelet was used as the source wavelet, with 40 sources and 40 detectors, all uniformly placed at the top of the model. The background medium model was a Gaussian-smoothed homogeneous isotropic medium model.

[0161] The signal-to-noise ratio is set to 4. The formula for calculating the signal-to-noise ratio is: Where d is the frequency data, w is the angular frequency, and E is the expected value.

[0162] Set the model inversion frequency range to 4Hz to 20Hz.

[0163] The inversion process proceeds sequentially from low to high frequencies. The model obtained at the previous frequency will serve as the initial model for the next frequency update. Following the set frequency range from low to high, the inversion results are as follows: Figure 6 As shown. The error between the maximum a posteriori solution and the true model is as follows. Figure 7 As shown, the standard deviation is as follows Figure 8As shown. The prior covariance matrices of any elastic modulus and all elastic moduli in the model are respectively as follows: Figure 9 and Figure 10 As shown, the posterior covariance matrix of the inversion result is as follows: Figure 11 As shown.

[0164] The inversion results show that even with a significant difference between the initial velocity model and the actual model, and a low signal-to-noise ratio, this invention can still accurately invert the five independent elastic moduli of anisotropic TTI media. Although some errors still exist in the deeper parts of the model and at the boundaries, the inversion results represent a significant improvement over the initial model, achieving a relatively ideal effect. Furthermore, the crosstalk relationship between the various parameters can be observed from the posterior covariance and the resolution of the inversion results. Parameter C 33 and C 55 The inversion results have the highest resolution, followed by C. 11 and C 66 Finally, it's C. 13 Because C 33 and C 55 These coefficients can be considered to be related only to the vertical velocities of P-waves and S-waves; therefore, they are more independent than other coefficients during the update process, mainly because surface data consists of seismic records with vertical components. Of the five elastic constants, only C... 66 Related to the SH wave, the other four elastic parameters C 11 C 55 C 33 C 13 It is related to the phase velocity and polarization of P-waves and SV-waves. However, C 11 and C 33 There is also some crosstalk, and to a certain extent, the inversion results of the two will influence each other. This is because C 11 and C 33 It approximates the so-called P-wave anisotropy, which is related to the velocity difference of the P-wave in the horizontal and vertical directions. Furthermore, the anisotropy intensity of the SH wave is relatively weaker compared to the other four elastic moduli. This is because C... 13 Not only will it be affected by C 55 It will also be affected by C 33 The influence of these three coefficients. These coefficients are closely related to the second derivative of the P-wave phase velocity at vertical incidence.

[0165] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A multi-parameter full-waveform inversion method for anisotropic media incorporating uncertainty analysis, characterized in that, The method includes: Step 1: Input the observed seismic data, extract the travel time information of the seismic data, perform tomographic imaging, and obtain the tomographic velocity model; Step 2: Using Fourier transform, the time-domain seismic data is converted to the frequency domain, its spectrum is plotted, the frequency band range of effective information is determined, and frequency-divided seismic data is obtained. ; Step 3: Use the Gardner formula to obtain the density model, combine it with the tomographic velocity model to establish the independent elastic modulus models in the anisotropic VTI medium, and use the Bond transformation to obtain the independent elastic modulus of the anisotropic TTI medium as the initial model for inversion. Step 4: Based on scattering theory, solve the frequency domain elastic dynamics wave equation of the anisotropic TTI medium for forward numerical simulation, and use third- and fourth-order Green's functions to construct the displacement-strain coupled integral equation in the form of the Lippmann-Swinger scattering integral equation. Step 5: Using the preprocessing generalized minimum residual method including discrete wavelet transform preprocessing operators, solve the strain scattering integral equation in the Krylov subspace, and then solve for the particle displacement field. To obtain a simulated earthquake dataset ; Step 6: Calculate the data residuals between observed seismic data and simulated seismic data. Based on the mapping relationship between the model and the data, obtain the sensitivity kernel according to the data residuals. Step 7: Based on frequency-division seismic data Based on the Bayesian probability system, the maximum a posteriori solution and a posteriori covariance of the full waveform inversion of anisotropic TTI media are constructed. Step 8: Use the iterative extended Kalman filter method to iteratively obtain the maximum posterior solution of the Bayesian full waveform inversion model parameters and the posterior covariance update. If the iteration stopping condition is met, output the model parameters and posterior covariance of this frequency as the input for the next frequency. Step 9: Repeat steps 4-8 to iteratively update the next frequency model parameters and posterior covariance matrix information until the calculated frequency value reaches the set maximum frequency value and the final data residual meets the accuracy requirements. The model parameters obtained at this time are the final anisotropic TTI medium inversion results.

2. The method for multi-parameter full waveform inversion of anisotropic media including uncertainty analysis according to claim 1, characterized in that, In step 3, the density model is obtained using the Gardner formula, and an independent elastic modulus model for the anisotropic VTI medium is established by combining it with the tomographic velocity model. The independent elastic modulus of the anisotropic TTI medium is obtained using the Bond transformation, which serves as the initial model for the inversion. Specifically, this includes: P-wave velocity from the initial velocity model and S-wave velocity ,density and the Thomsen parameters of the Thomsen anisotropic parameter model and Characterizing the elastic modulus of anisotropic VTI media, the VTI elastic modulus matrix is ​​constructed as follows: （1） Assume the angle between the z-axis and the symmetry axis of the anisotropic medium is... The elastic modulus matrix of the anisotropic VTI medium is transformed by Bond transformation relative to the axis of symmetry. A new set of stiffness tensors is obtained, which is the elastic modulus matrix of the anisotropic TTI medium expressed by formula (2), and serves as the initial model for the full-waveform inversion of the independent elastic moduli in the anisotropic TTI medium: （2）。 3. The method for multi-parameter full waveform inversion of anisotropic media including uncertainty analysis according to claim 1, characterized in that, In step 4, the frequency domain elastic dynamics wave equation of the anisotropic TTI medium is solved for forward numerical simulation. A displacement-strain coupled integral equation in the form of the Lippmann-Swinger scattering integral equation is constructed using third- and fourth-order Green's functions. Specifically, this includes: The particle displacement field is obtained by forward modeling using the frequency domain elastic dynamics wave equation of formula (3): （3） In the formula, and Representing density and angular frequency, respectively. This is the particle displacement. As the epicenter, For point The strain field at point I, where I represents the identity matrix; the elastic constants... Decomposed into background medium and disturbance medium Two components, introducing the background Green's function Solving equation (3), and then using the local integral symmetry of the elastic modulus and the displacement-strain relationship of the elastodynamic particle, the integral equations of the displacement-strain field are obtained as equations (4) and (5): （4） （5） The third and fourth Green's function tensors are written as follows: , .

4. The method for multi-parameter full waveform inversion of anisotropic media including uncertainty analysis according to claim 3, characterized in that, The displacement-strain equation can be written in the form of a pair of integral coupled operators according to formulas (4) and (5): （6） （7） Among them, among them, , These are third-order and fourth-order Green's function tensors, respectively. , For the background medium displacement and strain, For scattering potential operators, For strain field, and These are two points in the model region that are close but do not overlap. express The difference in elastic modulus before and after the iterative update at a point is used to obtain the simulated earthquake dataset by solving formula (6). .

5. The method for multi-parameter full waveform inversion of anisotropic media including uncertainty analysis according to claim 4, characterized in that, The strain field of formula (7) is solved by an iterative solution method. Formula (7) is written as: （8） make for and the matrix Divide into n column matrix blocks, and... ; Solving formula (8) is transformed into solving a matching problem. and Approximate solution that minimizes the relationship between Its mathematical expression is: （9） Constructing a matrix using the Arnoldi loop orthogonal Krylov subspace : （10） For any ,set up residual The calculation formula is: （11） in, yes An orthonormal basis that satisfies The least squares problem in formula (9) is simplified to: （12） in, Unit vector The quadratic norm of the residual.

6. The method for multi-parameter full waveform inversion of anisotropic media including uncertainty analysis according to claim 5, characterized in that, When the iteration step k is reached, let A new round of iterations will begin again until convergence, and then... As a preprocessing operator, it makes: （13）。 7. The method for multi-parameter full waveform inversion of anisotropic media including uncertainty analysis according to claim 1, characterized in that, Step 6 specifically includes the following steps: Step 6.1: Solve for the data residual of the i-th iteration. This is used to measure whether the current accuracy of the model meets the requirements for stopping the iteration. The difference between the elastic modulus of the current iteration and the previous iteration: （14） Step 6.2: Introduce the calculation formula for elastic modulus disturbance. Formula (14) can be rewritten as: （15） in, For model parameters, The number of elastic moduli; The structure matrix is ​​related to the model parameters, where p represents one of the 21 elasticity tensors. The structure matrix is ​​related to the model parameters; Step 6.3: Based on the relationship between the model and the data The sensitivity kernel is calculated as follows: （16） Step 6.4: Write the sensitivity kernel and perturbation update in the form of a block matrix: （17） （18）。 8. The method for multi-parameter full waveform inversion of anisotropic media including uncertainty analysis according to claim 7, characterized in that, Step 7 specifically includes: The maximum a posteriori solution and the a posteriori covariance of the model parameters based on the Bayesian probability system are calculated using formulas (19) and (20): （19） （20） in, , These are the data covariance and the model covariance, respectively. As a priori model, Represents a frequency-division seismic dataset The Group observation dataset, To calculate the data, These are the prior model parameters.

9. The method for multi-parameter full waveform inversion of anisotropic media including uncertainty analysis according to claim 8, characterized in that, Step 8 specifically includes: The iterative solution for the maximum a posteriori solution and the a posteriori covariance of the Bayesian probability system model is obtained using the iterative extended Kalman filter method: （21） （22） in I is the dynamic factor, and I is the identity matrix.

Images