Multi-mode high-precision surface wave dispersion inversion method based on shallow optimization

Abstract

The present application belongs to the technical field of surface wave dispersion inversion, and particularly provides a multi-mode high-precision surface wave dispersion inversion method based on shallow layer optimization, and the specific steps comprise: step one, establishing a geophysical model; a geophysical model containing three layers of medium is established; step two, forward calculation; an improved global matrix method is used to solve Rayleigh wave dispersion equation; step three, adaptive phase velocity search; according to frequency V, the search range is dynamically adjusted to V min ~ V max : step four, constructing a shallow layer optimization objective function; step five, simulated annealing global optimization. The present application can significantly improve the accuracy of shallow layer parameters, reduce the first layer thickness error from 5-8% to 1-2%, and reduce the first layer velocity error from 6-10% to 1-3%; the calculation efficiency is greatly improved, the numerical stability is enhanced, the multi-mode data is fully utilized, and the result is complete and clear.

Description

Technical Field

[0001] This invention belongs to the field of surface wave dispersion inversion technology, and specifically provides a multi-mode high-precision surface wave dispersion inversion method based on shallow optimization. Background Technology

[0002] Surface wave dispersion curve inversion is an important method for obtaining subsurface shear wave velocity structure, and it is widely used in engineering geological exploration, seismic safety assessment, resource exploration, and other fields. Traditional inversion methods have the following problems:

[0003] Insufficient accuracy of shallow parameter inversion: Conventional methods have large inversion errors for shallow geological parameters, which makes it difficult to meet the needs of high-precision engineering exploration.

[0004] Low computational efficiency: Traditional optimization algorithms have slow convergence speed and long computation time; solving the dispersion equation in forward modeling involves complex matrix operations, which is computationally intensive and difficult to meet the needs of efficient inversion in practical engineering.

[0005] Insufficient utilization of multimodal data: Failure to fully utilize the constraint capabilities of high-order modal data on deep structures;

[0006] Poor numerical stability: Numerical instability is prone to occur when calculating dispersion equations.

[0007] Insufficient shallow resolution: High-frequency surface waves are sensitive to shallow structures, but traditional inversion methods often lack sufficient constraints on shallow parameters, resulting in low accuracy in shallow thickness and velocity inversion.

[0008] Local extremum problem: Conventional local optimization algorithms (such as the least squares method) are highly dependent on the initial model and are prone to getting trapped in local extrema, which affects the stability of the inversion.

[0009] In summary, targeted optimizations are needed for the surface dispersion curve inversion method. Summary of the Invention

[0010] To solve the above-mentioned technical problems, the technical solution adopted by this invention is: a multi-mode high-precision surface wave dispersion inversion method based on shallow optimization, the specific steps of which include:

[0011] Step 1: Establish a geophysical model;

[0012] Establish a geophysical model that includes three layers of media;

[0013] Step 2: Forward modeling;

[0014] An improved global matrix method is used to solve the Rayleigh wave dispersion equation;

[0015] Step 3: Adaptive phase velocity search;

[0016] The search range is dynamically adjusted to V based on the frequency V. min ~V max :

[0017] V min (f)=max(100,min(V s )·γ1(f))

[0018] V max (f)=max(V s )·γ2(f)

[0019] γ is the scaling factor;

[0020] Step 4: Construct a shallow optimization objective function;

[0021] Step 5: Simulate annealing for global optimization;

[0022] An improved simulated annealing algorithm is used for inversion, and the temperature T update function is:

[0023] T k+1 =α·T k;

[0024] In the formula,

[0025] k is the number of iterations;

[0026] α is a constant that controls the rate of temperature decrease, 0 < α < 1;

[0027] Furthermore, in step one, the geophysical model parameter vector for the three-layer medium is defined as follows:

[0028] X=[H1,H2,V s1 V s2 V s3 ];

[0029] in,

[0030] H1 is the thickness of the first layer;

[0031] H2 is the thickness of the second layer;

[0032] V s1 The velocity of the first layer of shear waves;

[0033] V s2 This refers to the velocity of the second layer of shear waves;

[0034] V s3 This refers to the velocity of the third layer of transverse waves.

[0035] Furthermore, in step two, the dispersion function is defined as:

[0036] F(V,w)=Re[T surface (2)];

[0037] In the formula,

[0038] F is the dispersion index, which represents the relationship between the phase velocity of the surface wave and the frequency.

[0039] T surface The surface boundary condition vector is obtained through recursive calculation of the layer matrix:

[0040] T surface =M·L 3;

[0041] M=(A2 -1 ·B2)·(A1 -1 ·B1);

