A method and system for inversion of Rayleigh wave dispersion curve based on geostatistics

Abstract

The present invention relates to a Rayleigh wave dispersion curve inversion method and system based on geological statistics, and belongs to the technical field of seismic shear wave velocity inversion. First, an inversion equation is constructed based on a dispersion equation, and iterative optimization is performed using the inversion equation to obtain the first seismic shear wave velocity structure of each station. Then, multiple second initial seismic shear wave velocity structures of each inter-station point to be estimated are determined based on the first seismic shear wave velocity structure and position of all stations, and the second Rayleigh wave dispersion curve is calculated using the dispersion equation. Finally, all the second Rayleigh wave dispersion curves corresponding to the inter-station point to be estimated are used as input, and the optimal second Rayleigh wave dispersion curve is determined using the second objective function to determine the second seismic shear wave velocity structure of each inter-station point to be estimated. The present invention can simultaneously obtain the seismic shear wave velocity structure of stations and stations, and can greatly improve the lateral resolution.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic shear wave velocity inversion, and in particular to a Rayleigh wave dispersion curve inversion method and system based on geological statistics. Background Art

[0002] At present, the seismic shear wave velocity structure is mostly obtained by inverting the Rayleigh wave dispersion curve. However, this method can only obtain the seismic shear wave velocity structure of a station, but cannot obtain the seismic shear wave velocity structure between stations. In engineering applications, there is a problem of low lateral resolution of the inversion results.

[0003] Based on this, an inversion technology that can improve the lateral resolution is urgently needed. Summary of the invention

[0004] The purpose of the present invention is to provide a method and system for inverting Rayleigh wave dispersion curves based on geostatistics, which can simultaneously obtain the seismic shear wave velocity structure of stations and between stations, and can greatly improve the lateral resolution.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A Rayleigh wave dispersion curve inversion method based on geostatistics, the inversion method comprising:

[0007] For each station, a first initial seismic shear wave velocity structure and a first Rayleigh wave dispersion curve of the station are obtained; an inversion equation is constructed based on a dispersion equation; the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve are used as inputs, and the inversion equation is used for iterative optimization to obtain the first seismic shear wave velocity structure of the station; the Rayleigh wave dispersion curve is a curve showing the change of Rayleigh wave phase velocity with frequency;

[0008] Determine multiple second initial seismic shear wave velocity structures of each inter-station point to be estimated according to the first seismic shear wave velocity structures and positions of all the stations; the inter-station point to be estimated is a grid node between two adjacent stations;

[0009] For each of the second initial seismic shear wave velocity structures of each of the inter-station estimated points, taking the second initial seismic shear wave velocity structure as input, and using the dispersion equation to calculate a second Rayleigh wave dispersion curve;

[0010] For each of the inter-station points to be estimated, all the second Rayleigh wave dispersion curves corresponding to the inter-station points to be estimated are taken as input, and the optimal second Rayleigh wave dispersion curve is determined using the second objective function. The second initial seismic shear wave velocity structure corresponding to the optimal second Rayleigh wave dispersion curve is the second seismic shear wave velocity structure of the inter-station points to be estimated.

[0011] A Rayleigh wave dispersion curve inversion system based on geostatistics, the inversion system comprising:

[0012] The station velocity inversion module is used for obtaining, for each station, a first initial seismic shear wave velocity structure and a first Rayleigh wave dispersion curve of the station; constructing an inversion equation based on a dispersion equation; taking the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve as input, using the inversion equation for iterative optimization to obtain the first seismic shear wave velocity structure of the station; the Rayleigh wave dispersion curve is a curve of the change of Rayleigh wave phase velocity with frequency;

[0013] An inter-station velocity initial determination module is used to determine a plurality of second initial seismic shear wave velocity structures of each inter-station point to be estimated according to the first seismic shear wave velocity structures and positions of all the stations; the inter-station point to be estimated is a grid node between two adjacent stations;

