Young impedance and longitudinal and transverse wave velocity ratio nonlinear direct seismic inversion method based on precise YIRD-Zoeppritz equation

By using the precise YIRD-Zoeppritz equation and the ADA-MCMC-LpSC nonlinear inversion algorithm, Young's impedance and P-wave velocity ratio are directly inverted from seismic data, solving the problem that existing technologies cannot directly predict these parameters and achieving higher accuracy inversion results.

CN121679697APending Publication Date: 2026-03-17CHENGDU UNIVERSITY OF TECHNOLOGY
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2025-12-03
Publication Date
2026-03-17

AI Technical Summary

Technical Problem

Existing technologies cannot directly use the exact Zoeppritz equation to simultaneously predict Young's impedance and the P-wave/S-wave velocity ratio, and traditional methods suffer from accumulated errors and the inability to achieve direct inversion simultaneously due to linear approximation.

Method used

Based on the precise YIRD-Zoeppritz equation and combined with the ADA-MCMC-LpSC nonlinear inversion algorithm that provides sparse prior information using the Lp norm, Young's impedance and P-wave velocity ratio are directly nonlinearly inverted from pre-stack seismic data with multiple incident angles.

Benefits of technology

It provides a more accurate equation basis, realizes the direct and accurate inversion of Young's impedance and P-wave velocity ratio, solves key geophysical problems that cannot be directly predicted in existing technologies, and improves the accuracy and precision of the inversion results.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN121679697A_ABST
    Figure CN121679697A_ABST
Patent Text Reader

Abstract

The invention discloses a Young impedance and longitudinal and transverse wave velocity ratio nonlinear direct seismic inversion method based on an accurate YIRD-Zoeppritz equation. The Young impedance and longitudinal and transverse wave velocity ratio nonlinear direct seismic inversion method comprises the steps that the YIRD-Zoeppritz equation is obtained, and the Young impedance YI, the longitudinal and transverse wave velocity ratio PSR and the rock density are inversed through an ADA-MCMC-LpSC nonlinear inversion algorithm providing sparse prior information based on the Lp norm. Based on an accurate YIRD-Zoeppritz equation, the accurate prediction method for the Young impedance and the longitudinal and transverse wave velocity ratio under the ADA-MCMC-LpSC nonlinear optimization inversion algorithm is achieved, and the key geophysical problem that the accurate Zoeppritz equation cannot be directly used for predicting the Young impedance and the longitudinal and transverse wave velocity ratio parameters at the same time at present is solved.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of oil and gas field development engineering technology, and is based on a nonlinear direct seismic inversion method for Young's impedance and P-wave velocity ratio using the accurate YIRD-Zoeppritz equation. Background Technology

[0002] For predicting Young's impedance and P-wave velocity ratio of unconventional fractured gas-bearing shale reservoirs, there are two existing traditional methods: one is to indirectly calculate Young's impedance and P-wave velocity ratio using the inversion results of P-wave velocity, S-wave velocity, and rock density; the other is to derive a linear approximate equation containing Young's impedance or P-wave velocity ratio under the assumption of weak correlation of elastic parameters at the formation interface for direct inversion. Its limitations are mainly reflected in: (1) It requires the inversion results of P-wave velocity, S-wave velocity, and rock density to indirectly calculate Young's impedance and P-wave velocity ratio. This calculation process cannot avoid the accumulation of errors, resulting in inaccurate prediction of reservoir characterization parameters; (2) The linear approximation proposed in existing studies cannot simultaneously achieve direct inversion of Young's impedance and P-wave velocity ratio. Summary of the Invention

[0003] The purpose of this invention is to provide a nonlinear direct seismic inversion method based on the exact YIRD-Zoeppritz equation for Young's impedance and P-wave velocity ratio, in order to solve the key geophysical problem that it is currently impossible to directly use the exact Zoeppritz equation to simultaneously predict Young's impedance and P-wave velocity ratio parameters.

[0004] The embodiments of this application are implemented as follows: a nonlinear direct seismic inversion method based on Young's impedance and P-wave / S-wave velocity ratio using the exact YIRD-Zoeppritz equations, including: The YIRD-Zoeppritz equations are obtained as follows: In the formula, intermediate parameters , , , Meanwhile, here , , and Existing parameters and ; PSR Indicates the ratio of longitudinal to transverse wave velocity; R PP , T PP , R PS and T PSThese represent the P-wave reflection coefficient, P-wave transmission coefficient, converted S-wave reflection coefficient, and converted S-wave transmission coefficient, respectively; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. i 1 represents the angle of incidence of the reflected longitudinal wave; V p , V s and r These represent the longitudinal wave velocity, the transverse wave velocity, and the rock density, respectively. YI Indicates Young's impedance; The ADA-MCMC-LpSC nonlinear inversion algorithm, which uses sparse prior information provided by the Lp norm, inverts Young's impedance. YI S-wave velocity ratio PSR and rock density .

[0005] Optionally, in some embodiments of this application, the method for obtaining the YIRD-Zoeppritz equation includes: Based on the relationship between elastic parameters, both the longitudinal wave modulus and the transverse wave modulus can be expressed using Young's modulus and Poisson's ratio, as shown below: ; ; In the formula, M Represents the longitudinal wave modulus; m Indicates shear modulus; s Indicates Poisson's ratio; E This indicates Young's modulus; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. The relationship between Poisson's ratio and the ratio of P-wave velocity to S-wave velocity is nonlinear, as shown below: ; In the formula, s Indicates Poisson's ratio; PSR This represents the ratio of longitudinal to transverse wave speeds, i.e. PSR=V p / V s , V p and V s These represent the P-wave velocity and S-wave velocity, respectively; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. Young's impedance is defined as the product of Young's modulus and rock density, and its expression is as follows: ; In the formula, the parameters YI Indicates Young's impedance; E Indicates Young's modulus; rThis indicates the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. Formula and Substitute into the formula respectively and The following equation can be obtained: ; ; In the formula, M Represents the longitudinal wave modulus; m Describe the shear modulus; define two parameters related to... PSR Intermediate parameters: and ; YI Indicates Young's impedance; r This indicates the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. Longitudinal wave velocity and transverse wave velocity can be expressed in the following form: ; ; In the formula, V p , V s and r These represent P-wave velocity, S-wave velocity, and rock density, respectively; intermediate parameters and ; PSR Indicates the ratio of longitudinal to transverse wave velocity; YI Indicates Young's impedance; r This indicates the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. Based on the accurate Zoeppritz equations for the continuity of displacement and stress at the interface, which include the ratio of longitudinal to transverse wave velocities. PSR and Young's impedance YI The new exact equation is expressed as follows, and this new equation is called the YIRD-Zoeppritz equation: In the formula, intermediate parameters , , , Meanwhile, here , , and Existing parameters and ; PSR Indicates the ratio of longitudinal to transverse wave velocity; R PP ,T PP , R PS and T PS These represent the P-wave reflection coefficient, P-wave transmission coefficient, converted S-wave reflection coefficient, and converted S-wave transmission coefficient, respectively; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. i 1 represents the angle of incidence of the reflected longitudinal wave; V p , V s and r These represent the longitudinal wave velocity, the transverse wave velocity, and the rock density, respectively. YI This represents Young's impedance.

