A method for inversion of surface wave dispersion probability based on approximate Bayesian calculation

Abstract

The present invention relates to the technical field of surface wave dispersion inversion, and provides a surface wave dispersion probability inversion method based on approximate Bayesian calculation, comprising: 1) deriving a Bayesian equation for surface wave dispersion probability inversion under the framework of approximate Bayesian calculation ABC; 2) determining the constraint of parameters through surface wave dispersion data, thereby specifying the prior distribution of parameters; 3) assuming the distribution of observed data and specifying a likelihood function, eliminating the observed data error term in the likelihood function; 4) determining a distance function based on the specific form of the likelihood function; 5) deriving an analytical formula for quickly estimating the distance function value based on the distance function; 6) solving the Bayesian equation based on the Bayesian equation and the formula for analytically calculating the distance function, using ABC and combining subset simulation sampling; 7) quantifying the uncertainty of the inversion result based on the posterior sample. The present invention can better quantify the uncertainty of the shear wave velocity profile obtained by inverting the surface wave dispersion.

Description

Technical Field

[0001] The invention relates to the technical field of surface wave dispersion inversion, and in particular to a surface wave dispersion probability inversion method based on approximate Bayesian calculation. Background Art

[0002] The dispersion characteristics of surface waves refer to the phenomenon that when surface waves propagate on the surface, wavelengths of different frequencies are attenuated to different degrees and their speed changes during propagation. This phenomenon that the propagation speed (phase velocity) of surface waves changes with frequency is called dispersion, and the dispersion characteristics can be studied through frequency-wavenumber spectrum analysis.

[0003] When studying the dispersion characteristics of surface waves, the commonly used method is the surface wave analysis method. By deploying multiple seismographs on the surface, seismic records at different locations can be obtained. Then, the seismic records are processed using the frequency-wavenumber spectrum analysis method to obtain the phase velocity of the surface wave at each frequency, that is, the surface wave dispersion curve. Since there is a close relationship between the frequency, wavelength and stratigraphic structure of seismic waves, information about the stratigraphic structure can be inverted by analyzing the dispersion characteristics of surface waves.

[0004] Multi-channel surface wave analysis (MASW) is a method for estimating the shear wave velocity (v s ) method. Due to its economic, flexible and reliable characteristics, the MASW method has been widely used in various fields of geotechnical engineering. The field seismic records obtained by the MASW method include the test site v s The information about the variation with depth can be obtained by inverting the surface wave dispersion curve extracted from the seismic record. s Due to measurement errors, finite dispersion data and other factors, the inverted v s There is a lot of uncertainty in the profile.

[0005] Surface wave dispersion inversion methods can be divided into two categories: deterministic methods and probabilistic methods. Deterministic methods, such as genetic algorithms and differential search algorithms, are usually unable to quantify the uncertainty of the inversion results. Probabilistic methods assume that the parameters to be inverted are random variables, and the uncertainty of parameter identification can be quantified based on the inversion results. Bayesian inversion is the most commonly used probabilistic inversion method. In the Bayesian inversion problem of surface wave dispersion, evaluating the likelihood function is very time-consuming. Specifically, it is usually assumed that the observed phase velocity is equal to the sum of the theoretical phase velocity and the random error. When the data dimension is M (that is, there are M observed phase velocities at different frequencies), evaluating the likelihood function once requires calculating the theoretical phase velocities of M different frequencies. The theoretical phase velocity is often obtained by numerically solving the root of the dispersion curve characteristic equation G=0 by bisection. Calculating the G value itself is time-consuming because a large number of nonlinear equation calculations and matrix multiplication operations are required. Each time the theoretical phase velocity is calculated, the characteristic equation G is called N times. classic ≥(c max -cmin ) / dv+1+log2(dv / ΔE), where c max and c min are the upper and lower bounds of the search range, dv is the search step size, and ΔE is the phase velocity solution accuracy and is usually set to 0.001. For example, when assuming c max =300m / s,c min =100m / s, dv=4, N classic ≥50; when dv is reduced to 2, N classic ≥86. The number of times the characteristic equation G is called to calculate the likelihood function is approximately equal to M×N classic , the process is very time-consuming.