[0014] A forward modeling module, for obtaining a second Rayleigh wave dispersion curve by using the dispersion equation to calculate each of the second initial seismic shear wave velocity structures of each of the inter-station estimated points, taking the second initial seismic shear wave velocity structures as input;

[0015] The inter-station velocity estimation module is used to determine the optimal second Rayleigh wave dispersion curve for each inter-station point to be estimated, using all the second Rayleigh wave dispersion curves corresponding to the inter-station point to be estimated as input, and using a second objective function to determine the optimal second Rayleigh wave dispersion curve, wherein the second initial seismic shear wave velocity structure corresponding to the optimal second Rayleigh wave dispersion curve is the second seismic shear wave velocity structure of the inter-station point to be estimated.

[0016] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0017] The present invention is used to provide a method and system for inverting Rayleigh wave dispersion curves based on geological statistics. First, for each station, an inversion equation is constructed based on a dispersion equation. The first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve of the station are used as inputs, and the inversion equation is used for iterative optimization to obtain the first seismic shear wave velocity structure of the station. Then, multiple second initial seismic shear wave velocity structures of each inter-station to-be-estimated point are determined according to the first seismic shear wave velocity structures and positions of all stations. The second initial seismic shear wave velocity structure is used as input, and the dispersion equation is used to calculate the second Rayleigh wave dispersion curve. All second Rayleigh wave dispersion curves corresponding to the inter-station to-be-estimated point are used as inputs, and the second objective function is used to determine the optimal second Rayleigh wave dispersion curve. The second initial seismic shear wave velocity structure corresponding to the optimal second Rayleigh wave dispersion curve is the second seismic shear wave velocity structure of the inter-station to-be-estimated point, thereby obtaining the second seismic shear wave velocity structure of each inter-station to-be-estimated point. The present invention can simultaneously obtain the seismic shear wave velocity structure of the station and the inter-station, and can greatly improve the lateral resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0019] Figure 1 A method flow chart of the inversion method provided in Example 1 of the present invention;

[0020] Figure 2 This is a system block diagram of the inversion system provided in Example 2 of the present invention. DETAILED DESCRIPTION

[0021] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0022] The purpose of the present invention is to provide a Rayleigh wave dispersion curve inversion method and system based on geological statistics, which can simultaneously obtain the seismic shear wave velocity structure of stations and between stations, greatly improve the lateral resolution, and obtain the seismic shear wave velocity structure with higher lateral resolution.

[0023] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0024] Embodiment 1:

[0025] This embodiment is used to provide a Rayleigh wave dispersion curve inversion method based on geostatistics, such as Figure 1 As shown, the inversion method includes:

[0026] S1: For each station, obtain the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve of the station; construct an inversion equation based on the dispersion equation; take the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve as input, use the inversion equation to perform iterative optimization, and obtain the first seismic shear wave velocity structure of the station; the Rayleigh wave dispersion curve is a curve of the change of Rayleigh wave phase velocity with frequency;

[0027] The station in this embodiment is a seismic station, which refers to a basic seismic observation institution having at least one seismic observation site, observation facilities and implementation management functions.

[0028] In S1, obtaining the first initial seismic shear wave velocity structure of the station may include: first obtaining the rock type (or lithology) of each layer of the station, and then for each layer, determining the initial velocity of the layer according to the rock type of the layer, specifically, the rock type corresponds to the velocity range, such as the shear wave velocity range of sedimentary basalt is 3050-4500m / s, randomly selecting a velocity from the velocity range as the initial velocity of the layer, and the initial velocities of all layers constitute the first initial seismic shear wave velocity structure, so as to roughly estimate the first initial seismic shear wave velocity structure of the station according to the rock type of each layer of the station, and the first initial seismic shear wave velocity structure includes the initial velocity of each layer. Since the rock type and the seismic shear wave velocity structure only roughly correspond, it is still necessary to determine the accurate first seismic shear wave velocity structure of the station through subsequent inversion, and the first initial seismic shear wave velocity structure will be used as an initial model for subsequent inversion.