[0006] Optionally, in some embodiments of this application, the exact Zoeppritz equation based on the continuity of displacement and stress at the interface is as follows: ; In the formula, R PP , T PP , R PS and T PS These represent the longitudinal wave reflection coefficient, longitudinal wave transmission coefficient, converted transverse wave reflection coefficient, and converted transverse wave transmission coefficient, respectively; intermediate parameters. , , and The subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. i 1 represents the angle of incidence of the reflected longitudinal wave; V p , V s and r These represent the longitudinal wave velocity, the transverse wave velocity, and the rock density, respectively.

[0007] Optionally, some embodiments of this application include: The ADA-MCMC-LpSC nonlinear inversion algorithm, which provides sparse prior information based on the Lp norm, includes the following steps: Provide a likelihood function describing the probability density of seismic data: ; In the formula, (d|m) denotes the likelihood function; ||·||2 denotes the Euclidean 2 norm; the symbol exp denotes an exponential function with the natural constant e as the base; m The parameters to be inverted are represented by Young's impedance, the ratio of longitudinal to transverse wave velocity, and rock density.d Represents the seismic data matrix; R PP Represents the longitudinal wave reflection coefficient; W Represents the wavelet matrix; N For time sampling points; Here is the noise covariance matrix; symbol This is the value of the determinant, used for Gaussian distribution normalization; k The number of incident angles; To improve the accuracy of the inversion results, L is introduced into the prior probability density function. p The norm regularization term, with its prior distribution, is represented as: ; In the formula, Indicates having L The prior probability density function of the p-norm sparse constraint; ||·|2 represents the Euclidean 2-norm; The prior mean of the model is usually taken from the filter well curve; ∑ m Represents the prior covariance matrix with spatial correlation; symbol is the value of the determinant; the symbol exp represents an exponential function with the natural constant e as its base; m This represents the parameters to be inverted, consisting of Young's impedance, the ratio of longitudinal to transverse wave velocity, and the rock density. N The time sampling point is denoted by exp; the symbol exp represents an exponential function with the natural constant e as its base. R PP Represents the longitudinal wave reflection coefficient; ||·|| p Indicates the use of sparsity enhancement L p-norm; l The sparsity regularization parameter; 0 < p <1 is used to control the degree of sparsity; Based on the Bayesian framework, combining the likelihood function and the sparse-constrained prior distribution, the posterior probability distribution is formally represented as follows: ; In the formula, The proposed means having L The posterior probability distribution of the p-norm sparse constraint, where ||·|2 represents the Euclidean 2-norm; ||·| p Indicates the use of sparsity enhancement L p-norm; m This represents the parameters to be inverted, consisting of Young's impedance, the ratio of longitudinal to transverse wave velocity, and the rock density. d Represents the seismic data matrix; W The symbol exp represents the wavelet matrix; the symbol exp represents an exponential function with the natural constant e as the base. The prior mean of the model is usually taken from the filter well curve; lFor sparsity regularization parameters; Here is the noise covariance matrix; R PP Represents the longitudinal wave reflection coefficient; ||·|| p Indicates the use of sparsity enhancement L p-norm; 0 < p <1 is used to control the degree of sparsity; the symbol exp represents an exponential function with the natural constant e as the base; To further improve computational efficiency and accurately estimate the posterior probability distribution, a two-stage acceptance mechanism for the delayed acceptance Markov chain Monte Carlo algorithm is introduced. Finally, the Young's impedance YI, PSR, and rock density were obtained. The inversion value of the parameter.

[0008] Optionally, in some embodiments of this application, for a given angle of incidence... i k Its synthetic seismic data is represented as follows: ; In the formula, d This represents forward modeling seismic data after noise has been added. Indicates noise; subscript N Indicates the impedance of Young's method YI and PSR The number of sampling points with equal parameters; i k This represents the k-th angle of incidence; the subscript k indicates the index of the angle of incidence. W Indicates the relationship with the given first k An incident angle i k Related wavelet matrices; R PP Represents the longitudinal wave reflection coefficient; in W The definition is as follows: ; In the formula, W Indicates the relationship with the given first k An incident angle i k Related wavelet matrices; w Indicates the angle of incidence i k Seismic wavelet; subscript th Indicates the length of the seismic wavelet; matrix subscript N × N This indicates that the matrix contains N rows and N columns of elements.

[0009] Optionally, in some embodiments of this application, the current state is set to... And by candidate models To generate a new state, assuming a symmetric proposal distribution, the acceptance probability in the first stage is expressed as: ; In the formula, This represents the acceptance probability in the first stage. Represents the model parameter vector The prior probability distribution; This is the model parameter vector for the current (i-1)th iteration state; The candidate model parameters generated for the proposal distribution; symbol min{1, } indicates taking a value less than 1; d Represents the seismic data matrix; Let be the likelihood function, representing the candidate model parameter vector. The degree of matching between the seismic data obtained from forward modeling and the observed seismic data; Indicates the current model parameters The degree of matching between the seismic data obtained from forward modeling and the observed seismic data; For the current model Generate candidate models The distribution of proposals; To select candidate models Generate the current model The distribution of proposals, in the above formula = ; This represents the prior probability density of the candidate model; This represents the prior probability density of the current model.