[0042] In the formula,

[0043] M is the global matrix;

[0044] A i and B i Let be the layer matrix corresponding to the i-th layer of medium.

[0045] Furthermore, in step three, the scaling factor is defined as:

[0046] .

[0047] Furthermore, in step four, the objective function includes a data fitting term and a shallow penalty term:

[0048]

[0049] The data fitting term uses a weighted multi-mode mean squared error:

[0050]

[0051] Shallow penalty term constrains shallow parameters:

[0052]

[0053] The frequency weighting function is:

[0054] .

[0055] Furthermore, in step five, α takes the value 0.93.

[0056] The beneficial effects of using this invention are:

[0057] Significantly improved shallow layer parameter accuracy: Through a dedicated shallow layer penalty term, the first layer thickness error is reduced from 5-8% to 1-2%, and the first layer velocity error is reduced from 6-10% to 1-3%;

[0058] Significantly improved computational efficiency: Through algorithm optimization, computation time is reduced from 15-30 minutes to 5-10 minutes, improving efficiency by 3 times;

[0059] Enhanced numerical stability: Numerical instability is effectively avoided through matrix condition number checking and complex component detection;

[0060] Full utilization of multi-model data: By allocating reasonable weights, the constraining ability of data from each model at different depths can be fully utilized.

[0061] Results are displayed clearly and completely: The display range is automatically adjusted to ensure that all results are presented in their entirety. Attached Figure Description

[0062] Figure 1 This is a comparative view of the bottom cross-section of the present invention;

[0063] Figure 2 This is a comparison diagram of the dispersion curves of the present invention;

[0064] Figure 3 This is a comparison diagram of the velocity models of the present invention;

[0065] Figure 4 This is a residual distribution diagram of the present invention;

[0066] Figure 5 This is the residual distribution histogram of the present invention;

[0067] Figure 6 This is a data comparison table of the true values ​​and inverted values ​​of the present invention;

[0068] Figure 7 This is a comparison diagram of the depth and shear wave velocity profiles of the present invention;

[0069] Figure 8 This is a schematic diagram illustrating the variation of velocity error with depth according to the present invention;

[0070] Figure 9 This is a schematic diagram illustrating the variation of the relative velocity error with depth according to the present invention;

[0071] Figure 10 This is a parameter sensitivity analysis diagram of the present invention;

[0072] Figure 11 This is the convergence curve of the present invention. Detailed Implementation

[0073] The present invention will now be described in detail with reference to the accompanying drawings.

[0074] Reference Figures 1-9A multi-mode high-precision surface wave dispersion inversion method based on shallow optimization is proposed. By establishing a three-layer geophysical model and adopting a simulated annealing global optimization algorithm, a shallow parameter penalty term and a frequency weighting mechanism are introduced into the objective function, which significantly improves the accuracy of shallow parameter inversion.

[0075] The specific steps include:

[0076] Step 1: Establish a geophysical model;

[0077] A geophysical model incorporating three media layers is established, and the model parameter vector is defined as follows:

[0078] X=[H1,H2,V s1 V s2 V s3 ];

[0079] in,

[0080] H1 is the thickness of the first layer;

[0081] H2 is the thickness of the second layer;

[0082] V s1 The velocity of the first layer of shear waves;

[0083] V s2 This refers to the velocity of the second layer of shear waves;

[0084] V s3 The velocity of the third layer of shear waves;

[0085] Step 2: Forward modeling;

[0086] An improved global matrix method is used to solve the Rayleigh wave dispersion equation. The dispersion function is defined as:

[0087] F(V,w)=Re[T surface (2)];

[0088] In the formula,

[0089] F is the dispersion index, which represents the relationship between the phase velocity of the surface wave and the frequency.

[0090] T surface The surface boundary condition vector is obtained through recursive calculation of the layer matrix:

[0091] T surface =M·L 3;

[0092] M=(A2 -1 ·B2)·(A1 -1 ·B1);

[0093] In the formula,

[0094] M is the global matrix;

[0095] This global matrix is ​​a transformation matrix that connects the boundary conditions at the lowest level of the model with the boundary conditions at the surface. It contains the influence of all intermediate layers on wave propagation and is a description of the overall transmission characteristics of the entire multilayer medium system.

[0096] A i and B i Let i be the layer matrix corresponding to the i-th layer of medium;

[0097] These are determined by the elastic parameters of this layer (such as the shear wave velocity V). s Longitudinal wave velocity V p The stress and displacement fields are determined by the density ρ, thickness h, frequency ω, and wave number k, and are used to describe the transformation relationship of the stress and displacement fields from the top to the bottom of the i-th layer of medium when the wave propagates in the i-th layer.