[0029] In S1, obtaining the first Rayleigh wave dispersion curve of the station may include: collecting and acquiring the time domain seismic waveform data of the station; extracting the Rayleigh wave components of the time domain seismic waveform data, generally, the signals on both sides of the seismic trace except for the obvious reflection signals are regarded as Rayleigh waves, so the Rayleigh wave components can be obtained by directly cutting out this part of the signal, and the Rayleigh wave components can be understood as the curve of the change of Rayleigh wave phase velocity with time; using the cross-correlation method to process the Rayleigh wave components, and estimate the first Rayleigh wave dispersion curve, and the Rayleigh wave dispersion curve is the curve of the change of Rayleigh wave phase velocity with frequency.

[0030] Before constructing the inversion equation based on the dispersion equation, this embodiment first determines the form of the inversion solution of the Rayleigh wave dispersion curve.

[0031] The dispersion equation generally includes the frequency f j , Rayleigh wave phase velocity C Rj , seismic shear wave velocity V s , seismic longitudinal wave velocity V p , density ρ and thickness d, which are in the following form:

[0032] F(f j ,C Rj ,V s ,V p ,ρ,d)=0 (j=1,2,…,m);

[0033] Where m represents the number of frequencies; C Rj is a real number, V s ,V p,ρ are N-dimensional vectors, d is an N-1-dimensional vector, and N is the number of layers. In the subsequent inversion process, when the Rayleigh wave dispersion curve is known, the parameter V in the dispersion equation can be solved s ,V p ,ρ and d. Usually, the Rayleigh wave phase velocity is only sensitive to the seismic shear wave velocity, and is not sensitive to the seismic longitudinal wave velocity, density and thickness. Therefore, the seismic longitudinal wave velocity, density and thickness can all be used as known quantities. Furthermore, the dispersion equation can be regarded as j The Rayleigh wave phase velocity and seismic shear wave velocity V are described below: s An implicit function of , solving m dispersion equations, can obtain the Rayleigh wave dispersion curve corresponding to a certain seismic shear wave velocity structure. Therefore, the solution of the Rayleigh wave dispersion curve inversion is the seismic shear wave velocity structure. Based on this, this embodiment uses the roughly estimated first initial seismic shear wave velocity structure as the initial model for inversion, performs inversion processing on the initial model, and finally obtains a more accurate first seismic shear wave velocity structure.

[0034] This embodiment can collect relevant geological data of the work area to determine the seismic longitudinal wave velocity, density and thickness, so as to collect prior geological information in advance. Among them, the thickness of each layer is determined as d = (d1, d2, ..., d N-1 ), N represents the number of layers, and the last layer is assumed to be a uniform half-space with infinite thickness.

[0035] In this embodiment, constructing the inversion equation based on the dispersion equation may include:

[0036] (1) Perform Taylor expansion on the dispersion equation to obtain the first function;

[0037] Assume that the seismic shear wave velocity structure of the underground medium at the seismic station location is V s =(v s1 ,v s2 ,…,v sN ) T , the Rayleigh wave phase velocity is C R =(C R1 ,C R2 ,…,C Rm ) T , Taylor expansion is performed on the dispersion equation, and the first function obtained is:

[0038] JΔV s =ΔC R ;

[0039] Where J is the Jacobian matrix; ΔV s Indicates a small disturbance in the velocity structure of the first earthquake shear wave. In the subsequent iterative optimization process of inversion, this small disturbance represents the difference in the velocity structure of the first earthquake shear wave between two adjacent iterations; ΔC RIt represents a small perturbation of the Rayleigh wave phase velocity. In the subsequent iterative optimization process of inversion, this small perturbation represents the difference between the Rayleigh wave phase velocities in two adjacent iterations.

[0040] (2) constructing a first objective function based on the first function;

[0041] In this embodiment, the first objective function is as follows:

[0042]

[0043] Where Φ is the first objective function; J is the Jacobian matrix; ΔV s is the difference between the first seismic shear wave velocity structure in two consecutive iterations; ΔC R is the difference in the Rayleigh wave phase velocity between two adjacent iterations; W is a positive definite diagonal weight matrix used to balance the residual effects of different frequencies. W can be a known matrix given artificially or a covariance matrix of data errors. W = L T L, where L is a diagonal matrix and λ is the regularization parameter.