[0010] Optionally, in some embodiments of this application, from a uniform distribution Generate a random number ,if If the proposal is rejected, the Markov chain remains in its current state. ; if Then candidate models are needed. With probability The second phase of evaluation involves the final acceptance of a sample based on its... L The second-stage acceptance rate of the exact target probability distribution with p-norm sparse constraints The calculation method is as follows: ; From uniform distribution Generate a random number ,if If the proposal is rejected, the Markov chain remains in its current state. ;if If so, then accept the proposal, that is... ; This indicates the effect of introducing the Lp regularization term on the candidate model. Posterior probability density; For the current model The posterior probability density; For the current iterative inversion i The model parameter vector of the -1st order state; Candidate model parameters generated for the proposal distribution; d This represents the earthquake data matrix.

[0011] Optionally, in some embodiments of this application, an adaptive proposal distribution consisting of two parts, proposed by Roberts and Rosenthal in 2009, is introduced to obtain candidate models. Its expression is as follows: ; In the formula, Represents a multivariate normal distribution; Represents a unit diagonal matrix, used in ∑ m Even with a poor condition number, the existence of proposed samples can still be guaranteed; N For time sampling points; ∑ m Represents the prior covariance matrix with spatial correlation; This represents the candidate model parameter vector consisting of Young's impedance, the ratio of longitudinal to transverse wave velocity, and the rock density. For the current iterative inversion i The model parameter matrix for the -1st order state; parameters β The mixed weight coefficients are in the range [0,1], and their update method is as follows: ; In the formula, the parameters β The mixed weighting coefficients are in the range [0,1]. Q This indicates the acceptance rate of the MCMC Markov chain.

[0012] Optionally, in some embodiments of this application, after the above-described multiple iterative inversion process and evaluation and acceptance by the two-stage acceptance mechanism, the model parameters... After the iterative burning period, the Markov chain finally reaches its final state on the second iteration. i When the iteration reaches t0, it converges to a stationary phase. Finally, the mean of all model parameter values ​​retrieved during the iterative inversion from the stationary phase t0 to the total number of iterations is calculated to obtain the impedance from Young's equation. YI S-wave velocity ratio PSR and rock density Inversion parameters composed of parameters m .

[0013] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are: This application derives a new equation, the YIRD-Zoeppritz equation, that simultaneously incorporates Young's impedance, P-wave velocity ratio, and rock density, based on the exact YIRD-Zoeppritz equation. This provides a more accurate equational foundation for the direct and accurate inversion of Young's impedance and P-wave velocity ratio in shale reservoirs. Based on a Bayesian inversion framework, an ADA-MCMC-LpSC nonlinear inversion algorithm, which provides sparse prior information based on the Lp norm, is proposed. This algorithm can directly and nonlinearly invert Young's impedance and P-wave velocity ratio parameters from pre-stack seismic data with multiple incident angles. This application implements an accurate prediction method for Young's impedance and P-wave velocity ratio using the ADA-MCMC-LpSC nonlinear optimized inversion algorithm based on the exact YIRD-Zoeppritz equation, solving the key geophysical problem that currently cannot directly predict Young's impedance and P-wave velocity ratio parameters simultaneously using the exact Zoeppritz equation. Attached Figure Description

[0014] Figure 1 A comparison chart of the accuracy of different equations and the exact Zoeppritz equation in Model I provided in the embodiments of the present invention; Figure 2 A comparison chart of the accuracy of different equations and the exact Zoeppritz equation in Model II provided in the embodiments of the present invention; Figure 3 A comparison chart of the accuracy of different equations and the exact Zoeppritz equation in Model III provided in the embodiments of the present invention; Figure 4 A comparison chart of the accuracy of different equations and the exact Zoeppritz equation in Model IV provided in the embodiments of the present invention; Figure 5 The noiseless inversion results of Young's impedance and P-wave velocity ratio based on the ADA-MCMC-LpSC algorithm provided in this embodiment of the invention are shown in the figure. Figure 6 The inversion result diagram of Young's impedance and signal-to-noise ratio of P-wave velocity ratio of 5 based on the ADA-MCMC-LpSC algorithm provided in this embodiment of the invention; Figure 7 This is a superimposed seismic data at different incident angles along a survey line in an actual work area, provided in an embodiment of the present invention.

[0015] Figure 8 The inversion results of the longitudinal and transverse wave velocity ratio, Young's impedance, and rock density corresponding to a certain measuring line in an actual work area provided in this embodiment of the invention. Detailed Implementation

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

[0017] The technical solution of this application is as follows: This application provides a nonlinear direct seismic inversion method based on the exact YIRD-Zoeppritz equations for Young's impedance and P-wave / S-wave velocity ratio, including: S01. Obtain the YIRD-Zoeppritz equation as follows: In the formula, intermediate parameters , , , Meanwhile, here , , and Existing parameters and ; PSR Indicates the ratio of longitudinal to transverse wave velocity; R PP , T PP , R PS and T PS These represent the P-wave reflection coefficient, P-wave transmission coefficient, converted S-wave reflection coefficient, and converted S-wave transmission coefficient, respectively; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. i 1 represents the angle of incidence of the reflected longitudinal wave; V p , V s and r These represent the longitudinal wave velocity, the transverse wave velocity, and the rock density, respectively. YI Indicates Young's impedance; S02. The ADA-MCMC-LpSC nonlinear inversion algorithm based on Lp norm to provide sparse prior information inverts Young's impedance. YI S-wave velocity ratio PSR and rock density .