[0006] Approximate Bayesian Computation (ABC) method is often used to solve Bayesian equations where it is difficult or impossible to calculate the likelihood function. The basic principle is to sample from the joint posterior distribution of the simulated samples of the observed data and the model parameters to avoid directly calculating the likelihood function value; but it requires that simulated samples can be easily extracted from the likelihood function, and a "distance function" is used to evaluate the distance between the simulated samples and the observed data to determine whether the two are approximately equal. For the problem of surface wave dispersion inversion, it is often assumed that the observed phase velocity is equal to the sum of the theoretical phase velocity and the random error. Therefore, under the existing ABC method framework, it is still necessary to numerically solve the theoretical phase velocity (the number of times G is called is approximately equal to M×N classic ) is substituted into the distance function to determine whether the simulated sample is equal to the observed dispersion data. Although the ABC method provides a potential framework for solving the time-consuming likelihood function evaluation problem in the Bayesian inversion of surface wave dispersion, new distance functions are still needed to avoid or reduce the calculation of theoretical phase velocities.

[0007] On the other hand, Markov chain-Monte Carlo (MCMC) sampling is often used to solve the Bayesian equation under the ABC framework. For the surface wave dispersion inversion problem, because the dispersion curve characteristic equation G is highly nonlinear and the dimensions of the model parameters and the dimensions of the observed data are both high-dimensional, MCMC often faces problems such as the difficulty of the Markov chain to reach a steady state and low sample acceptance rate when sampling from the joint posterior distribution of simulated samples and parameters, resulting in low solution efficiency and inaccurate parameter estimation.

[0008] In summary, there are two main problems when applying the ABC method to the probability inversion of surface wave dispersion:

[0009] (1) In the framework of the ABC method, the theoretical phase velocity value needs to be numerically solved (the number of times G is called is approximately equal to M×N classic ) to evaluate the distance function, which is a very time-consuming process.

[0010] (2) For the inversion problem of surface wave dispersion, the MCMC method has problems such as the difficulty of the Markov chain to reach a steady state and the low sample acceptance rate when solving this type of Bayesian equation. Summary of the invention

[0011] The present invention provides a surface wave dispersion probability inversion method based on approximate Bayesian calculation, which can overcome certain defects of the prior art.

[0012] A surface wave dispersion probability inversion method based on approximate Bayesian calculation according to the present invention comprises the following steps:

[0013] (1) Derive the Bayesian equation for the probability inversion of surface wave dispersion under the framework of approximate Bayesian computing (ABC);

[0014] (2) Determine the parameter constraints through surface roll dispersion data to specify the parameter prior distribution;

[0015] (3) Assume the distribution of observed data and specify the likelihood function, and eliminate the observed data error term in the likelihood function by integration;

[0016] (4) based on the specific form of the likelihood function in step (3), determining a distance function for evaluating the distance between the observed data and the simulated samples;

[0017] (5) Based on the distance function in step (4), derive an analytical formula for quickly estimating the value of the distance function;

[0018] (6) Based on the Bayesian equation of step (1) and the formula for analytically calculating the distance function of step (5), ABC is used in combination with subset simulation sampling to solve the Bayesian equation;

[0019] (7) Quantify the uncertainty of the inversion results based on the parameter posterior samples in step (6).

[0020] As a preferred method, in step (1), the shear wave velocity v of the test site is approximated by a horizontal layered model. s The number of layers in the horizontal layered model is N, and the model parameters to be inverted are the thickness H and v of each layer. s ; Given the observed dispersion curve d obs Data, the probability of θ in the ABC framework is expressed by its posterior distribution p(θ|d obs ) quantized, solving the posterior distribution p(θ|d obs ) can be converted into solving the joint posterior distribution of θ and the simulated sample y of the observed dispersion curve data, that is:

[0021]