[0044] (3) Minimize the first objective function to obtain the inversion equation.

[0045] The first objective function is used to evaluate the inversion results and solve the drop. The first seismic shear wave velocity structure obtained by inversion should make the objective function value as small as possible. Minimize Φ to obtain the inversion equation. Minimization is to instruct the first objective function to minimize ΔV s The process of finding the minimum value of a function by taking the partial derivative and setting it to 0. The inverse equation is used to describe the amount of descent obtained by solving the descent direction of the objective function, as follows:

[0046] ΔV s =V(Λ 2 +λI) -1 ΛU T d;

[0047] Among them, A=LJ is subjected to singular value decomposition to obtain A=UΛV T , U, Λ, V are the left eigenvector matrix, diagonal matrix and right eigenvector matrix respectively; I is the unit matrix; d = LΔC R represents the weighted phase velocity vector.

[0048] Based on the above inversion equation, the following iterative format can be constructed:

[0049] (ΔV s ) n+1 =V(Λ 2 +λI) -1 ΛU T (d) n ;

[0050] (V s ) n+1 =(V s ) n +(ΔV s ) n+1 ;

[0051] Where (ΔV s ) n+1 represents the increment of seismic shear wave velocity structure in the n+1th iteration; (d) n represents the weighted phase velocity vector in the nth iteration; (V s ) n+1 and (V s ) n They represent the seismic shear wave velocity structures of the n+1th and nth iterations respectively.

[0052] This embodiment uses the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve as input, and uses the above inversion equation to perform continuous iterative optimization until (ΔV s ) n+1 When the modulus is less than a certain threshold, such as 0.001, the loop is exited to obtain the first earthquake shear wave velocity structure of the seismic station.

[0053] Since the first objective function contains a regularization parameter λ, this embodiment can use the L curve method to quantify the relationship between the regularization parameter and the first objective function, by setting the value range of the regularization parameter λ to [100, 0.01], and determining its value by the L curve method, and then outputting the inversion result. Specifically, a plurality of values ​​are randomly selected within the value range of the regularization parameter λ; for each value, the above-mentioned inversion process is performed once to obtain the seismic shear wave velocity structure, and then the seismic shear wave velocity structure is brought into the first objective function to obtain the first function value (i.e., the fitting difference before the plus sign) and the second function value (i.e., the roughness value after the plus sign) of the first objective function, thereby obtaining the fitting difference and roughness value corresponding to each value; with the fitting difference as the ordinate, the roughness value as the abscissa, and the fitting difference and roughness value corresponding to all values ​​as data points, L curve fitting is performed, and the regularization parameter at the inflection point is used as the optimal regularization parameter. The seismic shear wave velocity structure corresponding to the optimal regularization parameter is the first seismic shear wave velocity structure of the station output by S1.

[0054] Through S1, the inversion process of the seismic shear wave velocity structure of the underground medium at the station location can be completed, and the first seismic shear wave velocity structure of each station can be obtained.

[0055] This embodiment uses a geostatistical method to estimate the seismic shear wave velocity structure between stations. The geostatistical method is used to estimate the seismic shear wave velocity structure of the underground medium between two seismic stations, and obtain samples that have both spatial correlation and diversity. By optimizing the samples, an inversion result with a higher lateral resolution can be obtained, see subsequent S2-S4.

[0056] S2: determining a plurality of second initial seismic shear wave velocity structures of each inter-station point to be estimated according to the first seismic shear wave velocity structures and positions of all the stations; the inter-station point to be estimated is a grid node between two adjacent stations;

[0057] In this embodiment, the underground medium of the working area is discretized to obtain a plurality of grid nodes. The working area includes each station in S1, and the seismic shear wave velocity structure of the grid node corresponding to the station is calculated by S1. The grid nodes between stations are the inter-station estimated points, and their seismic shear wave velocity structure is calculated by S2-S4.