[0098] Among them, A i Related to the stress and displacement state at the top of the layer, B i Related to the stress and displacement state at the bottom of the layer;

[0099] Multiple layer matrices include:

[0100] Upward wave coefficient matrix E, downward wave attenuation matrix F, upward wave attenuation matrix G, and downward wave coefficient matrix H;

[0101] These matrices together constitute the layer propagation matrix, used to describe the propagation of elastic waves in a layered medium:

[0102] Wave decomposition: Decomposing the wave field into ascending waves and descending waves.

[0103] Boundary matching: ensuring the continuity of displacement and stress at interfaces between different media.

[0104] Propagation effects: describe the phase changes and attenuation of waves as they propagate within a layer.

[0105] In the global matrix method, these layer matrices are multiplied in the order of the layers to form a global matrix M, which is ultimately used to solve the surface wave dispersion equation.

[0106] Solving for the dispersion curve involves finding the combination of phase velocity V and frequency ω that satisfies the surface boundary conditions.

[0107] The i-th layer matrix is ​​constructed as follows:

[0108]

[0109] In the formula,

[0110] μ i Let μ be the shear modulus, the shear modulus of the i-th layer of medium.i =ρ i ·V si ², where ρ i It is density, V si It is the transverse wave velocity;

[0111] k is the horizontal wave number, k = ω / V, where ω is the angular frequency and V is the phase velocity;

[0112] β is the vertical wavenumber ratio. It is related to the vertical wavenumber of the transverse wave;

[0113] The propagation matrix includes P and Q, which describe the phase changes and attenuation of the wave as it propagates within the layer;

[0114] P is the negative sign in the exponential term describing the phase change of a wave as it propagates downwards, indicating that the phase delay of the wave increases with depth.

[0115] Q is the positive sign in the exponential term describing the phase change of a wave as it propagates upwards, indicating the phase change of the wave as the depth decreases.

[0116] ;

[0117] ;

[0118] Step 3: Adaptive phase velocity search;

[0119] Specifically, the search range is dynamically adjusted to V based on the frequency V. min ~V max :

[0120] V min (f)=max(100,min(V s )·γ1(f))

[0121] V max (f)=max(V s )·γ2(f)

[0122] The scaling factor is defined as:

[0123]

[0124] Step 4: Construct a shallow optimization objective function;

[0125] The objective function includes a data fitting term and a shallow penalty term:

[0126]

[0127] The data fitting term uses a weighted multi-mode mean squared error:

[0128]

[0129] Shallow penalty term constrains shallow parameters:

[0130]

[0131] The frequency weighting function is:

[0132]

[0133] Step 5: Simulate annealing for global optimization;

[0134] An improved simulated annealing algorithm is used for inversion, and the temperature T update function is:

[0135] T k+1 =α·T k;

[0136] In the formula,

[0137] k is the number of iterations;

[0138] α is a constant that controls the rate of temperature decrease, 0 < α < 1;

[0139] α close to 1: slow temperature decrease, slow convergence speed, but strong global search capability; temperature decay coefficient.

[0140] When α is close to 0, the temperature drops quickly and the convergence speed is fast, but it may get stuck in a local optimum.

[0141] Preferably, α is 0.93;

[0142] The reheating interval is dynamically adjusted based on the convergence status of the algorithm.

[0143] When the improvement of the objective function value within m consecutive iterations is less than a preset threshold ε, it is determined that the algorithm has fallen into a local extremum, and a reheating process is automatically triggered to ensure the algorithm's global search capability.

[0144] Preferably, the reheating interval is set to 150 iterations to ensure that local extrema are fully escaped.

[0145] Example 1

[0146] Step 1: Establish a geophysical model;

[0147] Parameter boundaries are set as follows:

[0148] Lower bound parameter vector: lb=[12,22,280,450,700];

[0149] Upper bound parameter vector: ub=[18,28,320,550,900];

[0150] Specifically, the lower bound and the upper bound correspond to the minimum and maximum values ​​of the parameter, respectively.

[0151] Step 2: Efficient forward modeling calculation;

[0152] In the forward modeling process, when constructing the layer matrix A and needing to invert it, a numerical stability constant ε = 1 × 10⁻¹² is introduced, and (A) is used. i +ε i ) -1 Instead of direct A i -1 calculate;

[0153] Specifically, a very small value ε is added to each diagonal element of the original matrix A. This operation effectively avoids numerical instability caused by matrix singularity or ill-conditionedness by perturbing the matrix eigenvalues. At the same time, since ε is extremely small, its impact on the accuracy of the physical forward modeling solution is negligible.