[0022] where p(y|θ) is the prior distribution of y given θ, i.e., d obsThe distribution it obeys; p(θ) is the prior distribution of θ; p ε (d obs |y) is the likelihood function; P ε (d obs ) is the normalization constant; in ABC, the likelihood function p ε (d obs |y) has an indicator function I S(ε) In the form of, when y is in the set S(ε)={L(d obs -y)≤ε} S(ε) (y)=1, then it is considered that y and d obs Approximately equal; otherwise, y and d obs Not equal; L is the calculation of y and d obs The distance function of the distance between .

[0023] Preferably, in step (2), the maximum identifiable depth and the minimum identifiable layer thickness of the shear wave velocity interface are determined based on the observed surface wave dispersion data to constrain the parameter value range and specify the prior distribution.

[0024] Preferably, in step (3), it is assumed that the observed surface wave dispersion data obeys a Gaussian distribution with a mean of 0 and an unknown variance, so the likelihood function can be specified as a multivariate Gaussian distribution. The unknown variance term in the likelihood function is eliminated by an integral method to reduce the number of inversion parameters.

[0025] Preferably, in step (4), the Euclidean distance is used to quantitatively evaluate the surface wave dispersion observation data d obs The distance between the simulated sample y.

[0026] Preferably, in step (5), the ratio G / G′ of the dispersion curve characteristic equation G and its first-order derivative G′ with respect to the phase velocity c is used to construct a fast estimation of the observed data d obs The Euclidean distance between the simulated sample y is used to estimate the distance function value.

[0027] Preferably, in step (6), the nested threshold decreasing sequence ε in the subset simulation is selected 1 =+∞>ε 2 >…>ε J =ε area as the observed dispersion data d in the observation space obs The middle region corresponding to the increasingly close approximation of S(ε t )={L(d obs -y)≤ε t}; Use the ABC method to sample from these intermediate regions to efficiently solve the posterior distribution p ε (θ,y|d obs ).

[0028] The beneficial effects of the present invention are:

[0029] The Bayesian method is used to invert the surface wave dispersion curve to obtain the shear wave velocity profile, and the uncertainty of the inversion result can be reasonably quantified based on the posterior samples; the Bayesian equation for surface wave dispersion inversion is evaluated based on the approximate Bayesian method, which can avoid the time-consuming likelihood function calculation in the traditional Bayesian method, and at the same time, a function for quickly estimating the distance between the observed phase velocity and the simulated sample is proposed, which reduces the time consumption for evaluating the Bayesian equation; the high-dimensional and nonlinear Bayesian equation in the surface wave dispersion inversion problem is efficiently evaluated using the approximate Bayes combined with the subset simulation sampling method, and the nested decreasing sequence regions in the subset simulation are selected as the regions corresponding to the increasingly close approximations of the observed dispersion data vector, so as to efficiently sample from the posterior distribution. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 It is a flow chart of a method for inverting the probability of surface wave dispersion based on approximate Bayesian calculation in an embodiment;

[0031] Figure 2(a) shows the virtual field v in case 1 of the embodiment. s Schematic cross-section;

[0032] FIG2( b ) is a schematic diagram of dispersion curve data synthesized in Case 1 of the embodiment;

[0033] Figure 3 The figure is a schematic diagram of the parameter posterior samples and normalized histogram obtained by the method proposed in Case 1 of the embodiment, as well as the results obtained by using the accurate distance function and the accurate Bayesian analysis results;

[0034] FIG4(a) is a diagram of the estimated distance function of case 1 in the embodiment. y Schematic diagram of sample and synthetic dispersion data;

[0035] Figure 4(b) shows the method for obtaining v in case 1 of the embodiment. s Schematic diagram of the cross-section results;

[0036] Figure 5 This is a schematic diagram of the actual observed dispersion data used in Case 2 in the embodiment;

[0037] Figure 6 The figure is a schematic diagram of the parameter posterior samples and normalized histogram obtained by the method proposed in Case 2 of the embodiment, as well as the results obtained by using the accurate distance function and the accurate Bayesian analysis results;