[0058] The first seismic shear wave velocity structures of the underground media of all seismic stations that have been obtained are placed at their corresponding grid nodes, and the velocity values ​​on the grid nodes are used as known conditions to perform random simulations on the second initial seismic shear wave velocity structures of the inter-station points to be estimated, and obtain multiple second initial seismic shear wave velocity structures of each inter-station point to be estimated. Specifically, the random simulation process may include:

[0059] (1) Randomly define the estimation order of multiple groups of all inter-station estimated points to obtain multiple random simulation orders;

[0060] (2) For each random simulation order, the variogram is obtained based on the first seismic shear wave velocity structure and position fitting of all stations, and the second initial seismic shear wave velocity structure of the first inter-station estimated point is determined using the variogram;

[0061] (3) A new variogram is obtained by fitting the first seismic shear wave velocity structure and position of all stations and the second initial seismic shear wave velocity structure and position of all inter-station estimated points with known velocities, and the second initial seismic shear wave velocity structure of the next inter-station estimated point is determined using the new variogram;

[0062] Before interpolation, it is necessary to establish a variogram to describe the statistical characteristics of the seismic shear wave velocity structure as it changes with space. The variogram value in geostatistics is obtained by the following formula:

[0063]

[0064] Where h is the variable range, M represents the total number of nodes with known velocities, including stations and inter-station estimated points; Z(x i ) is the i-th node x whose speed is knowni The value of the seismic shear wave velocity structure; Z(x i +h) is the same as x i The value of the seismic shear wave velocity structure of the node at a distance h. The above function v describes a quantitative relationship, and its continuous form needs to be fitted.

[0065] Using the above quantitative relationship, nonlinear fitting is performed on the distance of nodes with known speed (the distance is calculated based on the position of the node) and the speed, and a variogram that can take values ​​for any distance is obtained by fitting. Specifically, the value of v(h) is first obtained based on the nodes with known speed, and the base value, range, lag distance and other parameters are fitted with multiple values ​​of h, so that the function v(h) corresponding to any h can be obtained. The type of the variogram of this embodiment is a spherical model, an exponential function model or a power function model. Taking the spherical model as an example, the variogram can be expressed as:

[0066]

[0067] Among them, c represents the base value; h represents the hysteresis distance, that is, the distance; a represents the range.

[0068] Determining the second initial seismic shear wave velocity structure of the inter-station to-be-estimated point using the new variogram may include: estimating the mean and variance of the inter-station to-be-estimated point using the seismic shear wave velocity structure and variogram of the node with known velocity, the mean and variance are actually the Kriging mean and Kriging variance, that is, in the process of Kriging interpolation, the Kriging equation group is established through constraints such as variance minimization, and the weight coefficient required for interpolation is solved to further obtain the Kriging interpolation of the position to be estimated, and this is used as the Kriging mean, and the Kriging variance is then calculated by the Kriging mean. The local probability distribution of the inter-station to-be-estimated point can be determined based on the mean and variance. In this embodiment, it can be assumed that the seismic shear wave velocity structure of the inter-station to-be-estimated point satisfies the normal distribution, and the local probability distribution can be determined after determining the mean and variance. A random sample is extracted from the local probability distribution, and the value of the random sample is the second initial seismic shear wave velocity structure of the inter-station to-be-estimated point. It is recorded as a node with known speed, and the variogram is re-estimated. After the valuation of the inter-station points to be estimated is completed, it is equivalent to adding samples. The variogram can be re-estimated and the calculation of the next inter-station point to be estimated is carried out. The above process is repeated until the calculation of all inter-station points to be estimated is completed.

[0069] (4) Determine whether the second initial seismic shear wave velocity structure of all inter-station estimated points is obtained;

[0070] (5) If yes, then the iteration ends; the second initial seismic shear wave velocity structure of each inter-station to-be-estimated point under all random simulation orders constitutes a plurality of second initial seismic shear wave velocity structures of each inter-station to-be-estimated point;

