A surface wave inversion method for niche particle swarms without pattern recognition

Abstract

The present invention discloses a niche particle swarm surface wave inversion method that does not require pattern recognition. The method obtains measured surface wave dispersion curves from a study area based on field seismic shot records, establishes an initial inversion model based on geological data, uses the measured surface wave dispersion curves to determine a surface wave determinant mismatch function and uses it as the objective function. The objective function is inverted and a niche particle swarm algorithm is used to search for local and global optima. The individual optimal solutions from each inversion are output to form a model solution set. A theoretical surface wave dispersion curve is then determined, and the model solution set is filtered and sorted in ascending order based on a shortest distance fitness function. Cluster analysis is performed on the filtered model solution set based on an improved basic order algorithm. The individual optimal solutions in each cluster are statistically analyzed based on a given allowable error to obtain multiple inversion results. The present invention avoids pattern judgment and root-finding calculations on the surface wave dispersion curve, obtains the multi-solution feature of surface wave inversion, and exhibits excellent robustness and applicability.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic exploration, and in particular relates to a niche particle swarm surface wave inversion method that does not require pattern recognition. Background Art

[0002] Surface wave dispersion curves in layered media are multimodal, with higher-order modes being more sensitive. Furthermore, at the same wavelength, the detection depth of higher-order modes is deeper than that of the fundamental mode. Joint inversion of multi-order mode data can improve the accuracy and precision of shear wave velocity models and ensure the stability of the inversion process. However, in real-world applications, surface wave dispersion energy exhibits multimodal behavior, sometimes making it impossible to correctly distinguish between fundamental and higher-order modes, thus posing challenges to subsequent inversion.

[0003] The particle swarm algorithm (PSO), a collaborative random search method developed by simulating the social behavior of bird flocks, has been widely used to effectively search for optimization problems in various engineering fields. The basic idea of ​​the PSO is that each particle in a swarm shares its individual extreme value with other particles in the swarm, seeking the extreme value of the optimal individual particle as the current global optimal solution for the entire swarm. This allows all particles in the swarm to adjust their speed and position based on their current extreme value and the current global optimal solution shared by the swarm. The local optimal value guides the search direction of particles in the swarm, enhancing the swarm's ability to conduct refined searches.

[0004] Therefore, it is urgent to perform cluster analysis on the inversion results of the surface wave determinant mismatch function by optimizing the basic order of the particle swarm algorithm, and propose a niche particle swarm surface wave inversion method that does not require pattern recognition to find the model space solution closest to the measured dispersion curve from multiple local minima. Summary of the Invention

[0005] In order to solve the problem that the inversion of surface wave dispersion curves in layered media is difficult due to their multimodality, the present invention proposes a niche particle swarm surface wave inversion method that does not require pattern recognition. The dispersion function is used as the inversion objective function, which effectively avoids the mode judgment of the surface wave dispersion curve and the root calculation of the surface wave dispersion curve, improves the inversion speed of multi-mode surface wave inversion, obtains the multi-solution characteristics of surface wave inversion, and has good robustness and applicability.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] A surface wave inversion method for a niche particle swarm without pattern recognition, comprising the following steps:

[0008] Step 1: Select a study area and obtain the measured surface wave dispersion curve based on the field seismic shot gather records in the study area;

[0009] Step 2: Establish the initial inversion model. The initial inversion model is set as a multi-layer elastic horizontal layered medium. According to the geological data of the study area, the density and Poisson's ratio of each elastic horizontal layered medium in the initial inversion model are set. Then, the surface wave determinant mismatch function is determined based on the initial inversion model and the measured surface wave dispersion curve and used as the objective function, as shown in formula (1):

[0010]

[0011] Where obj(m) is the observed value of the inversion parameter of the elastic horizontal layered medium, and m is the model parameter of the elastic horizontal layered medium, including the shear wave velocity V S , longitudinal wave velocity V P , density ρ and formation thickness h, where is the observed value of the phase velocity at the i-th observation point, f i obs is the observed value of the frequency of the i-th observation point, w i is the weight of the i-th observation point, and N is the number of observation points;

[0012] Step 3: Perform multiple inversions on the objective function. In each inversion process, the niche particle swarm algorithm is used to search for the local optimal value and the global optimal value of the objective function, and multiple local optimal values ​​and multiple global optimal values ​​of the objective function are obtained. The individual optimal solution of each particle in the niche particle swarm is output in each inversion calculation to obtain the model solution set of the objective function.

[0013] Step 4: Calculate the shortest distance adaptation function based on the measured surface wave dispersion curve and the theoretical surface wave dispersion curve of the inversion model solution set, and screen and sort the individual optimal solutions in the objective function model solution set based on the shortest distance adaptation function to obtain the screened model solution set;

[0014] Step 5: Based on the improved basic sequential algorithm, cluster analysis is performed on the individual optimal solutions in the filtered model solution set, and the individual optimal solutions in the filtered model solution set are divided into multiple clusters. Combined with the preset allowable error value, the individual optimal solutions within each cluster that are less than the allowable error value are obtained and the average value and variance are calculated, and multiple surface wave inversion results are output.

[0015] Preferably, in step 2, a propagation matrix is ​​established for each elastic horizontal layered medium layer of the inversion initial model, as shown in formula (2):

[0016]

[0017] in,

[0018]

[0019]

[0020]

[0021] Where, T m is the propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the first sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the second sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the third sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the fourth sub-propagation matrix of the mth elastic horizontal layer in the initial inversion model, ρ m is the density of the mth elastic horizontal layered medium layer, V Sm is the shear wave velocity in the mth elastic horizontal layer, V Pm is the longitudinal wave velocity in the mth layer of elastic horizontal layered medium, k is the wave number, and f is the frequency;

[0022] By introducing the interface stiffness matrix S of the mth elastic horizontal layer m and the auxiliary stiffness matrix The recursive relationship is established as:

[0023]

[0024] Where, is the auxiliary stiffness matrix of the mth elastic horizontal layer As shown in formula (7):

[0025]

[0026] in,

[0027]

[0028] Where S m+1 is the interface stiffness matrix S of the m+1th elastic horizontal layer in the initial inversion model m+1 , I is the identity matrix, is the downgoing wave number matrix of the mth elastic horizontal layer in the initial inversion model, as shown in formula (9):