[0018] In S01: In some embodiments, the method for obtaining the YIRD-Zoeppritz equation includes: S011. Based on the relationship between elastic parameters, both the longitudinal wave modulus and the transverse wave modulus can be expressed using Young's modulus and Poisson's ratio, as shown below: ; ; In the formula, M Represents the longitudinal wave modulus; m Indicates shear modulus; s Indicates Poisson's ratio; E This indicates Young's modulus; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. S012. The relationship between Poisson's ratio and the P-wave / S-wave velocity ratio is non-linear, as shown below: ; In the formula, s Indicates Poisson's ratio; PSR This represents the ratio of longitudinal to transverse wave speeds, i.e. PSR=V p / V s , V p and V s These represent the P-wave velocity and S-wave velocity, respectively; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. S013. Young's impedance is defined as the product of Young's modulus and rock density, and its expression is as follows: ; In the formula, the parameters YI Indicates Young's impedance; E Indicates Young's modulus; r This indicates the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. S014, Formula and Substitute into the formula respectively and The following equation can be obtained: ; ; In the formula, M Represents the longitudinal wave modulus; m Describe the shear modulus; define two parameters related to... PSR Intermediate parameters: and ; YI Indicates Young's impedance; r This indicates the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. Longitudinal wave velocity and transverse wave velocity can be expressed in the following form: ; ; In the formula, V p , V s and r These represent P-wave velocity, S-wave velocity, and rock density, respectively; intermediate parameters and ; PSR Indicates the ratio of longitudinal to transverse wave velocity; YI Indicates Young's impedance; r This indicates the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. S015. Based on the accurate Zoeppritz equation for the continuity of displacement and stress at the interface, including the ratio of longitudinal to transverse wave velocities. PSR and Young's impedance YI The new exact equation is expressed as follows, and this new equation is called the YIRD-Zoeppritz equation: In the formula, intermediate parameters , , , Meanwhile, here , , and Existing parameters and ; PSR Indicates the ratio of longitudinal to transverse wave velocity; R PP , T PP , R PS and T PS These represent the P-wave reflection coefficient, P-wave transmission coefficient, converted S-wave reflection coefficient, and converted S-wave transmission coefficient, respectively; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. i 1 represents the angle of incidence of the reflected longitudinal wave; V p , V s and r These represent the longitudinal wave velocity, the transverse wave velocity, and the rock density, respectively. YI This represents Young's impedance.

[0019] In S015: Furthermore, the exact Zoeppritz equation based on the continuity of displacement and stress at the interface is as follows: ; In the formula, R PP , T PP , R PS and T PS These represent the longitudinal wave reflection coefficient, longitudinal wave transmission coefficient, converted transverse wave reflection coefficient, and converted transverse wave transmission coefficient, respectively; intermediate parameters. , , and The subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively. i 1 represents the angle of incidence of the reflected longitudinal wave; V p , V s and r These represent the longitudinal wave velocity, the transverse wave velocity, and the rock density, respectively.

[0020] In S012: It is understandable that Poisson's ratio is a rock parameter that is highly sensitive to fluid content.

[0021] Understandable. PSR The ratio of longitudinal to transverse wave velocity.

[0022] In S013: It is understandable that Young's impedance is defined as the product of Young's modulus and rock density. This parameter is widely used in lithological discrimination and brittleness analysis in the identification of sweet spots in unconventional reservoirs.

[0023] In S02: It is understandable that the nonlinear inversion problem is solved within a Bayesian framework, where the likelihood function is used to quantify the degree of mismatch between observed and predicted seismic data, and based on L... p The prior probability distribution of the norm is used to enhance the sparsity of the solution.

[0024] In some embodiments, the ADA-MCMC-LpSC nonlinear inversion algorithm, which provides sparse prior information based on the Lp norm, includes the following steps: S021. Provide a likelihood function describing the probability density of seismic data: ; In the formula, (d|m) denotes the likelihood function; ||·||2 denotes the Euclidean 2 norm; the symbol exp denotes an exponential function with the natural constant e as the base; m The parameters to be inverted are represented by Young's impedance, the ratio of longitudinal to transverse wave velocity, and rock density.d Represents the seismic data matrix; R PP Represents the longitudinal wave reflection coefficient; W Represents the wavelet matrix; N For time sampling points; Here is the noise covariance matrix; symbol This is the value of the determinant, used for Gaussian distribution normalization; k The number of incident angles; S022. To improve the accuracy of the inversion results, L is introduced into the prior probability density function. p The norm regularization term, with its prior distribution, is represented as: ; In the formula, Indicates having L The prior probability density function of the p-norm sparse constraint; ||·|2 represents the Euclidean 2-norm; The prior mean of the model is usually taken from the filter well curve; ∑ m Represents the prior covariance matrix with spatial correlation; symbol is the value of the determinant; the symbol exp represents an exponential function with the natural constant e as its base; m This represents the parameters to be inverted, consisting of Young's impedance, the ratio of longitudinal to transverse wave velocity, and the rock density. N The time sampling point is denoted by exp; the symbol exp represents an exponential function with the natural constant e as its base. R PP Represents the longitudinal wave reflection coefficient; ||·|| p Indicates the use of sparsity enhancement L p-norm; l The sparsity regularization parameter; 0 < p <1 is used to control the degree of sparsity; S023. Based on the Bayesian framework, combining the likelihood function and the sparse constraint prior distribution, the posterior probability distribution is formally represented as follows: ; In the formula, The proposed means having L The posterior probability distribution of the p-norm sparse constraint, where ||·|2 represents the Euclidean 2-norm; ||·| p Indicates the use of sparsity enhancement L p-norm; m This represents the parameters to be inverted, consisting of Young's impedance, the ratio of longitudinal to transverse wave velocity, and the rock density. d Represents the seismic data matrix; W The symbol exp represents the wavelet matrix; the symbol exp represents an exponential function with the natural constant e as the base. The prior mean of the model is usually taken from the filter well curve; l For sparsity regularization parameters; Here is the noise covariance matrix; R PP Represents the longitudinal wave reflection coefficient; ||·|| p Indicates the use of sparsity enhancement L p-norm; 0 < p <1 is used to control the degree of sparsity; the symbol exp represents an exponential function with the natural constant e as the base; S024. To further improve computational efficiency and accurately estimate the posterior probability distribution, a two-stage acceptance mechanism for the delayed acceptance Markov chain Monte Carlo algorithm is introduced. S025, finally obtaining Young's impedance YI, PSR (longitudinal and transverse wave velocity ratio), and rock density. The inversion value of the parameter.

[0025] It is understandable that the nonlinear inversion problem is solved within a Bayesian framework, where the likelihood function is used to quantify the degree of mismatch between observed and predicted seismic data, and based on L... p The prior probability distribution of the norm is used to enhance the sparsity of the solution. Specifically, the observation error is modeled as having a mean of zero and a variance of ∑ N Given Gaussian white noise, the likelihood function of the inversion model, under the given forward model conditions, is expressed as: ; In the formula, Let be the likelihood function; ||·|2 represents the Euclidean 2 norm; Here is the noise covariance matrix; N For time sampling points; k The number of incident angles; Here is the noise covariance matrix; symbol is the value of the determinant; the symbol exp represents an exponential function with the natural constant e as its base; d Represents the seismic data matrix; m This represents the parameters to be inverted, consisting of Young's impedance, the ratio of longitudinal to transverse wave velocity, and the rock density. d Represents the seismic data matrix; W Represents the wavelet matrix; R PP This represents the longitudinal wave reflection coefficient.