[0038] FIG. 7( a ) is a diagram of the estimated distance function of case 2 in the embodiment. y Schematic diagram of sample and observed dispersion data;

[0039] Figure 7(b) shows the method for obtaining v in case 2 of the embodiment. s Schematic diagram of the cross-section results. DETAILED DESCRIPTION

[0040] In order to further understand the content of the present invention, the present invention is described in detail in conjunction with the accompanying drawings and embodiments. It should be understood that the embodiments are only for explaining the present invention and are not intended to limit it.

[0041] Example

[0042] like Figure 1 As shown, this embodiment provides a surface wave dispersion probability inversion method based on approximate Bayesian calculation, which includes the following steps:

[0043] (1) The Bayesian equation for the probability inversion of surface wave dispersion is derived under the framework of approximate Bayesian computing (ABC), and the prior distribution and likelihood function are given;

[0044] In step (1), the horizontal layered model is used to approximate the shear wave velocity v of the test site. s The number of layers in the horizontal layered model is N, which needs to be specified in advance based on other site survey information. The model parameters to be inverted are the thickness H and v of each layer. s ; Use vectors H N =[H1,H2,…,H i ,…,H N-1 ]and Indicates that, where the subscript i represents the i-th layer, let Parameters are underlined to indicate vectors (dimension > 1).

[0045] The data in the Bayesian equation are the observed dispersion curves d obs =[d obs,1 ,d obs,2 ,…,d obs,M ],d obs,m Indicates frequency f m The value corresponding to the observed dispersion curve is the observed phase velocity, and M is the number of discrete points or data dimensions selected from the observed dispersion curve. d obs , model parameters i N The probability of can be expressed using the conditional probability distribution p( i N | d obs ) quantified. In the Bayesian framework, p( i N | d obs ) is called the posterior probability distribution. According to the principle of the ABC method, from the posterior distribution p( i N | d obs ) can be converted into i N and d obs The simulated sample ( y =[y1,y2,…,yM ]) of the joint posterior distribution p ε ( i N , y | d obs ), in the Bayesian framework, p( i N , y | d obs )equal:

[0046]

[0047] Where p( y | i N ) is given i N hour y The distribution of d obs It also follows the distribution; p ε ( d obs | y ) is the likelihood function, which usually has an indicator function form in ABC and its value is used to judge the observed data d obs Whether it is the same as its simulated sample y Approximately equal; P ε ( d obs )=∫p ε ( d obs | y )p( y | i N )p( i N )d( i N , y ) is the normalization constant. p( i N )yes i N The prior distribution of can be expanded into p( i N )=p( H N )p(v s | H N ), where p( H N )for H N The prior distribution of v s | H N ) is given H N hour v s The conditional prior distribution of .

[0048] (2) Determine the parameter constraints through surface roll dispersion data to specify the parameter prior distribution;

[0049] The maximum identifiable depth of the shear wave velocity interface and the minimum identifiable layer thickness are determined based on the observed surface wave dispersion data to constrain the parameter value range and specify the prior distribution.

[0050] The specific expressions of prior distribution of each parameter are:

[0051] (21) H N ) is assumed to be a distribution similar to the uniform Dirichlet distribution

[0052]

[0053] Among them, max H Yes s Maximum recognizable depth of the interface, min H is the minimum recognizable layer thickness.

[0054] (22) For p( v s | H N ), assuming that the shear wave velocity of each layer is v si Independent and truncated Gaussian distributed:

[0055]

[0056] Among them, μ si is the i-th layer v si The mean of σ i =CV×μ si Yes si The standard deviation of si The coefficient of variation of Φ is the Gaussian distribution. The corresponding cumulative probability density function; and They are v si The upper and lower bounds of possible values, Δv s The user needs to specify it according to the actual situation of the test site. si Independent of each other,

[0057]