[0029]

[0030] By inverting the interface stiffness matrix S of the m+1th elastic horizontal layer in the initial model m+1Substituting into formula (7) and formula (8), the auxiliary stiffness matrix of the mth elastic horizontal layered medium layer is calculated as Then the auxiliary stiffness matrix of the mth elastic horizontal layer is Substitute into the recursive relation and combine in the half space The interface stiffness matrix of each elastic horizontal layered medium layer in the inversion initial model is obtained by recursive calculation, and the surface wave dispersion function F(f,v,m) of the inversion initial model is obtained, as shown in formula (10):

[0031] F(f,v,m)=|det(S1(f,v,m))| (10)

[0032] Where S1 is the interface stiffness matrix of the first elastic horizontal layered medium in the initial inversion model, v is the phase velocity, and m is the model parameter of the elastic horizontal layered medium.

[0033] Preferably, the step 3 specifically includes the following sub-steps:

[0034] Step 3.1, set the number of inversions M and the population size N of the microhabitat particle swarm;

[0035] Step 3.2, initialize each particle in the niche particle swarm in the search space and randomly set the velocity value of each particle in the niche particle swarm;

[0036] Step 3.3, based on the niche particle swarm algorithm, the objective function is independently inverted M times. During each inversion process, each particle in the niche particle swarm searches in the search space to obtain a local optimal value and shares it in the niche particle swarm. The global optimal value searched by the current niche particle swarm is determined based on the local optimal value searched by each particle. Each particle in the niche particle swarm updates its own position and velocity based on the local optimal value searched by itself and the global optimal value shared by the niche particle swarm, thus obtaining the individual optimal solution of each particle in the niche particle swarm.

[0037] The position update formula of particles in a niche particle swarm is:

[0038]

[0039] Where, is the updated particle position, is the position of the particle before updating, V i d is the velocity of the particle before updating, d is the dimension of the inversion parameters contained in the particle;

[0040] The velocity update formula of particles in a niche particle swarm is:

[0041]

[0042] in,

[0043]

[0044] Where V i ′ is the updated particle velocity, V i d is the velocity of the particle before updating; ω is the inertia factor, which is used to balance the local search weight and global search weight of the niche particle swarm algorithm; nsize is the domain range of the niche particle swarm, which is used to weigh the diversity and convergence speed of the particle swarm in the niche particle swarm algorithm. As the number of inversions increases, nsize is gradually increased from 1 to 4; is a random number between [0, 4.1 / nsize], For all and, nbest j is the jth particle closest to the area where the optimal solution of the i-th particle is located, is the weight of the jth particle;

[0045] In step 3.3, the individual optimal solution of each particle in the inversion niche particle swarm is output to obtain the model solution set of the objective function.

[0046] Preferably, in step 3.1, the particle dimension D is set according to the total number of layers n of the multi-layer elastic horizontal layered medium in the inversion initial model, and the particle dimension D = 2n-1, and then the population number N of the habitat particle group is set according to the particle dimension D, and the population number N of the habitat particle group is N = 10×D.

[0047] Preferably, the step 4 specifically includes the following sub-steps:

[0048] Step 4.1, set the error value of the objective function model solution set;

[0049] Step 4.2: Based on the shortest distance adaptation function, the fitting error between the theoretical surface wave dispersion curve of each individual optimal solution in the objective function model solution set and the measured surface wave dispersion curve is obtained. The fitting error is used to characterize the matching degree E(m) between the theoretical surface wave dispersion curve of the individual optimal solution in the objective function model solution set and the measured surface wave dispersion curve, as shown in formula (14):

[0050]

[0051] Where K is the sequence number of the individual optimal solution in the objective function model solution set, is the observed value of the phase velocity corresponding to the Kth individual optimal solution, is the phase velocity value that is closest to the measured dispersion curve among all modes of the inverted surface wave curve;

[0052] In step 4.3, based on the fitting errors between the theoretical surface wave dispersion curves and the measured surface wave dispersion curves of the individual optimal solutions in the objective function model solution set and the set error value, the individual optimal solutions with fitting errors less than the error value in the objective function optimal solution set are screened out.

[0053] Preferably, the step 5 specifically includes the following sub-steps:

[0054] Step 5.1, set the distance radius and tolerance value;

[0055] Step 5.2: Arrange the individual optimal solutions in the filtered model solution set in ascending order, take the individual optimal solution at the top of the model solution set as the cluster head of the first population, and calculate the Euclidean distance between each individual optimal solution in the model solution set and the cluster head;

[0056] Step 5.3: If the Euclidean distance between the individual optimal solution and the cluster head is not greater than the distance radius, then the individual optimal solutions whose Euclidean distance is not greater than the distance radius are divided into the same population; if the Euclidean distance between the individual optimal solution and the cluster head is greater than the distance radius, then the first individual optimal solution whose Euclidean distance is greater than the distance radius is set as the cluster head of the next population, and the position of the new cluster head is obtained;

[0057] Step 5.4: Based on the position of the new cluster head, calculate the Euclidean distance between the optimal individual solutions of the remaining undivided populations and the new cluster head. Return to step 5.3 until all the optimal individual solutions in the model solution set are divided into their corresponding clusters. This completes the clustering of the optimal individual solutions in the model solution set and proceeds to step 5.5.

[0058] In step 5.5, based on the set allowable error value, the individual optimal solutions within each cluster that are not less than the allowable error value are screened out, and the individual optimal solutions within each cluster that are less than the allowable error value are obtained. The average value and variance of all the individual optimal solutions after screening within each cluster are calculated, and multiple surface wave inversion results are output.

[0059] The beneficial technical effects brought about by the present invention are:

[0060] This paper proposes a niche particle swarm surface wave inversion method that does not require pattern recognition. By improving the Haskell-Thomson transfer matrix algorithm, a niche particle swarm algorithm is used to invert the surface wave determinant mismatch function. Cluster analysis is then performed on the individual optimal solutions in the inverted model solution set, combined with an improved basic sequence algorithm. The model space solution that most closely matches the measured surface wave dispersion curve is searched from multiple local minima. This method uses the dispersion function as the objective function for surface wave inversion, avoiding the need for pattern judgment and root calculation of the surface wave dispersion curve. This method improves the applicability of multi-mode surface wave inversion and significantly accelerates computational speed. It also captures the multi-solution nature of surface wave inversion, demonstrating excellent robustness and applicability. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 This is a flow chart of a surface wave inversion method for a niche particle swarm without pattern recognition according to the present invention.

[0062] Figure 2 Schematic diagram of the surface wave dispersion function of the model containing a low-velocity interlayer. Figure 2 In the figure, (a) is a schematic diagram of the surface wave dispersion function generated based on the Haskell-Thomson method, and (b) is a schematic diagram of the surface wave dispersion function generated based on the improved Haskell-Thomson method.

[0063] Figure 3 Flowchart of the niche particle swarm inversion method.

[0064] Figure 4 This is a comparison chart for testing the effectiveness of the inversion algorithm of the present invention using a model containing a low-velocity interlayer. Figure 4 In the figure, (a) is the simulated shot gather record of the model containing a low-velocity interlayer, and (b) is the dispersion energy diagram of the model containing a low-velocity interlayer generated based on the Hankel transform. The black dots in the figure are the measured dispersion curves automatically picked according to the energy peak, and the solid line is the theoretical dispersion curve.

[0065] Figure 5 These are the inversion results of shear wave velocity and layer thickness when the number of layers is known. Figure 5 In the figure, (a) is the inversion dispersion curve of noise-free data, (b) is the inversion dispersion curve of noise-containing data, (c) is the inversion model of noise-free data, and (d) is the inversion model of noise-containing data.

[0066] Figure 6 The inversion results of shear wave velocity and layer thickness are sorted based on the dispersion function. Figure 6 In the figure, (a) is the inversion dispersion curve containing noisy data, and (b) is the inversion model containing noisy data.

[0067] Figure 7 To refine the inversion results of the layered method. Figure 7 In the figure, (a) is the inversion dispersion curve of noise-free data, (b) is the inversion dispersion curve of noise-containing data, (c) is the inversion model of noise-free data, and (d) is the inversion model of noise-containing data.

[0068] Figure 8 This is a schematic diagram of the layout of the L-shaped survey line.

[0069] Figure 9 This is a comparison chart for testing the effectiveness of the inversion algorithm of the present invention using the measured area 1. Figure 9 In the figure, (a) is a 37-channel vertical component shot gather record, (b) is a dispersion energy map generated based on Hankel transform, (c) is a dispersion energy map generated based on the extended SPAC method, and (d) is a dispersion energy map of the combined active and passive sources.

[0070] Figure 10 It is the P-wave reference model obtained based on the refraction wave time-distance curve and the acoustic wave time difference logging curve.

[0071] Figure 11 This is the inversion result of the L-shaped arrangement of measuring points. Figure 11 In the figure, (a) is V S Wave velocity structure diagram, (b) is a comparison diagram of the theoretical dispersion curve of the inversion model and the measured dispersion curve.

[0072] Figure 12 Schematic diagram of testing the inversion algorithm of the present invention using measured area 2. Figure 12 In the figure, (a) is a 24-channel shot collection record of a highway roadbed measurement point, and (b) is a dispersion energy diagram generated based on the Hankel transform. The energy value in the figure is the result of normalization and multiplication by the power of 3. The black dots are the measured dispersion curves automatically picked according to the energy peak.

[0073] Figure 13 This is the inversion result of 2 roadbed measurement points in the measured area. Figure 13 In the figure, (a) is V S Wave velocity structure diagram, the dotted line is the borehole shear wave logging curve, and (b) is the comparison diagram of the theoretical dispersion curve of the inversion model and the measured dispersion curve. DETAILED DESCRIPTION

[0074] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0075] This paper proposes a surface wave inversion method for a niche particle swarm without pattern recognition. Figure 1 As shown, the specific steps include:

[0076] Step 1: Select the study area and obtain the measured surface wave dispersion curve based on the field seismic shot gather records in the study area.

[0077] Step 2: Establish an initial inversion model. The initial inversion model is set as a multi-layer elastic horizontal layered medium. The density and Poisson's ratio of each elastic horizontal layered medium in the initial inversion model are set. The propagation matrix is ​​established for each elastic horizontal layered medium layer in the initial inversion model, as shown in formula (2):

[0078]

[0079] in,

[0080]

[0081]

[0082]

[0083] Where, T m is the propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the first sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the second sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the third sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the fourth sub-propagation matrix of the mth elastic horizontal layer in the initial inversion model, ρ m is the density of the mth elastic horizontal layered medium layer, V Sm is the shear wave velocity in the mth elastic horizontal layer, V Pm is the longitudinal wave velocity in the mth elastic horizontal layer, k is the wave number, and f is the frequency.

[0084] In order to solve the high-frequency numerical overflow problem of the Haskell-Thomson algorithm, based on the improved Haskell-Thomson algorithm, the exponential growth term can be effectively removed by introducing an auxiliary stiffness matrix, thereby fundamentally avoiding the high-frequency numerical overflow problem. In this embodiment, the interface stiffness matrix S of the mth elastic horizontal layered medium is introduced. m and the auxiliary stiffness matrix The recursive relationship is established as:

[0085]

[0086] Where, is the auxiliary stiffness matrix of the mth elastic horizontal layer As shown in formula (7):

[0087]

[0088] in,

[0089]

[0090] Where S m+1 is the interface stiffness matrix S of the m+1th elastic horizontal layer in the initial inversion model m+1 , I is the identity matrix, is the downgoing wave number matrix of the mth elastic horizontal layer in the initial inversion model, as shown in formula (9):

[0091]

[0092] By inverting the interface stiffness matrix S of the m+1th elastic horizontal layer in the initial model m+1 Substituting into formula (7) and formula (8), the auxiliary stiffness matrix of the mth elastic horizontal layered medium layer is calculated as Then the auxiliary stiffness matrix of the mth elastic horizontal layer is Substitute into the recursive relation and combine in the half space The interface stiffness matrix of each elastic horizontal layered medium layer in the inversion initial model is obtained by recursive calculation, and the surface wave dispersion function F(f,v,m) of the inversion initial model is obtained, as shown in formula (10):