[0071] (6) If not, return to the step of "obtaining a new variogram based on the first seismic shear wave velocity structure and position of all stations and the second initial seismic shear wave velocity structure and position of all inter-station estimated points with known velocities".

[0072] Based on the above process, this embodiment can generate multiple groups of second initial seismic shear wave velocity structures of inter-station estimated points under the same grid, and combine the first seismic shear wave velocity structures of the stations to form multiple groups of random realization models, which can be understood as multiple seismic shear wave velocity structure samples obtained in the same working area. In S3, these samples will be forward modeled, and in S4, the correlation between the forward modeling results of these samples and the known forward modeling results will be used to select the best matching random realization model as the final seismic shear wave velocity structure.

[0073] S3: for each of the second initial seismic shear wave velocity structures of each of the inter-station estimated points, using the second initial seismic shear wave velocity structure as input, and using the dispersion equation to calculate a second Rayleigh wave dispersion curve;

[0074] This step is the forward modeling process of multiple sets of random realization models. For each inter-station estimated point, it has multiple sets of random realization models. For multiple sets of random realization models, the dispersion equation is used to generate multiple sets of Rayleigh wave dispersion curves corresponding to them. Based on the multiple sets of Rayleigh wave dispersion curves, the (i.e. phase velocity), where j represents the jth group of random realization models and i represents the i-th inter-station point to be estimated.

[0075] S4: For each of the inter-station points to be estimated, all the second Rayleigh wave dispersion curves corresponding to the inter-station points to be estimated are taken as input, and the optimal second Rayleigh wave dispersion curve is determined using the second objective function. The second initial seismic shear wave velocity structure corresponding to the optimal second Rayleigh wave dispersion curve is the second seismic shear wave velocity structure of the inter-station points to be estimated.

[0076] The second objective function of this embodiment is:

[0077]

[0078] Among them, Φ i is the second objective function; PV i is the estimated value of the Rayleigh wave phase velocity of the i-th inter-station point to be estimated; COV represents the correlation coefficient; is the Rayleigh wave phase velocity of the jth second Rayleigh wave dispersion curve of the i-th inter-station estimated point.

[0079] The calculation formula for the estimated value of the Rayleigh wave phase velocity of the i-th inter-station estimated point is:

[0080]

[0081] Where, s = 1, 2, ..., N1, N1 is the N1 nodes closest to the i-th inter-station point to be estimated, and the size of N1 can be set by yourself; ω s is the inverse of the distance from the sth node to the i-th station to be estimated; PV s is the Rayleigh wave phase velocity of the sth node; nodes include stations and inter-station points to be estimated.

[0082] Through the second objective function, we can select i The Rayleigh wave phase velocity of the i-th inter-station to-be-estimated point with the largest correlation coefficient is taken as the optimal phase velocity, and the second initial seismic shear wave velocity structure corresponding to the optimal phase velocity is the optimal matching result of the inter-station to-be-estimated point, that is, the second seismic shear wave velocity structure of the inter-station to-be-estimated point.

[0083] This embodiment is different from the traditional technology of determining the shear wave velocity between seismic stations through linear interpolation method. The random simulation method adopted by S2 can increase the randomness and diversity of data under the premise of fully considering the spatial statistical correlation. S3 and S4 provide a preferred basis for a large number of random realization models, thereby ensuring that more reliable results are selected from the random realization models in S2, so as to obtain a seismic shear wave velocity structure with higher lateral resolution.

[0084] Embodiment 2:

[0085] This embodiment is used to provide a Rayleigh wave dispersion curve inversion system based on geostatistics, such as Figure 2 As shown, the inversion system comprises:

[0086] The station velocity inversion module M1 is used for obtaining the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve of each station; constructing an inversion equation based on the dispersion equation; taking the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve as input, using the inversion equation to perform iterative optimization to obtain the first seismic shear wave velocity structure of the station; the Rayleigh wave dispersion curve is a curve of the change of Rayleigh wave phase velocity with frequency;