[0058] (23) For the prior distribution p( y | i N ), usually assuming different frequencies f m The observed phase velocity d obs,m The errors ω are independent and identically distributed, and ω has a mean of zero and a variance of σ 2 Gaussian distribution. Because d obs,m =c R,m +ω m ,so d obs Or its simulated sample y The conditional distribution p( y |σ 2 , i N ) has a Gaussian distribution. Assume that σ -2 It follows the gamma distribution, Γ(σ -2 )~Gam(a,b), and integrate p( d obs |σ 2 , i N ) 2 Available

[0059]

[0060] Where n = 2a, It can be seen that the prior distribution p( y | i N ) can be assumed to be an M-dimensional student-T distribution, that is, y ~T( c R ,δ 2 I m ,n).

[0061] (3) Assume the distribution of observed data and specify the likelihood function, and eliminate the observed data error term in the likelihood function by integration;

[0062] Assuming that the observed surface wave dispersion data obeys a Gaussian distribution with a mean of 0 and an unknown variance, the likelihood function is specified as a multivariate Gaussian distribution; the unknown variance term in the likelihood function is eliminated by integration to reduce the number of inversion parameters. For the likelihood function p in formula (1.1) ε ( d obs | y ), which has the form of an indicator function

[0063] p ε ( d obs | y )=I S(ε) (1.6)

[0064] Among them, I S(ε) ( y ) is an indicator function, when y In the set S(ε) = {L( d obs - y )≤ε} S(ε) ( y )=1, then it is considered y and d obs Approximately equal, otherwise I S(ε) ( y )=0. L is calculated y and d obs A function that measures the distance or similarity between two objects, referred to as distance function.

[0065] (4) based on the specific form of the likelihood function in step (3), determining a distance function for quickly evaluating the distance between the observed data and the simulated samples;

[0066] In step (4), the distance function L( d obs - y ) is defined as

[0067]

[0068] (5) Based on the distance function in step (4), derive an analytical formula for quickly estimating the value of the distance function;

[0069] The ratio G / G′ of the dispersion curve characteristic equation G and its first-order derivative G′ with respect to the phase velocity c is used to construct a fast estimate of the observed data d obs The Euclidean distance between the simulated sample y is used to estimate the distance function value.

[0070]

[0071] The first item on the right side of the equal sign is the set Represents a set of observation data d obs Select M SUB The set of observed phase velocities under different dispersions, Δy m =d obs,m -y m ; In the second term on the right side of the equal sign, G is the characteristic equation of the dispersion curve, G′ represents the first-order derivative of G with respect to the phase velocity c, G(d obs,m -ω m ) / G′(d obs,m -ω m ) represents the value of the characteristic equation G projected on the phase velocity c axis, which is used to estimate Δy m The value of ω m is the observed phase velocity d obs,m The error, ω m Subject to T(0,δ 2 I m ,n).

[0072] (6) Based on the Bayesian equation in step (1) and the formula for analytically calculating the distance function in step (5), ABC is used in combination with subset simulation sampling to efficiently solve the Bayesian equation;

[0073] Nested Threshold Decreasing Sequences ε in Selected Subsetting Simulations 1 =+∞>ε 2 >…>ε J =ε area as the observed dispersion data d in the observation space obs The middle region corresponding to the increasingly close approximation of S(ε t )={L(dobs -y)≤ε t}; Use the ABC method to sample from these intermediate regions to efficiently solve the posterior distribution p ε (θ,y|d obs ).

[0074] In step (6), ABC-SubSim is used to solve the Bayesian equation under the ABC framework, that is, formula (1.1).

[0075] Let F represent z=( i N , y ) Prior space The area of ​​interest is

[0076] F={ z ∈Z:L( d obs - y )<ε} (1.9)