[0026] It is understandable that, in order to improve the accuracy of the inversion results, L is introduced into the prior probability density function. p The norm regularization term can be viewed as a practical relaxation of the sparse constraint of the ideal L0 norm.

[0027] It is understandable that the aforementioned posterior distribution forms the theoretical basis of the inversion strategy, enabling simultaneous processing of uncertainty quantification, noise suppression, and automatic sparsity constraints. This framework provides an inversion method for obtaining Young's impedance (YI) and PSR profiles with higher geological plausibility and resolution.

[0028] In S021: In some embodiments, for a given angle of incidence i k Its synthetic seismic data is represented as follows: ; In the formula, d This represents forward modeling seismic data after noise has been added. Indicates noise; subscript N Indicates the impedance of Young's method YI and PSR The number of sampling points with equal parameters; i k This represents the k-th angle of incidence; the subscript k indicates the index of the angle of incidence. W Indicates the relationship with the given first k An incident angle i k Related wavelet matrices; R PP Represents the longitudinal wave reflection coefficient; in W The definition is as follows: ; In the formula, W Indicates the relationship with the given first k An incident angle i k Related wavelet matrices; w Indicates the angle of incidence i k Seismic wavelet; subscript th Indicates the length of the seismic wavelet; matrix subscript N × N This indicates that the matrix contains N rows and N columns of elements.

[0029] It is understandable that, firstly, the convolution model assumes that seismic data related to the incident angle can be generated by convolving the reflection coefficient with wavelets at different incident angles. Therefore, for a given incident angle... i k Its synthetic seismic data can be represented as: ; In the formula, d This represents forward modeling seismic data after noise has been added. i k Indicates the k-th angle of incidence;W Indicates the relationship with the given first k An incident angle i k Related wavelet matrices; R PP Represents the longitudinal wave reflection coefficient; Represents noise. Matrix subscripts. N × 1 This indicates that the matrix contains N rows and 1 column of elements.

[0030] In S024: This is understandable, using the Delayed Acceptance Markov Chain Monte Carlo algorithm (DA-MCMC).

[0031] Furthermore, let the current state be... And by candidate models To generate a new state, assuming a symmetric proposal distribution, the acceptance probability in the first stage is expressed as: ; In the formula, This represents the acceptance probability in the first stage. Represents the model parameter vector The prior probability distribution; This is the model parameter vector for the current (i-1)th iteration state; The candidate model parameters generated for the proposal distribution; symbol min{1, } indicates taking a value less than 1; d Represents the seismic data matrix; Let be the likelihood function, representing the candidate model parameter vector. The degree of matching between the seismic data obtained from forward modeling and the observed seismic data; Indicates the current model parameters The degree of matching between the seismic data obtained from forward modeling and the observed seismic data; For the current model Generate candidate models The distribution of proposals; To select candidate models Generate the current model The distribution of proposals, in the above formula = ; This represents the prior probability density of the candidate model; This represents the prior probability density of the current model.

[0032] It is understood that the distribution of the proposals in this application is a symmetrical Gaussian distribution, therefore in the above formula... = ; This represents the prior probability density of the candidate model; This represents the prior probability density of the current model.

[0033] Furthermore, from a uniform distribution Generate a random number ,if If the proposal is rejected, the Markov chain remains in its current state. ; if Then candidate models are needed. With probability The second phase of evaluation involves the final acceptance of a sample based on its... L The second-stage acceptance rate of the exact target probability distribution with p-norm sparse constraints The calculation method is as follows: ; From uniform distribution Generate a random number ,if If the proposal is rejected, the Markov chain remains in its current state. ;if If so, then accept the proposal, that is... ; This indicates the effect of introducing the Lp regularization term on the candidate model. Posterior probability density; For the current model The posterior probability density; For the current iterative inversion i The model parameter vector of the -1st order state; Candidate model parameters generated for the proposal distribution; d This represents the earthquake data matrix.

[0034] Furthermore, an adaptive proposal distribution consisting of two parts, proposed by Roberts and Rosenthal in 2009, is introduced to obtain candidate models. Its expression is as follows: ; In the formula, Represents a multivariate normal distribution; Represents a unit diagonal matrix, used in ∑ m Even with a poor condition number, the existence of proposed samples can still be guaranteed; N For time sampling points; ∑ m Represents the prior covariance matrix with spatial correlation; This represents the candidate model parameter vector consisting of Young's impedance, the ratio of longitudinal to transverse wave velocity, and the rock density. For the current iterative inversion i The model parameter matrix for the -1st order state; parameters β The mixed weight coefficients are in the range [0,1], and their update method is as follows: ; In the formula, the parameters β The mixed weighting coefficients are in the range [0,1]. Q This indicates the acceptance rate of the MCMC Markov chain.

[0035] It is understandable that the acceptance rate is calculated as the percentage of accepted iterations out of every 1000 iterations.

[0036] Understandable. This represents a diagonal matrix (usually the identity matrix).

[0037] It is understandable that the adaptive proposal distribution form, consisting of two parts, proposed by Roberts and Rosenthal (2009), is introduced, and its expression is as follows: ; In the formula, Represents a multivariate normal distribution; Denotes a diagonal matrix, used in ∑ m Even with a poor condition number, the existence of proposed samples can still be guaranteed; N For time sampling points; ∑ m Represents the prior covariance matrix with spatial correlation; m This represents the parameters to be inverted, consisting of Young's impedance, the ratio of longitudinal to transverse wave velocity, and the rock density. For the current iterative inversion i The model parameter matrix for the -1st order state; The candidate model parameters generated for the proposal distribution.

[0038] In S025: Furthermore, after the aforementioned multiple iterative inversion processes and evaluation and acceptance through a two-stage acceptance mechanism, the model parameters... After the iterative burning period, the Markov chain finally reaches its final state on the second iteration. i When the iteration reaches t0, it converges to a stationary phase. Finally, the mean of all model parameter values ​​retrieved during the iterative inversion from the stationary phase t0 to the total number of iterations is calculated to obtain the impedance from Young's equation. YI S-wave velocity ratio PSR and rock density Inversion parameters composed of parameters m .