[0093] F(f,v,m)=|det(S1(f,v,m))| (10)

[0094] Where S1 is the interface stiffness matrix of the first elastic horizontal layered medium in the initial inversion model, v is the phase velocity, and m is the model parameter of the elastic horizontal layered medium.

[0095] Since the surface wave determinant mismatch function is equal to the minimum value 0 only when all points (v, f) are surface wave eigenvalues, the surface wave determinant mismatch function can simultaneously perform multi-mode inversion without the need for pattern recognition. That is, the method of the present invention obtains the theoretical surface wave dispersion curve of the inversion initial model by obtaining all zero points of the surface wave dispersion function of the inversion initial model, without the need to solve the dispersion curve, thereby greatly improving the speed of the inversion calculation.

[0096] According to the initial inversion model, the theoretical surface wave dispersion curve and the measured surface wave dispersion curve of the field seismic shot gather records in the study area are obtained to determine the surface wave determinant mismatch function and use it as the objective function, as shown in formula (1):

[0097]

[0098] Where obj(m) is the observed value of the inversion parameter of the elastic horizontal layered medium, and m is the model parameter of the elastic horizontal layered medium, including the shear wave velocity VS , longitudinal wave velocity V P , density ρ and formation thickness h, where is the observed value of the phase velocity at the i-th observation point, f i obs is the observed value of the frequency of the i-th observation point, w i is the weight of the i-th observation point, and N is the number of observation points. In this embodiment, weights are used to measure the uncertainty of data, and the weights are set to 1, which means that the weights of all observation data are equal.

[0099] Taking the surface wave dispersion function of the model with low-velocity interlayer as an example, the setting parameters of the model with low-velocity interlayer are shown in Table 1.

[0100] Table 1 Formation parameters of the model with low-velocity interlayer

[0101]

[0102] The zero point of the surface wave dispersion function corresponds to its dispersion curve. Figure 2 The dispersion function diagrams generated by the Haskell-Thomson method before and after the improvement are shown. The dispersion function diagram generated by the Haskell-Thomson method begins to experience numerical overflow at approximately 45 Hz, while the dispersion function diagram generated by the improved Haskell-Thomson method effectively avoids this numerical overflow problem. Furthermore, the dispersion function diagram generated by the Haskell-Thomson method includes not only the dispersion curves of surface wave modes with displacement on the surface, but also the dispersion curves of channel waves with displacement only in the low-velocity layer, i.e., a full-mode solution, while the improved Haskell-Thomson method only includes the dispersion curves of surface wave modes with displacement on the surface. Therefore, for methods such as surface wave exploration and background noise imaging that deploy detectors on the surface, the improved Haskell-Thomson method of the present invention is more suitable as a forward algorithm.

[0103] Sensitivity analysis of the surface wave dispersion function shows that several local minima exist depending on the shallow shear wave velocity. To avoid falling into local minima, the inversion generally starts with a model where the first layer velocity is greater than the true velocity. However, when no high-frequency information is available, such as in the presence of a low-velocity interlayer, the true velocity of the first layer in the inversion model cannot be accurately obtained. The inversion of the surface wave determinant mismatch function is essentially a multi-peak function optimization problem. To reduce the solution's dependence on the initial inversion model, a niche particle swarm algorithm is used to search for multiple local extremes of the objective function and has strong fine-grained search capabilities. Therefore, without restricting the velocity value of the first layer of the model, the model space solution closest to the measured dispersion curve can be found from multiple local minima.

[0104] Step 3: Use the niche particle swarm inversion method to perform multiple inversions on the objective function, such as Figure 3 As shown in the figure, in each inversion process, the niche particle swarm algorithm is used to search for the local optimal value and the global optimal value of the objective function, and multiple local optimal values ​​and multiple global optimal values ​​of the objective function are obtained. The individual optimal solution of each particle in the niche particle swarm is output for each inversion calculation, and the model solution set of the objective function is obtained. Specifically, the following steps are included:

[0105] In step 3.1, set the number of inversions M = 50 and the population number N of the microhabitat particle swarm. Set the particle dimension D according to the total number of layers n of the multi-layer elastic horizontal layered medium in the inversion initial model, and the particle dimension D = 2n-1. Then, set the population number N of the microhabitat particle swarm according to the particle dimension D, and the population number N of the microhabitat particle swarm = 10×D.

[0106] In step 3.2, each particle in the microhabitat particle swarm is initialized in the search space, and the velocity value of each particle in the microhabitat particle swarm is randomly set.

[0107] In step 3.3, the objective function is inverted 50 times independently based on the niche particle swarm algorithm. During each inversion process, each particle in the niche particle swarm searches in the search space to obtain the local optimal value and shares it in the niche particle swarm. The global optimal value searched by the current niche particle swarm is determined based on the local optimal value searched by each particle. Each particle in the niche particle swarm updates its own position and velocity based on the local optimal value searched by itself and the global optimal value shared by the niche particle swarm, and obtains the individual optimal solution of each particle in the niche particle swarm.

[0108] The position update formula of particles in a niche particle swarm is:

[0109]

[0110] Where, is the updated particle position, is the position of the particle before updating, V i d is the velocity of the particle before updating, and d is the dimension of the inversion parameters contained in the particle.

[0111] The velocity update formula of particles in a niche particle swarm is:

[0112]

[0113] in,

[0114]

[0115] Where V i ′ is the updated particle velocity, V id is the velocity of the particle before the update; ω is the inertia factor, which is used to balance the local search weight and the global search weight of the niche particle swarm algorithm. In this embodiment, ω = 0.729843788; nsize is the domain range of the niche particle swarm, which is used to balance the diversity and convergence speed of the particle swarm in the niche particle swarm algorithm. As the number of inversions increases, nsize is gradually increased from 1 to 4; is a random number between [0, 4.1 / nsize], For all and, nbest j is the jth particle closest to the area where the optimal solution of the i-th particle is located, is the weight of the jth particle.

[0116] In step 3.3, the individual optimal solution of each particle in the inversion niche particle swarm is output to obtain the model solution set of the objective function. In this embodiment, the model solution set of the objective function contains 50×10×D individual optimal solutions.