[0087] The inter-station velocity initial determination module M2 is used to determine a plurality of second initial seismic shear wave velocity structures of each inter-station point to be estimated according to the first seismic shear wave velocity structures and positions of all the stations; the inter-station point to be estimated is a grid node between two adjacent stations;

[0088] A forward modeling module M3 is used for obtaining a second Rayleigh wave dispersion curve by using the dispersion equation to calculate each of the second initial seismic shear wave velocity structures of each of the inter-station estimated points, taking the second initial seismic shear wave velocity structures as input;

[0089] The inter-station velocity estimation module M4 is used to determine the optimal second Rayleigh wave dispersion curve for each inter-station point to be estimated by using the second objective function, taking all the second Rayleigh wave dispersion curves corresponding to the inter-station point to be estimated as input, and the second initial seismic shear wave velocity structure corresponding to the optimal second Rayleigh wave dispersion curve is the second seismic shear wave velocity structure of the inter-station point to be estimated.

[0090] Each embodiment in this specification focuses on the differences from other embodiments, and the same or similar parts between the embodiments can be referred to each other. For the system disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0091] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for inverting Rayleigh wave dispersion curves based on geostatistics, characterized in that: The inversion method includes: For each station, a first initial seismic shear wave velocity structure and a first Rayleigh wave dispersion curve of the station are obtained; an inversion equation is constructed based on the dispersion equation; the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve are used as input, and an iterative optimization is performed using the inversion equation to obtain the first seismic shear wave velocity structure of the station; the Rayleigh wave dispersion curve is a curve showing the change of Rayleigh wave phase velocity with frequency; Determine multiple second initial seismic shear wave velocity structures of each inter-station point to be estimated based on the first seismic shear wave velocity structures and positions of all the stations; the inter-station point to be estimated is a grid node between two adjacent stations; For each of the second initial seismic shear wave velocity structures of each inter-station point to be estimated, using the second initial seismic shear wave velocity structure as input, the dispersion equation is used to calculate a second Rayleigh wave dispersion curve; For each of the inter-station to-be-estimated points, all the second Rayleigh wave dispersion curves corresponding to the inter-station to-be-estimated points are used as input, and an optimal second Rayleigh wave dispersion curve is determined using a second objective function, wherein the second initial seismic shear wave velocity structure corresponding to the optimal second Rayleigh wave dispersion curve is the second seismic shear wave velocity structure of the inter-station to-be-estimated point; The second objective function is: Among them, Φ i is the second objective function; PV i is the estimated value of the Rayleigh wave phase velocity of the i-th inter-station point to be estimated; is the Rayleigh wave phase velocity of the jth second Rayleigh wave dispersion curve of the i-th inter-station estimated point; The calculation formula for the estimated value of the Rayleigh wave phase velocity of the i-th inter-station estimated point is: Where s = 1, 2, ..., N1, N1 is the N1 nodes closest to the i-th inter-station point to be estimated; ω s The reciprocal of the distance from the sth node to the i-th station to be estimated; PV s is the Rayleigh wave phase velocity of the sth node; the nodes include stations and inter-station points to be estimated.

2. The inversion method according to claim 1, characterized in that: Obtaining the first initial seismic shear wave velocity structure of the station specifically includes: Obtain the rock type of each layer at the station; For each of the horizons, the initial velocity of the horizon is determined according to the rock type of the horizon; the initial velocities of all the horizons constitute a first initial seismic shear wave velocity structure.

3. The inversion method according to claim 1, characterized in that: Obtaining the first Rayleigh wave dispersion curve of the station specifically includes: Acquiring time-domain seismic waveform data of the station; extracting Rayleigh wave components of the time-domain seismic waveform data; The Rayleigh wave component is processed using a cross-correlation method to obtain a first Rayleigh wave dispersion curve.

4. The inversion method according to claim 1, characterized in that: The inversion equation constructed based on the dispersion equation specifically includes: Taylor expansion is performed on the dispersion equation to obtain the first function; A first objective function is constructed based on the first function; the first objective function is Where Φ is the first objective function; J is the Jacobian matrix; ΔV s is the difference in the shear wave velocity structure of the first earthquake between two adjacent iterations; ΔC R is the difference in Rayleigh wave phase velocity between two adjacent iterations; W is the positive definite diagonal weight matrix; λ is the regularization parameter; The first objective function is minimized to obtain an inversion equation.