[0039] In traditional methods, the first approach involves using the predicted P-wave velocity, S-wave velocity, and rock density obtained through inversion from the Aki-Richards equations, and then indirectly calculating Young's impedance. YI and the ratio of longitudinal to transverse wave speeds PSR The second method directly incorporates Young's impedance, Poisson's ratio, and rock density into the linear approximation equation (called the YIPD equation), enabling direct inversion of Young's impedance. The third method directly inverts the P-wave velocity ratio (PSR) using a linear approximation equation composed of P-wave velocity, PSR, and rock density (called the PRD equation). The specific forms of the linear approximation equations mentioned in these three traditional methods are expressed as follows: ① Aki-Richards equation: ; In the formula, R PP Represents the longitudinal wave reflection coefficient; parameter k It represents the square of the ratio of the transverse wave velocity to the longitudinal wave velocity; i Indicates the incident angle of the reflected longitudinal wave; V p , V s and r These represent P-wave velocity, S-wave velocity, and rock density, respectively; the overline "-" indicates the average value of the parameters of the upper and lower strata; the symbol Δ indicates the amount of variation between the parameters of the lower and upper strata.

[0040] ② YIPD equation: ; In the formula, R PP Represents the longitudinal wave reflection coefficient; parameter k It represents the square of the ratio of the transverse wave velocity to the longitudinal wave velocity; i Indicates the incident angle of the reflected longitudinal wave; YI , s and r These represent Young's impedance, Poisson's ratio, and rock density, respectively; the overline "-" indicates the average value of the parameters of the upper and lower strata; the symbol Δ indicates the amount of variation between the parameters of the lower and upper strata.

[0041] ③ PRD equation: ; In the formula, R PP The value represents the longitudinal wave reflection coefficient; the parameter k represents the square of the ratio of the transverse wave velocity to the longitudinal wave velocity. i Indicates the incident angle of the reflected longitudinal wave; PSR , Vp and r These represent the ratio of P-wave velocity to S-wave velocity, P-wave velocity, and rock density, respectively; the underline "-" indicates the average value of the corresponding parameters of the upper and lower strata; the symbol Δ indicates the amount of variation between the parameters of the lower and upper strata.

[0042] To verify the accuracy advantage of the proposed YIRD-Zoeppritz equation compared to the three traditional linear approximation equations mentioned above, this application comprehensively considers different types of AVO models, corresponding to model numbers I, II, III, and IV in the table below. Under different incident angles, the computational accuracy and error of these equations are compared with the standard Zoeppritz equation. This application designs four two-layer models to characterize different AVO features, with model II further subdivided into cases A and B. The relevant elastic parameters and rock density are listed in Table 1.

[0043] Table 1. Relevant parameters of four types of two-level models

[0044] exist Figure 1 In the study, several approximate equations (Aki-Richards, YIPD, and PRD) showed similar trends and small deviations from the exact Zoeppritz equation (black curve) when the incident angle was less than approximately 30°, demonstrating good consistency. However, when the incident angle exceeded 30°, the curves of the linear approximation equations began to deviate significantly from the exact Zoeppritz equation. Figure 2 In the case of Type II-A AVO, the variation trends of all approximate equations within a 50° incident angle range are roughly consistent with the exact Zoeppritz equations. However, the errors between the linear approximation equations and the exact Zoeppritz equations still exist and cannot be avoided. Figure 3 In the Type III AVO shown, the response error of the PRD approximation equation is most pronounced at large incident angles, which may reduce the accuracy of inverting PSR parameters using the PRD approximation equation. Figure 4 In Type IV AVO models, the PRD approximation equation did not demonstrate an accuracy advantage over other approximation equations when the incident angle exceeded 30°. Overall, the accuracy of different approximation equations varied significantly across different AVO types, while the novel equation proposed in this application exhibited a clear advantage in both highest accuracy and lowest error. In contrast, the YIRD-Zoeppritz equation proposed in this application achieved accuracy closer to the exact Zoeppritz equation in all four AVO types, with the error approaching zero (e.g., ...). Figure 1 (b)–4(b)). Therefore, the YIRD-Zoeppritz equation proposed in this application is Young's impedance. YIThe direct inversion of the PSR (Positioning-to-Side Wave Ratio) lays the foundation for the equations.

[0045] Figure 1 In Figure (a), the different colored curves represent the curves of the reflection coefficient of different equations as a function of the incident angle. Figure (b) shows the curves of the error between the three traditional linear approximation equations (Aki-Richards equation, YIPD equation and PRD equation) and the new equation YIRD-Zoeppritz equation and the exact Zoeppritz equation as a function of the incident angle. Figure 2 In Figure (a), the different colored curves represent the curves of the reflection coefficient of different equations as a function of the incident angle. Figure (b) shows the curves of the error between the three traditional linear approximation equations (Aki-Richards equation, YIPD equation and PRD equation) and the new equation YIRD-Zoeppritz equation and the exact Zoeppritz equation as a function of the incident angle. Figure 3 In Figure (a), the different colored curves represent the curves of the reflection coefficient of different equations as a function of the incident angle. Figure (b) shows the curves of the error between the three traditional linear approximation equations (Aki-Richards equation, YIPD equation and PRD equation) and the new equation YIRD-Zoeppritz equation and the exact Zoeppritz equation as a function of the incident angle. Figure 4 In Figure (a), the different colored curves represent the curves of the reflection coefficient of different equations as a function of the incident angle. Figure (b) shows the curves of the error between the three traditional linear approximation equations (Aki-Richards equation, YIPD equation and PRD equation) and the new equation YIRD-Zoeppritz equation and the exact Zoeppritz equation as a function of the incident angle.

[0046] from Figure 5 and Figure 6 The logging curve inversion test results show that, under both noise-free and signal-to-noise ratio (SNR) conditions of 5, the P-wave velocity ratio, Young's impedance, and rock density can all be effectively solved by the ADA-MCMC-LpSC nonlinear optimization inversion algorithm. The inversion results (red) of each parameter are in good agreement with their corresponding logging curves (black).