[0117] Step 4: Calculate the shortest distance adaptation function based on the measured surface wave dispersion curve and the theoretical surface wave dispersion curve of the inversion model solution set, and screen and sort the individual optimal solutions in the objective function model solution set based on the shortest distance adaptation function to obtain the screened model solution set. This specifically includes the following sub-steps:

[0118] Step 4.1, set the error value of the objective function model solution set.

[0119] Step 4.2: Based on the shortest distance adaptation function, the fitting error between the inverted surface wave dispersion curve of each individual optimal solution in the objective function model solution set and the measured surface wave dispersion curve is obtained. The fitting error is used to characterize the matching degree E(m) between the inverted surface wave curve of the individual optimal solution in the objective function model solution set and the measured surface wave dispersion curve, as shown in formula (14):

[0120]

[0121] Where K is the sequence number of the individual optimal solution in the objective function model solution set, is the observed value of the phase velocity corresponding to the Kth individual optimal solution, is the phase velocity value that is closest to the measured dispersion curve among all modes of the inverted surface wave curve.

[0122] In step 4.3, based on the fitting errors between the theoretical surface wave dispersion curves and the measured surface wave dispersion curves of the individual optimal solutions in the objective function model solution set and the set error value, the individual optimal solutions with fitting errors less than the error value in the objective function optimal solution set are screened out.

[0123] Step 5, taking into account the presence of error data in the actual observation data, the best fit data may not necessarily be able to obtain the best inversion result. Therefore, for the measured dispersion curve containing errors, after inversion optimization, the individual optimal solution in the model solution set is not necessarily the true model solution, and the true model solution often corresponds to a local extreme value. Therefore, relative to directly selecting the best fit solution as the final inversion result, the present embodiment carries out cluster analysis on all individual optimal solutions and then launches further inversion interpretation. Due to the lack of prior information on population types in the cluster analysis, the present embodiment adopts a kind of improved basic sequence algorithm to carry out cluster analysis on the individual optimal solutions in the objective function model solution set.

[0124] Based on the improved basic sequential algorithm, cluster analysis is performed on the individual optimal solutions in the filtered model solution set. The individual optimal solutions in the filtered model solution set are divided into multiple clusters. Combined with the preset allowable error value, the individual optimal solutions within each cluster that are less than the allowable error value are obtained and the average and variance are calculated. Multiple surface wave inversion results are output. The specific sub-steps include the following:

[0125] Step 5.1, set the distance radius (in this embodiment, the distance radius is set to ) and the allowable error value.

[0126] In step 5.2, the individual optimal solutions in the filtered model solution set are arranged in ascending order, the individual optimal solution at the first position in the model solution set is used as the cluster head of the first population, and the Euclidean distance between each individual optimal solution in the model solution set and the cluster head is calculated respectively.

[0127] In step 5.3, if the Euclidean distance between the individual optimal solution and the cluster head is not greater than the distance radius, the individual optimal solutions whose Euclidean distance is not greater than the distance radius are divided into the same population; if the Euclidean distance between the individual optimal solution and the cluster head is greater than the distance radius, the first individual optimal solution whose Euclidean distance is greater than the distance radius is set as the cluster head of the next population, and the position of the new cluster head is obtained.

[0128] In step 5.4, based on the position of the new cluster head, calculate the Euclidean distance between the individual optimal solutions of the remaining undivided populations and the new cluster head, and return to step 5.3 until all individual optimal solutions in the model solution set are divided into their corresponding clusters, completing the cluster division of the individual optimal solutions in the model solution set and proceeding to step 5.5.

[0129] In step 5.5, based on the set allowable error value, the individual optimal solutions within each cluster that are not less than the allowable error value are screened out, and the individual optimal solutions within each cluster that are less than the allowable error value are obtained. The average value and variance of all the individual optimal solutions after screening within each cluster are calculated, and multiple surface wave inversion results are output.

[0130] Example 1

[0131] This example uses the low-velocity interlayer model shown in Table 1 as an example to test the effectiveness of the inversion algorithm. In the low-velocity interlayer model, the first layer is a surface high-velocity layer, and the second layer contains saturated soil with a high Poisson's ratio. This typical model will cause the splitting of the Rayleigh wave mode. The discrete wave number method based on the stiffness matrix method is used to simulate 200 shot records. Figure 4 As shown in (a), the simulation parameters are as follows: the time step is set to 0.001s, the offset is set to 1m, the trace spacing is set to 1m, the maximum offset is set to 200m, the source wavelet is a Ricker wavelet with a main frequency of 20Hz, a delay of 0.1s, and a maximum frequency of 500Hz. The actual calculation frequency band ranges from 0 to 100Hz (conventional seismic exploration frequency band), and the part greater than 100Hz is filled with zero values. Figure 4 (b) shows the dispersion energy diagram of a model containing a low-velocity interlayer generated based on the Hankel transform. The 200m array shown in the figure reliably recovers surface wave information of approximately the same wavelength. The picking range is set to 2-30Hz, and the frequency spacing is set to 0.5Hz. Comparing the inverted dispersion curve with the picked dispersion curve reveals that the picked dispersion curve contains six modes. Starting at 15Hz, the picked dispersion curve smoothly transitions from the third to the fifth highest-order mode, making it difficult to correctly identify each mode from the dispersion energy diagram.

[0132] Therefore, two layered methods are used to perform inversion tests:

[0133] (1) Inversion of S-wave velocity and layer thickness with known number of layers

[0134] Assuming that the number of layers, P-wave velocity and density of the initial inversion model are known, the lower limit of the S-wave velocity of each layer is set to 100 m / s, and the upper limit is set to the highest, that is, The layer thickness range is set to 1 to 40 m, which is twice the actual value. Under the condition of a very wide model range, the dispersion curves without noise and with 5% random noise interference are used to test the proposed method for surface wave inversion of niche particle swarms without pattern recognition.