[0077] Therefore, p(z|F) is equivalent to the posterior distribution p to be solved ε ( i N , y | d obs ). Furthermore, F can be represented by a series of decreasing threshold sequences ε 1 =+∞>ε 2 >…>ε J =ε is replaced by the intermediate target area determined by where F t ={z t ∈Z t :L( d obs - y )<ε t In the ABC-SubSim method, ABC is used to evaluate a series of intermediate regions F t Corresponding distribution

[0078]

[0079] where p(z)=p( y |θ N )p(θ N ). It should be noted that the ε value should be given in advance, and other intermediate thresholds ε t It can be determined adaptively during the sampling process. Solve the posterior distribution p based on the ABC-SubSim method ε ( i N , y | d obs ) include:

[0080] (61) When t = 1, the Monte Carlo (MC) method is used to obtain p( z |F1 ) is used for sampling, and the sampling number is set to N s ; Then define the conditional probability P(F 2 |F 1 )=p0, and the N s Samples are classified according to their L( d obs - y ) values ​​in ascending order, let That is F 2 The threshold value; finally, the front p0N s samples as the next level (t = 2) from p( z |F 2 ) is the “seed” for conditional sampling in .

[0081] (62) When t ≥ 2, the Markov chain Monte Carlo (MCMC) method is used to obtain the “seed” from p( z |F t ) extract N s These samples are processed in the same way as in step S6.1 and the conditional probability P(F t+1 |F t )=p0, then the next level (t+1) target area F can be determined t+1 Threshold and p0N s A "seed" sample.

[0082] (63) Repeat step S4.2 until ε t+1 ≤ε, then let ε t+1 = ε and use MCMC method to extract the target region F t+1 = Extract N from F s samples, that is, the posterior distribution p ε ( i N , y | d obs ) as the numerical solution of formula (1-2).

[0083] By setting a suitable p0 value (usually 0.1-0.2), it is easy to obtain the t-1 Sample z t-1 Find the one that falls in F t These "seeds" are used as initial samples using the MCMC method from p(z|F t ) There is no need to go through a long “burn-in” process when sampling, which improves the sampling efficiency of the MCMC method.

[0084] (7) Quantify the uncertainty of the inversion results based on the parameter posterior samples in step (6).

[0085] Case 1: Three-layer virtual venue

[0086] Figure 2(a) shows a three-layer v s The virtual site of the section, the site material density and Poisson's ratio are uniform, that is, ρ = 2000kg / m 3 and ν = 0.3. The observed dispersion curve of the site is given by the formula d obs,m =c R,m +ω m Simulation, where c R,m The error ω is obtained by the stiffness matrix method. m Assumed to be Gaussian distribution N(0,2.5 2 ). From frequency f1 = 10 Hz to f 71 = 80 Hz, and one observed phase velocity is simulated every 1 Hz, and a total of M = 71 observed phase velocities are synthesized into a set of observed dispersion data (as shown in Figure 2(b)). Finally, assuming that the inverted v s The profile has N = 3 layers, so the model parameters to be inverted are H 3=[H1,H2] and

[0087] In order to verify the rationality and efficiency of the proposed method, in addition to using the proposed method to invert the dispersion data, two other sets of simulations were performed: 1. The distance function in the ABC method was replaced with the Euclidean distance formula to accurately evaluate the d obs and y The distance between (referred to as the accurate distance function), other parameter settings are consistent with the proposed method. At this time, each time the accurate distance function value is calculated, it is necessary to calculate the distance function from T( c R ,δ 2 I m ,n) to extract a set of observed data simulation values ​​with dimension M = 71, so it is necessary to calculate M = 71 theoretical phase velocities with different dispersions c R The results of this example are compared with the proposed method to evaluate the computational efficiency of the proposed method; 2 Based on the exact Bayesian method (that is, the traditional Bayesian method, in order to distinguish it from the approximate Bayesian method), the MCMC method is used to directly obtain the posterior distribution P( θ3 | d obs ) and the sampling results were used as reference values ​​to verify the rationality of the proposed method.

[0088] Table 1 shows the number and time of calling the dispersion curve characteristic equation G in the ABC-SusSim sampling process using the proposed distance function and the accurate distance function respectively. The computational efficiency of the proposed distance function is more than 11 times that of the accurate distance function, which greatly alleviates the time-consuming problem of evaluating the distance function under the framework of the ABC method.

[0089] Table 1 Virtual site calculation results

[0090] <![CDATA[Desired c R Dimensions]]> The number of times the characteristic equation G is called time Exact distance function 71 <![CDATA[2.35ⅹ10 8 ]]> 3718 The proposed distance function 3 <![CDATA[2.06ⅹ10 7 ]]> 277 ratio 23.67 11.40 13.40

[0091] Figure 3 The comparison of three sets of simulation results is shown. The marginal posterior histograms of each parameter obtained by the proposed method can quantify the uncertainty of the inversion results. The marginal posterior histograms are consistent with the ABC simulation results based on the accurate distance function and the accurate Bayesian results, which shows that the results of the proposed method are consistent with the posterior distribution P( θ3 | d obs ) has a better approximation effect. Figure 4(a) shows the estimated y Value, that is, y m ≈d obs,m -G(d obs,m -ω m ) / G′(d obs,m -ω m ), y Randomly distributed on both sides of the synthetic dispersion data, and d obs is very close, further emphasizing the accuracy of the proposed method. Figure 4(b) shows the v obtained by the proposed method. s Profile sample (gradient solid line) and correct v s The profile (dashed line) shows the same trend with depth, v s Profile samples are randomly distributed at the correct v s Both sides of the profile illustrate the rationality of the inversion results.

[0092] Case 2: Actual site

[0093] In order to further investigate and evaluate the applicability of the proposed method, the observed dispersion data obtained by Xia et al. (2002) (Xia, J., Miller, RD, Park, CB, Hunter, JA, Harris, JB, & Ivanov, J. (2002). Comparing shear-wave velocity profiles inverted from multichannel surface wave with borehole measurements. Soil dynamics and earthquake engineering, 22(3), 181-190.) in a MASW test in Wyoming were used. Figure 5 ) for analysis. The site material density and Poisson’s ratio are assumed to be known, i.e., 2000kg / m 3 and 0.3. Finally, assuming that the inverted v s The profile has N = 3 layers, and the model parameters to be inverted are H 3=[H1,H2] and