[0047] Application examples To verify the effectiveness of the method protected in this application in practical applications, and to demonstrate its ability to effectively extract reservoir parameters, specifically Young's impedance, from seismic data. YI S-wave velocity ratio PSR and rock density rFor a specific shale field in China, a longitudinal survey line passing through well A was selected and angles were superimposed according to different incident angles, such as... Figure 7 As shown, this specifically includes incident angles of 3° (e.g.) Figure 7 a) 15° (e.g.) Figure 7 b) and 27° (as shown in b) Figure 7 c). Then, based on the precise YIRD-Zoeppritz equation proposed above, and combined with the proposed ADA-MCMC-LpSC nonlinear inversion algorithm based on sparse prior information provided by the Lp norm, seismic data is used as input data, and the parameter prediction results corresponding to the seismic profile are obtained through iterative inversion, such as... Figure 8 As shown, the proposed method can effectively and accurately retrieve the PSR, Young's impedance, and rock density parameter profiles from actual seismic data. Notably, low PSR values ​​can be used to characterize good gas-bearing properties in shale reservoirs, as Young's impedance and density are also low. The retrieved PSR, Young's impedance, and rock density results are consistent with the corresponding curve values ​​shown in Well A. Figure 7 The black solid curves shown in the profiles of each parameter are the corresponding well logging curves, indicating that the proposed method is effective and reliable.

[0048] 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 nonlinear direct seismic inversion method for Young's impedance and P-S wave velocity ratio based on the exact YIRD-Zoeppritz equation, characterized in that, Comprise: Obtaining YIRD-Zoeppritz equation, as follows: where the intermediate parameters , , , ; while here , , and there are the parameters and ; PSR denotes the ratio of P- to S-wave velocities; R PP , T PP , R PS and T PS denote the P-wave reflection coefficient, the P-wave transmission coefficient, the converted S-wave reflection coefficient and the converted S-wave transmission coefficient, respectively; the subscripts "1" and "2" denote the parameters of the upper and lower strata, respectively; Theta 1 denotes the incidence angle of the reflected P-wave; V p , V s and Rho denote the P-wave velocity, the S-wave velocity and the rock density, respectively; YI denotes the Young's impedance; The ADA-MCMC-LpSC nonlinear inversion algorithm based on Lp norm provides sparse prior information to invert the Young's impedance YI , the ratio of P-wave and S-wave velocity PSR , and the rock density .

2. The non-linear direct seismic inversion method based on the precise YIRD- Zoeppritz equation of Young's impedance and P-S wave velocity ratio according to claim 1, characterized in that, Obtaining method of YIRD-Zoeppritz equation comprises: According to the relationship between elastic parameters, both longitudinal wave modulus and transverse wave modulus can be expressed by Young's modulus and Poisson's ratio, as shown below: ; ; wherein M represents the longitudinal wave modulus; Mu represents the shear modulus; Sigma represents the Poisson's ratio; E represents the Young's modulus; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively; Poisson's ratio and longitudinal-transverse wave velocity ratio are in nonlinear relationship, as shown below: ; wherein Sigma represents the Poisson's ratio; PSR represents the ratio of P-wave velocity to S-wave velocity, i.e. PSR=V p V s V p and V s represents the P-wave velocity and S-wave velocity, respectively; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively.​​ Young's impedance is defined as the product of Young's modulus and rock density, and its expression is as follows: ; where the parameters YI represents the Young's impedance; E represents the Young's modulus; Rho represents the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively; Substituting the equations and into the equations and respectively, the following equations are obtained: ; ; wherein M represents the longitudinal wave modulus; Mu represents the shear modulus; defines two intermediate parameters of PSR and ; YI represents the Young's impedance; Rho represents the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively;​ Longitudinal wave velocity and transverse wave velocity can be expressed in the following form: ; ; wherein, V p , V s and Rho represent the P-wave velocity, the S-wave velocity and the rock density, respectively; the intermediate parameters and ; PSR represent the P-S wave velocity ratio; YI represent the Young's impedance; Rho represent the rock density; the subscripts "1" and "2" represent the parameters of the upper and lower strata, respectively; Based on the exact Zoeppritz equation for displacement and stress continuity at the interface, which contains the ratio of P- to S-wave velocities PSR and Young's impedance YI The new exact equation is given by YIRD-Zoeppritz equation where the intermediate parameters , , , ; while here , , and the parameters and ; PSR represent the ratio of P- to S-wave velocities; R PP , T PP , R PS and T PS represent the P-wave reflection coefficient, the P-wave transmission coefficient, the converted S-wave reflection coefficient and the converted S-wave transmission coefficient, respectively; the subscripts "1" and "2" represent the parameters of the overburden and the underlying formation, respectively; Theta 1 represents the angle of incidence of the reflected P-wave; V p , V s and Rho represent the P-wave velocity, the S-wave velocity and the rock density, respectively; YI represents the Young's impedance.

3. The nonlinear direct seismic inversion method based on the exact YIRD- Zoeppritz equation for Young's impedance and P-S wave velocity ratio according to claim 2, characterized in that, Based on the continuity of displacement and stress at the interface, the accurate Zoeppritz equation is as follows: ; wherein R PP , T PP , R PS and T PS respectively represent the reflection coefficient of the longitudinal wave, the transmission coefficient of the longitudinal wave, the reflection coefficient of the converted transverse wave and the transmission coefficient of the converted transverse wave; the intermediate parameters , , and ; the subscripts "1" and "2" respectively represent the parameters of the upper and lower strata; Theta 1 represents the incident angle of the reflected longitudinal wave; V p , V s and Rho respectively represent the longitudinal wave velocity, the transverse wave velocity and the rock density.