5. The inversion method according to claim 4, characterized in that: The regularization parameter is calculated using the L-curve method.

6. The inversion method according to claim 1, characterized in that: The step of determining a plurality of second initial seismic shear wave velocity structures of each inter-station estimated point based on the first seismic shear wave velocity structures and positions of all the stations specifically includes: Randomly define the estimation order of all inter-station estimated points in multiple groups to obtain multiple random simulation orders; For each of the random simulation sequences, a variogram is obtained based on the first seismic shear wave velocity structures and positions of all the stations, and the second initial seismic shear wave velocity structure of the first inter-station point to be estimated is determined using the variogram; A new variogram is obtained by fitting the first seismic shear wave velocity structure and position of all the stations and the second initial seismic shear wave velocity structure and position of all the inter-station estimated points with known velocities, and the second initial seismic shear wave velocity structure of the next inter-station estimated point is determined using the new variogram; determining whether the second initial seismic shear wave velocity structure of all the inter-station estimated points is obtained; If yes, then the iteration ends; the second initial seismic shear wave velocity structures of each of the inter-station to-be-estimated points under all the random simulation orders constitute multiple second initial seismic shear wave velocity structures of each of the inter-station to-be-estimated points; If not, return to the step of "obtaining a new variogram by fitting the first seismic shear wave velocity structure and position of all the stations and the second initial seismic shear wave velocity structure and position of all the inter-station estimated points with known velocities".

7. The inversion method according to claim 6, characterized in that: The type of the variogram is a spherical model, an exponential function model or a power function model.

8. A Rayleigh wave dispersion curve inversion system based on geostatistics, characterized by: The inversion system comprises: A station velocity inversion module is configured to obtain, for each station, a first initial seismic shear wave velocity structure and a first Rayleigh wave dispersion curve of the station; construct an inversion equation based on a dispersion equation; and use the first initial seismic shear wave velocity structure and the first Rayleigh wave dispersion curve as inputs to iteratively optimize the inversion equation to obtain the first seismic shear wave velocity structure of the station; the Rayleigh wave dispersion curve is a curve showing the variation of the Rayleigh wave phase velocity with frequency; an inter-station velocity initial determination module, configured to determine a plurality of second initial seismic shear wave velocity structures of each inter-station point to be estimated based on the first seismic shear wave velocity structures and positions of all the stations; the inter-station point to be estimated is a grid node between two adjacent stations; a forward modeling module configured to calculate, for each of the second initial seismic shear wave velocity structures of each of the inter-station estimated points, a second Rayleigh wave dispersion curve using the dispersion equation with the second initial seismic shear wave velocity structure as input; an inter-station velocity estimation module, configured to determine, for each inter-station to-be-estimated point, an optimal second Rayleigh wave dispersion curve using all second Rayleigh wave dispersion curves corresponding to the inter-station to-be-estimated point as input, and to use a second objective function to determine the optimal second Rayleigh wave dispersion curve, wherein the second initial seismic shear wave velocity structure corresponding to the optimal second Rayleigh wave dispersion curve is the second seismic shear wave velocity structure of the inter-station to-be-estimated point; The second objective function is: Among them, Φ i is the second objective function; PV i is the estimated value of the Rayleigh wave phase velocity of the i-th inter-station point to be estimated; is the Rayleigh wave phase velocity of the jth second Rayleigh wave dispersion curve of the i-th inter-station estimated point; The calculation formula for the estimated value of the Rayleigh wave phase velocity of the i-th inter-station estimated point is: Where s = 1, 2, ..., N1, N1 is the N1 nodes closest to the i-th inter-station point to be estimated; ω s The reciprocal of the distance from the sth node to the i-th station to be estimated; PV s is the Rayleigh wave phase velocity of the sth node; the nodes include stations and inter-station points to be estimated.

Images