[0094] Table 2 shows the computational efficiency comparison results of the proposed distance function and the accurate distance function in the actual case. Because the dimension of the observed dispersion data (M=17) in the actual case used is relatively low, the computational efficiency of the proposed distance function is more than 4 times that of the accurate distance function. Compared with the results in Case 1, it can be seen that the relative efficiency of the proposed method increases with the increase of the dimension of the data.

[0095] Table 2 Actual site calculation results

[0096] <![CDATA[Dimensions required c R > The number of times the characteristic equation G is called Time(s) The proposed distance function 17 <![CDATA[5.27ⅹ10 7 ]]> 975 Exact distance function 3 <![CDATA[1.18ⅹ10 7 ]]> 201 ratio 5.67 4.48 4.83

[0097] Figure 6 The results of the proposed method in the case are compared with the other two groups of results. The three groups of results have high consistency; Figure 7(a) shows y The sample and d obs Very close. This shows that the results of the proposed method are very close to the posterior distribution P( θ3 | d obs ) has a better approximation effect. Figure 7(b) also shows the v obtained by the proposed method s The profile sample (solid line with gradient color) and the inversion results of Xia et al. (2002) (solid lines with solid dots) and v measured in hanging logging s The data (solid squares and dashed lines) show a consistent trend with depth, which shows that the proposed method has good applicability to the actual field.

[0098] In summary, the proposed method is proved to be efficient through virtual site (Case 1) and actual site (Case 2), and can reasonably quantify the uncertainty of the inversion results. Comparison of the two examples proves that the relative efficiency of the proposed method will increase with the increase of data dimension.

[0099] The above is a schematic description of the present invention and its implementation methods, which is not restrictive. The drawings show only one implementation method of the present invention, and the actual structure is not limited thereto. Therefore, if a person skilled in the art is inspired by it and does not deviate from the purpose of the invention, he or she can creatively design a structure and an embodiment similar to the technical solution, which should fall within the protection scope of the present invention.

Claims

1. A surface wave dispersion probability inversion method based on approximate Bayesian calculation, characterized by: The following steps are involved: (1) Derive the Bayesian equation for the probability inversion of surface wave dispersion under the framework of approximate Bayesian computing (ABC); (2) Determine the parameter constraints through surface roll dispersion data to specify the parameter prior distribution; (3) Assume the distribution of observed data and specify the likelihood function, and eliminate the observed data error term in the likelihood function by integration; (4) based on the specific form of the likelihood function in step (3), determining a distance function for evaluating the distance between the observed data and the simulated samples; (5) Based on the distance function in step (4), derive an analytical formula for quickly estimating the value of the distance function; (6) Based on the Bayesian equation of step (1) and the formula for analytically calculating the distance function of step (5), ABC is used in combination with subset simulation sampling to solve the Bayesian equation; (7) Quantify the uncertainty of the inversion results based on the parameter posterior samples in step (6).

2. The method for inverting surface wave dispersion probability based on approximate Bayesian calculation according to claim 1, characterized in that: In step (1), the horizontal layered model is used to approximate the shear wave velocity v of the test site. s The number of layers in the horizontal layered model is N, and the model parameters to be inverted are the thickness H and v of each layer. s ; Given the observed dispersion curve d obs Data, the probability of θ in the ABC framework is expressed by its posterior distribution p(θ|d obs ) quantized, solving the posterior distribution p(θ|d obs ) can be converted into solving the joint posterior distribution of θ and the simulated sample y of the observed dispersion curve data, that is: where p(y|θ) is the prior distribution of y given θ, i.e., d obs The distribution it obeys; p(θ) is the prior distribution of θ; p ε (d obs |y) is the likelihood function; P ε (d obs ) is the normalization constant; in ABC, the likelihood function p ε (d obs |y) has an indicator function I S(ε) In the form of, when y is in the set S(ε)={L(d obs -y)≤ε} S(ε) (y)=1, then it is considered that y and d obs Approximately equal; otherwise, y and d obs Not equal; L is the calculation of y and d obs The distance function of the distance between .

3. The method for inverting surface wave dispersion probability based on approximate Bayesian calculation according to claim 1, characterized in that: In step (2), the maximum identifiable depth of the shear wave velocity interface and the minimum identifiable layer thickness are determined based on the observed surface wave dispersion data to constrain the parameter value range and specify the prior distribution.

4. The method for inverting surface wave dispersion probability based on approximate Bayesian calculation according to claim 1, characterized in that: In step (3), it is assumed that the observed surface wave dispersion data obeys a Gaussian distribution with a mean of 0 and an unknown variance, so the likelihood function is specified as a multivariate Gaussian distribution; the unknown variance term in the likelihood function is eliminated by the integration method to reduce the number of inversion parameters.

5. The method for inverting surface wave dispersion probability based on approximate Bayesian calculation according to claim 1, characterized in that: In step (4), the Euclidean distance is used to quantitatively evaluate the surface wave dispersion observation data d obs The distance between the simulated sample y.

6. The method for inverting surface wave dispersion probability based on approximate Bayesian calculation according to claim 1, characterized in that: In step (5), the ratio G / G′ of the dispersion curve characteristic equation G and its first-order derivative G′ with respect to the phase velocity c is used to construct a fast estimation of the observed data d obs The Euclidean distance between the simulated sample y is used to estimate the distance function value.

7. The method for inverting surface wave dispersion probability based on approximate Bayesian calculation according to claim 1, characterized in that: In step (6), the nested threshold decreasing sequence ε in the subset simulation is selected 1 =+∞>ε 2 >…>ε J =ε area as the observed dispersion data d in the observation space obs The middle region corresponding to the increasingly close approximation of S(ε t )={L(d obs -y)≤ε t }; Use the ABC method to sample from these intermediate regions to solve the posterior distribution p ε (θ,y|d obs ).

Images