4. The nonlinear direct seismic inversion method based on the exact YIRD- Zoeppritz equation of Young's impedance and P-S wave velocity ratio of claim 3, wherein, Comprise: ADA-MCMC-LpSC nonlinear inversion algorithm based on Lp norm providing sparse prior information comprises the following steps: Providing a likelihood function describing the probability density of seismic data: ; wherein, (d | m) denotes the likelihood function; || · ||2denotes the Euclidean 2-norm; the symbol exp denotes the exponential function with base the natural constant e; m denotes the parameters to be inverted consisting of the Young's impedance, the P-S wave velocity ratio and the rock density; d denotes the seismic data matrix; R PP denotes the P-wave reflection coefficient; W denotes the wavelet matrix; N is the time sampling point; is the noise covariance matrix; the symbol is the value of the determinant used for the Gaussian distribution normalization; k is the number of incidence angles; To improve the precision of inversion results, L p norm regularization term is introduced into the prior probability density function, and the prior distribution is expressed as: ; wherein, represents the prior covariance matrix with spatial correlation; the symbol L is the prior probability density function with sparsity constraint of p-norm; ||·||2represents the Euclidean 2-norm; is the model prior mean, usually taken from the filtered well curve; ∑ m represents the prior covariance matrix with spatial correlation; the symbol is the value of the determinant; the symbol exp represents the exponential function with base of the natural constant e; m represents the parameter to be inverted composed of Young's impedance, P-wave to S-wave velocity ratio and rock density; N is the time sampling point; the symbol exp represents the exponential function with base of the natural constant e; R PP represents the P-wave reflection coefficient; ||·||p p represents the sparsity enhancing function L p-norm; Lambda is the sparsity regularization parameter; 0 p <1 is used to control the degree of sparsity; Based on the Bayesian framework, the likelihood function and the sparse constraint prior distribution are integrated, and the posterior probability distribution is expressed as follows: ; wherein, denotes the proposed posterior probability distribution with L p-norm sparsity constraint; ‖·‖2denotes the Euclidean 2-norm; ‖·‖ p denotes the p-norm used to enhance sparsity L p-norm; m denotes the parameters to be inverted consisting of Young's impedance, P / S velocity ratio and rock density; d denotes the seismic data matrix; W denotes the wavelet matrix; the symbol exp denotes the exponential function with base the natural number e; is the model prior mean, usually taken from filtered well curves; Lambda is the sparsity regularization parameter; is the noise covariance matrix; R PP denotes the P-wave reflection coefficient; ‖·‖ p denotes the p-norm used to enhance sparsity L p-norm; 0 p <1 controls the degree of sparsity; the symbol exp denotes the exponential function with base the natural number e; In order to further improve the calculation efficiency and accurately estimate the posterior probability distribution, a two-stage acceptance mechanism of delayed acceptance Markov chain Monte Carlo algorithm is introduced; The final results of Young's impedance YI, the ratio of longitudinal and transverse wave velocities PSR and rock density The inversion values of the parameters.

5. The non-linear direct seismic inversion method based on the exact YIRD- Zoeppritz equation for Young's impedance and P-S wave velocity ratio as claimed in claim 4, characterized in that, For a given angle of incidence Theta k whose synthetic seismic data is represented by: ; wherein d denotes the forward seismic data after adding noise; denotes noise; subscript N denotes the Young's impedance YI and PSR the number of sampling points equal to the parameter; Theta k denotes the kth incident angle; subscript k denotes the serial number of the incident angle; W denotes the wavelet matrix k related to the given kth incident angle Theta k ; subscript k denotes the serial number of the incident angle; R PP denotes the P-wave reflection coefficient; wherein W are defined as follows: ; wherein W represents the wavelet matrix associated with a given incidence angle k of the nth trace; Theta k of the nth trace; w represents the seismic wavelet for an incidence angle Theta of the nth trace; k of the nth trace; Th represents the length of the seismic wavelet; the matrix subscript N × N represents that the matrix contains N rows and N columns of elements.

6. The non-linear direct seismic inversion method for Young's impedance and P-S wave velocity ratio based on the exact YIRD-Zoeppritz equation of claim 4, wherein The current state is and a new state is generated by the candidate model Under the assumption that the proposal distribution is symmetric, the acceptance probability of the first stage is given by ; where is the acceptance probability of the first stage; is the prior probability distribution of the model parameter vector ; is the model parameter vector of the current i-1th iteration state; is the candidate model parameter generated by the proposal distribution; the symbol min{1, } means taking the smaller value than 1; d is the seismic data matrix; is the likelihood function, indicating the matching degree of the seismic data obtained by the forward modeling of the candidate model parameter vector ; is the matching degree of the seismic data obtained by the forward modeling of the current model parameter ; is the proposal distribution for generating the candidate model from the current model ; is the proposal distribution for generating the current model from the candidate model , and in the above formula = ; is the prior probability density of the candidate model; is the prior probability density of the current model.

7. The nonlinear direct seismic inversion method based on the precise YIRD- Zoeppritz equation of Young's impedance and P-S wave velocity ratio according to claim 6, characterized in that, from a uniform distribution generating a random number , if , reject the offer, the Markov chain remains in the current state, i.e. ; If , then the candidate model is accepted with probability L , the sample enters the second phase, and its final acceptance or rejection is determined based on the acceptance rate of the second phase with the exact target probability distribution under the p-norm sparsity constraint , which is calculated as follows: ; From uniform distribution Generate a random number ,if If the proposal is rejected, the Markov chain remains in its current state. ;if If so, then accept the proposal, that is... ; This indicates that after introducing the Lp regularization term, the candidate model... Posterior probability density; For the current model The posterior probability density; For the current iterative inversion i The model parameter vector of the -1st order state; Candidate model parameters generated for the proposal distribution; d This represents the earthquake data matrix.

8. The nonlinear direct seismic inversion method based on the precise YIRD- Zoeppritz equation of Young's impedance and P-S wave velocity ratio of claim 7, wherein, The form of the adaptive proposal distribution consisting of two parts introduced by Roberts and Rosenthal in 2009 is brought in to obtain the candidate model whose expression is as follows: ; where denotes a multivariate normal distribution; denotes an identity matrix, which is used to make m ensure the existence of the proposed sample even when the condition number is poor; N is the time sampling point;∑ m denotes the prior covariance matrix with spatial correlation; denotes the candidate model parameter vector consisting of Young's impedance, P-wave to S-wave velocity ratio, and rock density; is the model parameter matrix of the current iteration i -1th state; parameter β is the mixing weight coefficient between [0, 1], which is updated as follows: ; where the parameters β are mixing weight coefficients between [0, 1]; Q denotes the acceptance rate of the MCMC Markov chain.

9. The nonlinear direct seismic inversion method based on the precise YIRD- Zoeppritz equation of Young's impedance and P-S wave velocity ratio of claim 8, wherein, After the aforementioned multiple iterative inversion processes and the evaluation and acceptance through a two-stage acceptance mechanism, the model parameters... After the iterative burning period, the Markov chain finally reaches its final state on the second iteration. i When the iteration reaches t0, it converges to a stationary phase. Finally, the mean of all model parameter values ​​retrieved during the iterative inversion from the stationary phase t0 to the total number of iterations is calculated to obtain the impedance from Young's equation. YI S-wave velocity ratio PSR and rock density Inversion parameters composed of parameters m .