[0154] Step 3: Adaptive phase velocity search;

[0155] Dynamically adjust the phase velocity search range and number of scan points based on frequency characteristics:

[0156] High frequency band (f>5Hz): 150 scan points, narrowed search range

[0157] Mid-to-low frequency band: 250 scan points, with a slightly wider search range.

[0158] An improved bisection method was used for root location, with a convergence tolerance of 1×10⁻⁶. -5 The maximum number of iterations is 50.

[0159] Step 4: Construct a shallow optimization objective function;

[0160] The objective function includes a data fitting term and a shallow penalty term, where:

[0161] The mode weights are set to [1.0, 0.7, 0.5], increasing the weights of the base mode.

[0162] Shallow penalty coefficients: α1=50, α2=30, α3=20;

[0163] Prior parameter: H1 prior =15m, V s1 prior =300m / s, V s prior =200m / s;

[0164] Step 5: Simulate annealing for global optimization;

[0165] This is implemented using MATLAB's `simulannealbnd` function. Key parameter settings are as follows:

[0166] Maximum number of iterations: 100;

[0167] Temperature update function: temperaturefast (a temperature scheduling function in the simulated annealing algorithm, used to control the rate at which the temperature decreases during the algorithm's iteration process);

[0168] Reheating interval: 150 iterations;

[0169] Function tolerance: 1×10 -4 ;

[0170] Maximum number of function calculations: 2000.

[0171] Reference Figure 5 Inversion results statistics:

[0172] Objective function value: 99.833039;

[0173] Calculation time: 1.26 minutes;

[0174] Shallow H1 error: 3.19%;

[0175] Shallow Vs1 error: 0.35%;

[0176] Deep Vs3 error: 2.81%.

[0177] The above content is only a preferred embodiment of the present invention. For those skilled in the art, many changes can be made in the specific implementation and application scope based on the concept of the present invention. As long as these changes do not depart from the concept of the present invention, they all fall within the protection scope of the present invention.

Claims

1. A multi-mode high-precision surface wave dispersion inversion method based on shallow optimization, characterized in that: The specific steps include: Step 1: Establish a geophysical model; Establish a geophysical model that includes three layers of media; Step 2: Forward modeling; An improved global matrix method is used to solve the Rayleigh wave dispersion equation; Step 3: Adaptive phase velocity search; The search range is dynamically adjusted to V based on the frequency V. min ~V max : V min (f)=max(100,min(V s )·γ1(f)) V max (f)=max(V s )·γ2(f) γ is the scaling factor; Step 4: Construct a shallow optimization objective function; Step 5: Simulate annealing for global optimization; An improved simulated annealing algorithm is used for inversion, and the temperature T update function is: T k+1 =α·T k; In the formula, k is the number of iterations; α is a constant that controls the rate of temperature decrease, 0 < α < 1.

2. The multi-mode high-precision surface wave dispersion inversion method based on shallow optimization as described in claim 1, characterized in that: In step one, the parameter vector of the geophysical model for the three-layer medium is defined as follows: X=[H1,H2,V s1 ,V s2 ,V s3 ]; in, H1 is the thickness of the first layer; H2 is the thickness of the second layer; V s1 The velocity of the first layer of shear waves; V s2 This refers to the velocity of the second layer of shear waves; V s3 This refers to the velocity of the third layer of transverse waves.

3. The multi-mode high-precision surface wave dispersion inversion method based on shallow optimization as described in claim 1, characterized in that: In step two, the dispersion function is defined as: F(V,w)=Re[T surface (2)]; In the formula, F is the dispersion index, which represents the relationship between the phase velocity of the surface wave and the frequency. T surface The surface boundary condition vector is obtained through recursive calculation of the layer matrix: T surface =M·L 3; M=(A2 -1 ·B2)·(A1 -1 ·B1); In the formula, M is the global matrix; A i and B i Let be the layer matrix corresponding to the i-th layer of medium.

4. The multi-mode high-precision surface wave dispersion inversion method based on shallow optimization as described in claim 1, characterized in that: In step three, the scaling factor is defined as: 。 5. The multi-mode high-precision surface wave dispersion inversion method based on shallow optimization as described in claim 1, characterized in that: In step four, the objective function includes a data fitting term and a shallow penalty term: The data fitting term uses a weighted multi-mode mean squared error: Shallow penalty term constrains shallow parameters: The frequency weighting function is: 。 6. The multi-mode high-precision surface wave dispersion inversion method based on shallow optimization as described in claim 1, characterized in that: In step five, α takes the value 0.93.

Images