[0135] For ease of analysis, only the first cluster was considered, and the tolerance was set to the individual optimal value of the cluster head of the second cluster. Individual optimal values ​​in the first cluster exceeding the tolerance were removed. For noise-free data, due to the interference of non-plane wave propagation and body wave energy, there was a slight difference between the theoretical dispersion curve and the shortest distance mismatch function. The shortest distance mismatch function between the two was 1.63 m / s, and the number of inversions was set to 200. The tolerance was set to 5.13 m / s, which is equal to the optimal value of the cluster head of the second cluster. The inversion results showed that the first 140 individual optimal values ​​of the first cluster were below this tolerance. For noisy data, the shortest distance mismatch function between the theoretical dispersion curve and the first 36 individual optimal values ​​of the first cluster were below this tolerance. The number of inversions was set to 400 and the tolerance was set to 5.94 m / s.

[0136] The statistical analysis results of the inversion results are shown in Table 2. Figure 5 The figure shows the inversion results of shear wave velocity and layer thickness when the number of layers is known. The average value of the individual optimal solutions that meet the conditions is taken as the final inversion result.

[0137] Table 2 Statistics of dispersion curve inversion results of three-layer model without noise and with noise interference

[0138]

[0139] The analysis shows that, no matter under noise-free conditions or under noise-containing conditions, the final inverted dispersion curve (such as Figure 5 (a) and Figure 5 (b) and the measured dispersion curve (shown as the solid line in Figure 5 (a) and Figure 5 (b) shows a good fit), and the inversion models obtained in both cases are in good agreement with the true model (e.g. Figure 5 (c) and Figure 5 (d)). In the noise-free case, the average relative error between the inversion model and the true model is 4.11%, while in the noise-containing case, the average relative error between the inversion model and the true model is 0.50%. Of course, the fitting error of the noise-free inversion result is smaller than that of the noisy inversion result because the noisy inversion result only contains 36 individuals, so its average error and mean square error are smaller. Then, a cluster analysis was performed on all individuals based on the determinant misfit function value, as shown in Figure 2. Figure 6 As shown, Figure 6 The results of the first cluster analysis are shown. It can be seen that the model is consistent with the theoretical model (such as Figure 6 (b)) is inconsistent, resulting in a large difference between its corresponding dispersion curve and the measured value (such as Figure 6(a)). After analyzing the clustering results in sequence, it is found that the previous inversion result containing noise is only ranked in the 9th cluster, which means that its corresponding determinant mismatch function value is not a global minimum, but a local minimum. For the noise-free inversion result, this is not the case. Its determinant mismatch function value and the shortest distance mismatch function value are both global minima.

[0140] (2) S-wave velocity inversion using the refined layering method

[0141] In surface wave exploration, refined layering is often used to approximate subsurface structures. Commonly used refined layering schemes include equal-thickness layering and variable-thickness layering. Assuming the P-wave velocity and density of the inversion model are known, each layer is subdivided into four layers with a thickness of 5 m. The S-wave velocity range of each layer is consistent with the S-wave velocity range listed in Table 1. The proposed method for surface wave inversion based on a niche particle swarm without pattern recognition was tested using dispersion curves with and without noise and with 5% random noise.

[0142] For noise-free data, the number of inversions is set to 200, the allowable error value is set to 2.93m / s, and the optimal values ​​of the first 73 individuals in the first cluster in the inversion results are lower than the allowable error value. For noisy data, the number of inversions is set to 400, the allowable error value is set to 5.92m / s, and the optimal values ​​of the first 42 individuals in the first cluster in the inversion results are lower than the allowable error value. The statistical analysis results of the inversion results are shown in Figure 7 , the average value of the individual optimal solutions that meet the conditions is taken as the final inversion result. The analysis shows that, whether under noise-free conditions or under noise-containing conditions, the final inverted dispersion curve (such as Figure 7 (a) and Figure 7 (b) and the measured dispersion curve (shown as the solid line in Figure 7 (a) and Figure 7 The inversion models obtained in the two cases are basically consistent with the true model (as shown by the solid points in (b)). Figure 7 (c) and Figure 7 (d)). In the noise-free case, the average relative error between the inversion model and the true model is 5.67%, while in the noise-containing case, the average relative error between the inversion model and the true model is 6.47%, thus proving the effectiveness of the refined layering strategy in the proposed method for surface wave inversion of niche particle swarms without pattern recognition.

[0143] Application Experiment

[0144] The surface wave inversion method of a niche particle swarm without pattern recognition proposed in the present invention is applied to two measured areas respectively. The measured area 1 is detected by combining active source and micro-motion detection, and the measured area 2 is the detection of a highway roadbed.

[0145] Measurement area 1

[0146] An L-shaped survey line was laid out in Shijiazhuang Village in the north of Zhuozhou City, Hebei Province. Figure 8 As shown, the detector frequency was set to 4 Hz, the trace spacing was set to 2 m, the sampling interval was set to 1 ms, and a drop-hammer seismic source was used, with the shot point located at the intersection of two vertically aligned lines. The micro-seismic sampling interval was set to 1 ms, and the acquisition time was 2 hours. The processing parameters of the extended SPAC method were set as follows: downsampling to 200 Hz, data segment length of 60 s, overlap rate of 50%, and bandpass filtering of 0.5 to 30 Hz.

[0147] Depend on Figure 9 (a) and Figure 9 As can be seen from (b), the energy of the high-order mode of Rayleigh wave starts to dominate from 33 Hz. Figure 9 (b) and Figure 9 (c) shows that the two are consistent in the common frequency band of 8 to 30 Hz, which further proves the reliability of the results obtained by the extended SPAC method. The active source and passive source dispersion energy maps are spliced ​​together, as shown in Figure 9 (d) shows that 3-8 Hz is the result of the extended SPAC method, and 8-60 Hz is the result of the active source. The measured dispersion curve automatically picked up according to the energy peak in the figure is shown in Figure 9 As shown in the black dots in (d).

[0148] The dispersion curve is much more sensitive to S-wave velocity and layer thickness than to P-wave velocity and density. However, if the P-wave velocity model deviates too much from the true model, it will cause large errors or non-convergence of the inversion results. Therefore, we establish an approximate P-wave velocity model based on the borehole data near the measurement point and the refraction wave information on the shot gather record, such as Figure 10 As shown, it is assumed that the P-wave velocity and layer thickness remain unchanged during the inversion process.

[0149] The maximum detection depth is set to half of the maximum wavelength, that is, the maximum detection depth is set to 70m. In order to avoid the appearance of velocity models that are inconsistent with reality, the lower limit of the S-wave velocity of each layer is set to 100m / s, and the upper limit is set to 500m / s and 500m / s, except for the half-space. The smaller value between the two, the range of half space is set to 500~1000m / s, and the inversion initial model parameters are shown in Table 3.

[0150] Table 3 Initial model of L-shaped arrangement measuring points

[0151]

[0152] Since the surface wave inversion is non-unique, after 200 inversion calculations and cluster analysis of the inversion results, several different optimal model solutions were obtained, and the corresponding dispersion curves all fit well with the measured dispersion curves. Figure 11 (a) The dotted line shows the inversion result. The tolerance value is set to 4 m / s and the first 47 optimal solutions of the third cluster are selected as the final inversion results ( Figure 11 (a) The solid line shows that the inversion model profile is in good agreement with the drilling data. Figure 11 (b) shows a comparison of the theoretical dispersion curves from the inversion model and the measured dispersion curves, showing a good fit. The inversion successfully identified a low-velocity clayey silt layer between 10 and 20 m and a low-velocity sandstone interlayer between 60 and 70 m. Below 70 m, near the bottom of the borehole, the reliability of the shear wave logging curve deteriorates due to shallow mud deposition. Furthermore, below 70 m, the maximum depth of surface wave detection is exceeded, resulting in significant deviations from the suspended shear wave logging curve.

[0153] Measurement area 2

[0154] Taking a highway subgrade as an example, the proposed microhabitat particle swarm surface wave inversion method without pattern recognition is used. The goal is to use Rayleigh waves to detect whether there is a low-speed soft interlayer in the subgrade. The acquisition parameters include: the instrument uses a domestic SE2404 seismograph, the detector is a 4.5Hz vertical component low-frequency detector, the excitation method is hammer excitation, the number of channels is 24, the recording time is set to 0.5s, the sampling interval is set to 0.5ms, the offset is set to 7.0m, and the channel spacing is set to 1.0m. Figure 12 It can be seen that since the measured data are almost in the same curve position and there is no obvious mode discrimination mark, it is easy to uniformly regard all data points as the fundamental wave mode.

[0155] Throughout the surface wave inversion process, the Poisson's ratio is assumed to be constant. The P-wave velocity can be calculated from the Poisson's ratio and the S-wave velocity. The layer ratio is set to 1.5, and the parameters for establishing the initial inversion model are shown in Table 4.

[0156] Table 4 Initial model of highway roadbed measurement points

[0157]

[0158] After 200 inversion calculations and cluster analysis of the inversion results, several different optimal model solutions were obtained, and the corresponding dispersion curves all fitted well with the measured dispersion curves. Figure 13 (a) (shown by the dotted line), we set the tolerance to 3 m / s and selected the first 170 optimal solutions of the fourth cluster as the final inversion results ( Figure 13(a) The solid line shows that the inversion model profile is in good agreement with the drilling data. Figure 13 (b) shows a comparison of the theoretical dispersion curve of the inversion model and the measured dispersion curve, which shows a good fit. The inversion successfully identified a low-velocity interlayer between 3 and 6 m.

[0159] The test results based on noise-free and noise-containing theoretical data verify that the method of the present invention has good robustness. At the same time, based on two typical case surveys of measured area 1 and measured area 2, the effectiveness and applicability of the method of the present invention are demonstrated.

[0160] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the scope of protection of the present invention.

Claims

1. A surface wave inversion method for a niche particle swarm without pattern recognition, characterized in that: The specific steps include: Step 1: Select a study area and obtain the measured surface wave dispersion curve based on the field seismic shot gather records in the study area; Step 2: Establish the initial inversion model. The initial inversion model is set as a multi-layer elastic horizontal layered medium. According to the geological data of the study area, the density and Poisson's ratio of each elastic horizontal layered medium in the initial inversion model are set. Then, the surface wave determinant mismatch function is determined based on the initial inversion model and the measured surface wave dispersion curve and used as the objective function, as shown in formula (1): Where obj(m) is the observed value of the inversion parameter of the elastic horizontal layered medium, and m is the model parameter of the elastic horizontal layered medium, including the shear wave velocity V S , longitudinal wave velocity V P , density ρ and formation thickness h, where is the observed value of the phase velocity at the i-th observation point, f i obs is the observed value of the frequency of the i-th observation point, w i is the weight of the i-th observation point, and N is the number of observation points; Step 3: Perform multiple inversions on the objective function. In each inversion process, the niche particle swarm algorithm is used to search for the local optimal value and the global optimal value of the objective function, and multiple local optimal values ​​and multiple global optimal values ​​of the objective function are obtained. The individual optimal solution of each particle in the niche particle swarm is output in each inversion calculation to obtain the model solution set of the objective function. Step 4: Calculate the shortest distance adaptation function based on the measured surface wave dispersion curve and the theoretical surface wave dispersion curve of the inversion model solution set, and screen and sort the individual optimal solutions in the objective function model solution set based on the shortest distance adaptation function to obtain the screened model solution set; Step 5: Based on the improved basic sequential algorithm, cluster analysis is performed on the individual optimal solutions in the filtered model solution set, and the individual optimal solutions in the filtered model solution set are divided into multiple clusters. Combined with the preset allowable error value, the individual optimal solutions within each cluster that are less than the allowable error value are obtained and the average value and variance are calculated, and multiple surface wave inversion results are output.

2. The method for surface wave inversion in a niche particle swarm without pattern recognition according to claim 1, characterized in that: In step 2, a propagation matrix is ​​established for each elastic horizontal layered medium layer of the inversion initial model, as shown in formula (2): in, Where, T m is the propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the first sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the second sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the third sub-propagation matrix of the mth elastic horizontal layer in the inversion initial model, is the fourth sub-propagation matrix of the mth elastic horizontal layer in the initial inversion model, ρ m is the density of the mth elastic horizontal layered medium layer, V Sm is the shear wave velocity in the mth elastic horizontal layer, V Pm is the longitudinal wave velocity in the mth layer of elastic horizontal layered medium, k is the wave number, and f is the frequency; By introducing the interface stiffness matrix S of the mth elastic horizontal layer m and the auxiliary stiffness matrix The recursive relationship is established as: Where, is the auxiliary stiffness matrix of the mth elastic horizontal layer As shown in formula (7): in, Where S m+1 is the interface stiffness matrix S of the m+1th elastic horizontal layer in the initial inversion model m+1 , I is the identity matrix, is the downgoing wave number matrix of the mth elastic horizontal layer in the initial inversion model, as shown in formula (9): By inverting the interface stiffness matrix S of the m+1th elastic horizontal layer in the initial model m+1 Substituting into formula (7) and formula (8), the auxiliary stiffness matrix of the mth elastic horizontal layered medium layer is calculated as Then the auxiliary stiffness matrix of the mth elastic horizontal layer is Substitute into the recursive relation and combine in the half space The interface stiffness matrix of each elastic horizontal layered medium layer in the inversion initial model is obtained by recursive calculation, and the surface wave dispersion function F(f,v,m) of the inversion initial model is obtained, as shown in formula (10): F(f,v,m)=|det(S1(f,v,m))| (10) Where S1 is the interface stiffness matrix of the first elastic horizontal layered medium in the initial inversion model, v is the phase velocity, and m is the model parameter of the elastic horizontal layered medium.

3. The method for surface wave inversion in a niche particle swarm without pattern recognition according to claim 1, characterized in that: The step 3 specifically includes the following sub-steps: Step 3.1, set the number of inversions M and the population size N of the microhabitat particle swarm; Step 3.2, initialize each particle in the niche particle swarm in the search space and randomly set the velocity value of each particle in the niche particle swarm; Step 3.3, based on the niche particle swarm algorithm, the objective function is independently inverted M times. During each inversion process, each particle in the niche particle swarm searches in the search space to obtain a local optimal value and shares it in the niche particle swarm. The global optimal value searched by the current niche particle swarm is determined based on the local optimal value searched by each particle. Each particle in the niche particle swarm updates its own position and velocity based on the local optimal value searched by itself and the global optimal value shared by the niche particle swarm, thus obtaining the individual optimal solution of each particle in the niche particle swarm. The position update formula of particles in a niche particle swarm is: Where, X i ' d is the updated particle position, is the position of the particle before updating, V i d is the velocity of the particle before updating, d is the dimension of the inversion parameters contained in the particle; The velocity update formula of particles in a niche particle swarm is: in, Where V i ′ is the updated particle velocity, V i d is the velocity of the particle before updating; ω is the inertia factor, which is used to balance the local search weight and global search weight of the niche particle swarm algorithm; nsize is the domain range of the niche particle swarm, which is used to weigh the diversity and convergence speed of the particle swarm in the niche particle swarm algorithm. As the number of inversions increases, nsize is gradually increased from 1 to 4; is a random number between [0, 4.1 / nsize], For all and, nbest j is the jth particle closest to the area where the optimal solution of the i-th particle is located, is the weight of the jth particle; In step 3.3, the individual optimal solution of each particle in the inversion niche particle swarm is output to obtain the model solution set of the objective function.

4. The method for surface wave inversion in a niche particle swarm without pattern recognition according to claim 3, characterized in that: In step 3.1, the particle dimension D is set according to the total number of layers n of the multi-layer elastic horizontal layered medium in the inversion initial model, and the particle dimension D=2n-1. Then, the population number N of the habitat particle group is set according to the particle dimension D, and the population number N of the habitat particle group is set according to the particle dimension D, and the population number N of the habitat particle group is 10×D.

5. The method for surface wave inversion in a niche particle swarm without pattern recognition according to claim 1, characterized in that: The step 4 specifically includes the following sub-steps: Step 4.1, set the error value of the objective function model solution set; Step 4.2: Based on the shortest distance adaptation function, the fitting error between the theoretical surface wave dispersion curve of each individual optimal solution in the objective function model solution set and the measured surface wave dispersion curve is obtained. The fitting error is used to characterize the matching degree E(m) between the theoretical surface wave dispersion curve of the individual optimal solution in the objective function model solution set and the measured surface wave dispersion curve, as shown in formula (14): Where K is the sequence number of the individual optimal solution in the objective function model solution set, is the observed value of the phase velocity corresponding to the Kth individual optimal solution, is the phase velocity value that is closest to the measured dispersion curve among all modes of the inverted surface wave curve; In step 4.3, based on the fitting errors between the theoretical surface wave dispersion curves and the measured surface wave dispersion curves of the individual optimal solutions in the objective function model solution set and the set error value, the individual optimal solutions with fitting errors less than the error value in the objective function optimal solution set are screened out.

6. The method for surface wave inversion in a niche particle swarm without pattern recognition according to claim 1, characterized in that: The step 5 specifically includes the following sub-steps: Step 5.1, set the distance radius and tolerance value; Step 5.2: Arrange the individual optimal solutions in the filtered model solution set in ascending order, take the individual optimal solution at the top of the model solution set as the cluster head of the first population, and calculate the Euclidean distance between each individual optimal solution in the model solution set and the cluster head; Step 5.3: If the Euclidean distance between the individual optimal solution and the cluster head is not greater than the distance radius, then the individual optimal solutions whose Euclidean distance is not greater than the distance radius are divided into the same population; if the Euclidean distance between the individual optimal solution and the cluster head is greater than the distance radius, then the first individual optimal solution whose Euclidean distance is greater than the distance radius is set as the cluster head of the next population, and the position of the new cluster head is obtained; Step 5.4: Based on the position of the new cluster head, calculate the Euclidean distance between the optimal individual solutions of the remaining undivided populations and the new cluster head. Return to step 5.3 until all the optimal individual solutions in the model solution set are divided into their corresponding clusters. This completes the clustering of the optimal individual solutions in the model solution set and proceeds to step 5.

5. In step 5.5, based on the set allowable error value, the individual optimal solutions within each cluster that are not less than the allowable error value are screened out, and the individual optimal solutions within each cluster that are less than the allowable error value are obtained. The average value and variance of all the individual optimal solutions after screening within each cluster are calculated, and multiple surface wave inversion results are output.

Images