Rayleigh surface wave polarization curve calculation method based on joint algorithm

Abstract

The invention discloses a Rayleigh surface wave polarization curve calculation method based on a joint algorithm. The method comprises the following steps: preprocessing a three-component original waveform record; a forward greedy selection mechanism module utilizing population search-migration is used for searching a parameter optimal solution; parameters in the variational mode decomposition algorithm are optimized, so that the variational mode decomposition algorithm adaptively determines a signal decomposition scale, and sub-waveform data with different orthogonalities are obtained; calculating the phase relation of each component of the signal by using an adaptive algorithm with fewer input parameters, and finding out a signal segment meeting the polarization characteristics of the Rayleigh surface wave; calculating polarizability of different frequency bands; and calculating a Rayleigh surface wave polarization curve recorded by the three-component ground pulsation waveforms with different scales. According to the method, the stable Rayleigh surface wave components can be extracted from the complex tremor pulsation waveform in a multi-scale and multi-level mode, calculation of the polarization curve is completed, and technical support is provided for subsequent multi-level and multi-view interpretation of the physical characteristics of the geologic body.

Description

Technical Field

[0001] The present invention relates to the field of Rayleigh wave calculation, and in particular to a method for calculating a Rayleigh wave polarization curve based on a joint algorithm. Background Art

[0002] With the rapid development of urbanization, the demand for underground space resources is increasing, and high-resolution and high-precision exploration of near-surface geological structures has become the safety basis for engineering construction and operation. As a relatively safe and convenient geological non-destructive exploration method, ground pulse monitoring has been increasingly widely used in urban underground space exploration, geological disaster body exploration, and earthquake site classification. Ground pulses, also known as constant micro-motion (in Japan) or background noise (in Europe and the United States), are wave fields with rich environmental seismic noise as the source. Its main seismic phases are surface waves (Rayleigh surface waves and Love surface waves) propagating along the near surface. Among them, Rayleigh surface waves are more widely used in multi-scale geological-geophysical exploration because of their dispersion and polarization characteristics. The polarization characteristics of Rayleigh surface waves refer to the ratio of its horizontal vibration component (H) to its vertical vibration component (V) in different frequency bands, and are therefore expressed in the form of a frequency (X-axis)-polarizability-H / V (Y-axis) curve, while the morphological changes of the polarization curve are highly correlated with the elastic wave velocity, density, layer thickness and Poisson's ratio of the formation. Therefore, how to extract Rayleigh surface waves in the in-situ monitoring of ground pulsation at multiple levels and calculate its polarization characteristic curve (H / V curve) is extremely important for inverting the optimal geological body physical parameter model (elastic wave velocity, density, layer thickness, etc.) and ultimately interpreting the characteristics of the geological body.

[0003] At present, the core of the calculation of the polarization characteristic curve of Rayleigh surface wave is to convert the original three-component (east-west / north-south / vertical) observed vibration waveform (time domain) into a three-component spectral amplitude curve (frequency domain) through numerical integral transformation, and then synthesize the horizontal two-component amplitude curve into a horizontal amplitude curve, and finally compare it with the vertical amplitude curve to obtain the polarization characteristic curve. The X-axis of the curve is frequency and the Y-axis is H / V value (polarizability). The peak of the polarization curve represents the more obvious change in the physical parameters of the geological body, such as the inherent properties such as wave velocity and stiffness. However, the shortcomings of the current algorithm are mainly in the following two aspects: 1. In terms of gain Rayleigh surface wave components, the continuous wavelet transform algorithm is currently used. Although the algorithm can take into account the signal energy resolution of the time domain and the frequency domain at the same time, the wavelet basis functions adopted by the algorithm are of various types, and the parameters of the wavelet basis functions need to be set manually. Therefore, the parameters of the algorithm are semi-empirically set, which has a certain impact on the calculation results of the final polarization curve. 2. The scale of natural geological bodies responds differently to the vibration wave fields of noise sources in different frequency bands. Therefore, the surface wave polarization curves presented in different areas of the geological body are different in form. The current extraction and calculation of polarization curves are based on aliased waveform signals. Although the algorithm can separate the interference signals of specific frequency bands through filtering modules, the remaining signal energy still cannot exist orthogonally. Therefore, the polarization curves also present a mixed state of multiple forms and cannot be separated into multi-scale polarization curves, reflecting the physical property characteristics of multi-scale geological bodies. Summary of the invention

[0004] In order to overcome the shortcomings and deficiencies of the prior art, the present invention provides a method for calculating the Rayleigh surface wave polarization curve based on a joint algorithm. This invention mainly improves the current extraction technology of the Rayleigh surface wave polarization curve for the ground pulsation source. The purpose of this invention is as follows: 1. By optimizing the population algorithm, the in-situ observed waveform data is separated into multi-scale independent orthogonal time series components to minimize signal correlation; 2. By using the improved zero-crossing phase analysis method, waveform time series of different scales are calculated to obtain Rayleigh surface wave polarization curves of different scales;

[0005] Compared with the overall polarization curve obtained by the traditional algorithm, the present invention can obtain multi-scale polarization curves, and then carry out multi-resolution geological structure physical characteristics inversion, reduce the inversion error rate, and then interpret the geological body structure characteristics more comprehensively and accurately.

[0006] A method for calculating a Rayleigh wave polarization curve based on a joint algorithm, the method comprising:

[0007] Step S1, pre-processing the original waveform records of the three components of the Rayleigh surface wave polarization characteristic curve by removing instrument response, removing mean, removing potential terms, and filtering;

[0008] Step S2, using a population search-migration forward greedy selection mechanism module to find the optimal solution for parameters, providing support for subsequent decomposition of waveform data;

[0009] Step S3, combining the population optimization algorithm in step S2 with the fitness function of minimum envelope entropy, optimizing the parameters K and α in the variational mode decomposition algorithm, making the variational mode decomposition algorithm adaptively determine the signal decomposition scale, decomposing the original data waveform in step S1, and obtaining sub-waveform data of different orthogonality;

[0010] Step S4, for the sub-waveform data separated in step S3, using an adaptive algorithm with fewer input parameters, calculate the phase relationship of each signal component, and find a signal segment that meets the polarization characteristics of Rayleigh surface waves;

[0011] Step S5, obtaining the signal component satisfying the Rayleigh wave characteristics through the method in step S4, and calculating the polarizability of different frequency bands;

[0012] Step S6: Combine steps S2-S4 to calculate the Rayleigh surface wave polarization curves recorded by the three-component ground pulsation waveforms at different scales.

[0013] Furthermore, the step S2 comprises:

[0014] Step S2.1, randomly generate N initial solutions, which represent the population size, set the variable state dimension dim of each solution, and the upper boundary L of the dimension b , the lower boundary of dimension U b , set the maximum number of iterations Max_inter, set the function for evaluating the solution fitness, and apply the fitness function of minimum envelope entropy;

[0015] Step S2.2, through the loop, calculate the fitness value of each solution in turn, and make a greedy selection, the judgment basis is the fitness function of the minimum envelope entropy;

[0016] Step S2.3, updating the state position of the population according to the current iteration state of the biological population;

[0017] Step S2.4, loop to evaluate the fitness of each solution in the new group, and update the solution and fitness value through the forward greedy selection mechanism. If the fitness of the new solution is better than the original solution, replace the original solution with the new solution and update the group fitness. At the same time, if the new solution is the current optimal solution, update the global optimal solution;

[0018] Step S2.5, through a large number of iterations, the optimal fitness value Best_R_rate, the optimal solution Best_R and the fitness value convergence of the optimal solution in each iteration are finally calculated through iteration.

[0019] Furthermore, in step S2.3, update the state position of the population, and the expression is:

[0020] M new,j = M best,j + random 1 cosγ·δ·γ·(U bj - L bj ) + L bj , random 2 < T

[0021]

[0022] where M new,j represents the new state position of the j-th dimension of the population; M best,j represents the optimal state position of the j-th dimension of the population; random 1 : a random number in the range (-1, 1), which controls the growth direction of the population; γ represents the relationship used to simulate the state position of the population and the current iteration state; δ represents the environmental factor, which simulates the influence of the external environment on the convergence of the algorithm. As the number of iterations increases, the value gradually decreases to ensure the gradual convergence of the algorithm; U bj and L bj represent the upper and lower boundaries of the j-th dimension of the population, which are used to ensure that the updated state position does not exceed the allowed range; random 2 represents a random number. If random2 < T, then update the particle position; this condition ensures a certain degree of randomness and probability during position update; T represents the convergence factor depending on the current number of iterations, which controls the possibility of update as the number of iterations increases; t represents the current calculation iteration number.

[0023] Furthermore, step S4 includes:

[0024] Step S4.1, set the initial parameters required for this calculation, mainly including the total number of data points K 0 , the number of time windows N win , the time interval tau, the maximum time window width DT max , the Nyquist frequency f nqy , the starting point f 1 and the ending point f 2 of the frequency range, and the frequency list constant constlog. Among these parameters, except for the number of time windows N win and the maximum time window width DT max which need to be defined by oneself, the rest of the parameters are time series data attributes, and the parameters follow the following calculation formulas:

[0025] tau = t 2 - t 1

[0026] f nqy =1 / (2*tau)

[0027] f start =max(f 1 ,1 / DT max )

[0028] f end =min(f 2 ,1 / DT max )

[0029] Where: f 1 、f 2 Indicates the starting and ending points of the frequency range of the original signal, and f start 、f end The start and end frequencies of the time window function.

[0030] Step S4.2: embed a Chebyshev bandpass filter Chv_fliter, whose cutoff frequency ω 0 Determined by the current frequency, the filter is as follows:

[0031]

[0032] ε represents the ripple factor, which determines the ripple size in the stop band; R s represents the stopband ripple, set to 0.5; T n Indicates the filter order, set to 2; ω 0 represents the cut-off frequency, which is [f 0 / 2 / f nyq f 0 *2 / f nyq ] loop; S indicates the default parameter;

[0033] Step S4.3: Use a Chebyshev filter to perform bandpass filtering on the signal within a certain time window, and use a cone window function to reduce the edge effect. The calculation is as follows:

[0034] S(t) 2 =Chv_filter*(taper*S(t) 1 )

[0035] Among them, S(t) 1 is the original time waveform sequence within a certain time window, S(t) 2 is the time waveform sequence after being filtered by the Chebyshev filter, and taper represents the cone window function, which is used to reduce the edge energy leakage phenomenon when the signal is truncate.

[0036] Step S4.4, calculate the zero crossing point (the zero crossing point from negative to positive) of the vertical component Vert(t) in the three-component signal, which can be calculated in the following way:

[0037]

[0038] Where derive[i] represents the zero crossing point of the vertical component signal detected by the sign change function at a certain time node;

[0039] When Vert(t) changes from negative to positive, sign(Vert(t+1)) is 1, sign(Vert(t)) is -1, the difference is 2, and dividing by 2 gives 1; when Vert(t) changes from positive to negative, sign(Vert(t+1)) is -1, sign(Vert(t)) is 1, the difference is -2, and dividing by 2 gives -1; when Vert(t) does not change its sign, the difference is 0, and dividing by 2 still gives 0.

[0040] Step S4.5, loop through all possible zero crossing points, and mark the east-west and south-north components that meet the zero crossing points at the same time, and obtain the three components E_applied(t), N_applied(t), and V_applied(t) that meet the requirements, but V_applied(t) is 1 / (4*f*tau) data points ahead of E_applied(t) and N_applied(t).

[0041] Step S4.6: Based on the three-component data in steps S4.4-S4.5, calculate the horizontal azimuth θ and correlation corr. The azimuth is calculated as follows:

[0042]

[0043] integral 1 =ΣV_applied(t).*E_applied(t)

[0044] integral 2 =ΣV_applied(t).*N_applied(t)

[0045] Among them, V_applied(t), E applied (t), N_applied(t) represent the three-component waveform data that meets the zero-crossing point detection in step S4.5; integral 1 The above V_applied(t), E applied (t) summation and integration operation; integral 2 The above V_applied(t), N applied(t) is the summation and integration operation.

[0046] Here, the azimuth needs to be converted between 0 and 360°.

[0047] The correlation value of the horizontal component and the vertical component is calculated by correlation. First, E_applied(t) and N_applied(t) need to be synthesized into H_applied(t) by horizontal components, and then the correlation coefficient corr between the horizontal components H_applied(t) and V_applied(t) is calculated, ranging from -1 to 0, as follows:

[0048] H_applied(t)=sin(Θ)*E_applied(t)+cos9Θ)*N applied (t)

[0049]

[0050] Among them, H_applied(t) represents the horizontal waveform component formed by horizontal synthesis of E_applied(t) and N_applied(t) in step S4.6; V_applied(t) is the vertical waveform component.

[0051] Furthermore, in step S5, the polarizability (H / V) of different frequency bands is calculated, and the expression is:

[0052]

[0053] Among them, H_applied(t) represents the horizontal waveform component formed by horizontal synthesis of E_applied(t) and N_applied(t) in step S4.6; V_applied(t) is the vertical waveform component.

[0054] Beneficial effects:

[0055] The present invention proposes a method for calculating the polarization curve of Rayleigh surface waves based on a joint algorithm. The present invention designs an optimization population algorithm, utilizes the forward greedy selection mechanism module of the search-migration in the algorithm, and continuously iterates to obtain the required optimal parameters. The algorithm is applied to the VMD-variational mode decomposition method (better than the filter effect), and the optimal decomposition parameters are found. Finally, the original ground pulsation record is multi-scale decomposed into sub-waveform records with high orthogonality. The waveform of each scale carries the physical property information of the geological body that does not interfere with each other, ensuring the multi-level and multi-perspective of the geological body exploration. In the method of the present invention, the propagation of Rayleigh surface waves presents a special elliptical shape, that is, the vertical vibration component has a 1 / 4 phase lead compared to the horizontal component (the body wave and the Love surface wave in the ground pulsation do not have this feature). By designing a zero-crossing phase algorithm, the data segment that meets this phenomenon is extracted; at the same time, the azimuth and correlation calculation of the horizontal and vertical component waveforms are used to ensure that the three-component waveforms are all from the Rayleigh surface wave, and finally a stable Rayleigh surface wave polarization curve is calculated. The present invention can extract stable Rayleigh surface wave components at multiple scales and levels from complex ground pulsation waveforms, complete the calculation of polarization curves, and provide technical support for subsequent multi-level and multi-perspective interpretation of the physical properties of geological bodies. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 is a flow chart of the method steps of the present invention;

[0057] Figure 2 This is an example diagram of optimization of the decomposition parameters of the signal waveform by the optimization population algorithm of the present invention;

[0058] Figure 3 A schematic diagram of signal decomposition of the optimized population algorithm of the present invention;

[0059] Figure 4 It is a schematic diagram of the polarization curve of the Rayleigh surface wave of the present invention. DETAILED DESCRIPTION

[0060] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application may be combined with each other. The present application is further described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0061] like Figure 1As shown in the figure, a Rayleigh surface wave polarization curve calculation method based on a joint algorithm is provided. The core modules of the method are mainly composed of two parts. First, the original three-component waveform data is decomposed into mutually orthogonal sub-waveforms by optimizing the population algorithm. The orthogonality ensures that each waveform carries the physical properties of geological bodies of different scales. Then, an adaptive zero-crossing point phase algorithm is used to calculate the phase relationship between the vertical component and the horizontal component of each three-component quantum waveform, extract the phase relationship satisfied by the Rayleigh surface wave, and retain all signal segments that satisfy the relationship, and finally calculate the Rayleigh surface wave polarization curves of different scales. The method includes:

[0062] Step S1, pre-processing the original waveform records of the three components of the Rayleigh surface wave polarization characteristic curve by removing instrument response, removing mean, removing potential terms, and filtering;

[0063] Step S2, using a population search-migration forward greedy selection mechanism module to find the optimal solution for parameters, providing support for subsequent decomposition of waveform data;

[0064] Step S2.2, through the loop, calculate the fitness value of each solution in turn, and perform greedy selection (that is, compare the fitness value of each solution with the current optimal value. If it is better than the current optimal value, update the optimal value and the optimal solution). The judgment here is based on the fitness function of the minimum envelope entropy.

[0065] Step S2.3: Update the state position of the population according to the current iteration state of the biological population, using the following update rules:

[0066] M new,j =M best,j +random 1 cosγ·δ·γ·(Ub j -L bj )+L bj ,random 2 <T

[0067]

[0068] Where: M new,j Represents the new state position of the jth dimension of the population. best,j Indicates the optimal state position of the j-th dimension of the population. 1 A random number in the range of (-1,1) that controls the direction of population growth. γ represents the relationship between the position of the population state and the current iteration state. δ represents the environmental factor, which simulates the influence of the external environment on the convergence of the algorithm. As the number of iterations increases, the value gradually decreases to ensure that the algorithm gradually converges. U bj With L bjDenote the upper and lower bounds of the j-th dimension of the population, which are used to ensure that the updated state position does not exceed the allowed range. random 2 Denote a random number. If random 2 < T, then update the particle position. This condition ensures a certain degree of randomness and probability during position update. T represents a convergence factor depending on the current iteration number, which controls the possibility of update as the iteration number increases. t represents the current calculation iteration number.

[0069] Specifically, for each dimension of each solution, according to the probability random 1 Judge whether to execute the population evolution search strategy. If the condition is met, update the corresponding dimension of the solution, and consider the upper and lower bounds of the solution and the current optimal solution. According to the probability random 2 Judge whether to execute the population migration mechanism. If the condition is met, directly set a certain dimension of the solution to the corresponding dimension of the current optimal solution, and perform boundary processing on each new solution to ensure that the value of the solution does not exceed the upper and lower bounds.

[0070] Step S2.4: Loop to evaluate the fitness of each solution in the new population, and update the solution and fitness value through the forward greedy selection mechanism. If the fitness of the new solution is better than the original solution, replace the original solution with the new solution and update the population fitness. At the same time, if the new solution is the current optimal solution, update the global optimal solution.

[0071] Step S2.5: Through a large number of iterations, finally iteratively calculate the optimal fitness value (Best_R_rate), the optimal solution (Best_R), and the convergence of the fitness value of the optimal solution in each iteration.

[0072] Step S3: Combine the population optimization algorithm in Step S2 with the fitness function of the minimum envelope entropy to optimize the parameters K and α in the variational mode decomposition (VMD), so that the VMD algorithm can adaptively determine the signal decomposition scale, and then use this method to decompose the original data waveform in Step S1 to obtain sub-waveform data with different orthogonality (carrying different original signal features).

[0073] The process of Step S2 is relatively abstract. Therefore Figure 2 And Figure 3 Shows an example of signal decomposition (Step S3) using Step S2. Figure 2 a is a time series signal with a sampling interval of 0.01 s. Considering that the optimized parameters in the decomposition algorithm are K (the number of decompositions) and α (the smoothness of the signal), the variable state dimension dim is 2; the upper and lower bounds L b And U bUsually it is taken as 100 and 2500, the decomposition number K can be given a range of 1 to 20, that is, the original signal can be decomposed into 1 to 20 self-signals; the maximum number of iterations Max_iter is 10; the number of populations is 10. Figure 2 b shows the randomly generated initial population range (initial K and α) and the optimal population range (optimal K and α). Figure 2 c shows the convergence curve under the number of iterations to ensure the update of the forward greedy selection mechanism. The results of this example indicate that for this signal, the optimal K value is 4, that is, the optimal sub-signals of decomposition are 4; and the optimal α value is 174, that is, the smoothness of the signal is 174.

[0074] Figure 3 a shows the waveform of the sub-signals decomposed from the signal. The waveforms of the sub-signals are different, indicating that they carry different characteristics of the original total signal. Figure 3 b shows the frequency characteristics of each sub-signal. The energy carried by each sub-signal is different, showing orthogonal incoherent characteristics.

[0075] Step S4: for the sub-signals separated in step S3, design an adaptive algorithm with fewer input parameters to calculate the phase relationship between the components of the signal, so as to find a signal segment that meets the polarization characteristics of the Rayleigh wave.

[0076] Step S4.1: Set the initial parameters required for the calculation, mainly including the total number of data points (K 0 ), the number of time windows (N win ), time interval (tau) and maximum time window width (DT max ), Nyquist frequency (f nqy ), the starting point of the frequency range (f 1 ) and end point (f 2 ) and frequency list constants (constlog). These parameters include the number of time windows (N win ) and the maximum time window width (DT max ) needs to be defined by yourself, and the other parameters are time series data attributes. The parameters follow the following calculation formula:

[0077] tau=t 2 -t 1

[0078] f nqy =1(2*tau)

[0079] f start =max(f 1 ,1 / DT max )

[0080] f end =min(f2 ,1 / DT max )

[0081] Step S4.2, embed a Chebyshev bandpass filter Chv_fliter, whose cutoff frequency ω0 is determined by the current frequency, and the filter is as follows:

[0082]

[0083] ε represents the ripple factor, which determines the ripple size in the stop band; R s Represents the stopband ripple, usually 0.5; T n Indicates the filter order, usually 2; ω 0 represents the cut-off frequency, which is [f 0 / 2 / f nyq f 0 *2 / f nyq ] loop; S indicates the default parameter;

[0084] Step S4.3: Use a Chebyshev filter to perform bandpass filtering on the signal within a certain time window, and use a cone window function to reduce the edge effect. The calculation is as follows:

[0085] S(t) 2 =Chv_filter*(taper*S(t) 1 )

[0086] Step S4.4, calculate the zero crossing point (the zero crossing point from negative to positive) of the vertical component Vert(t) in the three-component signal, which can be calculated in the following way:

[0087]

[0088] When Vert(t) changes from negative to positive, sign(Vert(t+1)) is 1, sign(Vert(t)) is -1, the difference is 2, and dividing by 2 gives 1; when Vert(t) changes from positive to negative, sign(Vert(t+1)) is -1, sign(Vert(t)) is 1, the difference is -2, and dividing by 2 gives -1; when Vert(t) does not change its sign, the difference is 0, and dividing by 2 still gives 0.

[0089] Step S4.5, loop through all possible zero crossing points, and mark the east-west and south-north components that meet the zero crossing points at the same time, and obtain the three components E_applied(t), N_applied(t), and V_applied(t) that meet the requirements, but V_applied(t) is 1 / (4*f*tau) data points ahead of E_applied(t) and N_applied(t).

[0090] Step S4.6: Based on the three-component data in steps S4.4-S4.5, calculate the horizontal azimuth (θ) and correlation (corr). The azimuth is calculated as follows:

[0091]

[0092] integral 1 =ΣV_applied(t).*E_applied(t)

[0093] integral 2 =ΣV_applied(t).*N_applied(t)

[0094] Here, the azimuth needs to be converted between 0-360°.

[0095] The correlation value of the horizontal component and the vertical component is calculated by correlation. First, E_applied(t) and N_applied(t) need to be synthesized into H_applied(t) by horizontal components, and then the correlation coefficient corr between the horizontal components H_applied(t) and V_applied(t) is calculated, ranging from -1 to 0, as follows:

[0096] H_applied(t)=sin(Θ)*E_applied(t)+cos(Θ)*N applied (t)

[0097]

[0098] Step S5: By using the method in step S4, the signal components satisfying the Rayleigh surface wave characteristics are obtained, and the polarization ratios (H / V) of different frequency bands are calculated as follows:

[0099]

[0100] Among them, H_applied(t) represents the horizontal waveform component formed by horizontal synthesis of E_applied(t) and N_applied(t) in step S4.6; V_applied(t) is the vertical waveform component.

[0101] Step S6, combining the two algorithms of steps S2 to S4 to calculate the Rayleigh surface wave polarization curves recorded by the three-component ground pulsation waveforms of different scales.

[0102] like Figure 4 As shown, the zero-crossing algorithm in steps S4 and S5 is used to extract and calculate the Rayleigh wave polarization curve. Figure 4a is a typical three-component ground pulsation original waveform record (bandpass filter 0.5-20Hz), Figure 4 b is the three-component waveform after filtering by the designed Chebyshev bandpass filter; calculate the zero crossing point of the vertical component V, loop through all possible zero crossing points, and find the data segment that satisfies the vertical component to be 1 / (4*f*tau) ahead of the horizontal two components; Figure 4 c. Figure 4 d shows the horizontal azimuth (θ) and correlation (corr) of all data segments; Figure 4 As shown in e, the data segments that meet the correlation requirements are used to produce polarization rates (H / V) of different frequency bands to form a polarization curve of a single data segment; Figure 4 As shown in Fig. 5, the polarization curves of all data segments are averaged and smoothed using a smoothing function to obtain the final Rayleigh surface wave polarization curve; Figure 4 As shown in Figure g, the presence of multiple peaks in the polarization curve indicates that the geological body has a complex structure and physical stratification at different depths.

[0103] Therefore, the population search-migration forward greedy selection mechanism module designed in step S2 is combined to optimize the decomposition function parameters, and the original three-component ground pulsation waveform is decomposed into multiple scales (step S3); then the improved zero-crossing phase analysis method designed in step S4 is used to extract the Rayleigh surface wave components for each scale waveform; and the polarization curve of the multi-scale Rayleigh surface wave is finally calculated using steps S5 and S6.

[0104] From the perspective of analytical scale, current algorithms usually filter the original ground pulsation records, and then perform an integral transformation from the time domain to the frequency domain to obtain a polarization curve. The source of the ground pulsation is background noise, which covers a complex combination of noise sources. Simple filters cannot completely separate vibrations in different frequency bands, so the polarization curve finally obtained is a corresponding aliasing of various geological bodies. For this reason, the first improved technology of the present invention is to design an optimized population algorithm, use the search-migration forward greedy selection mechanism module in the algorithm, and continuously iterate to obtain the required optimal parameters. This algorithm is applied to the VMD-variational mode decomposition method (better than the filter effect) to find the optimal decomposition parameters, and finally decompose the original ground pulsation records into sub-waveform records with high orthogonality at multiple scales. The waveforms of each scale carry non-interfering physical property information of the geological body, ensuring multi-level and multi-perspective exploration of the geological body.

[0105] Unlike the characteristic of Love surface waves that only have horizontal vibrations, the vibrations of Rayleigh surface waves exist in both horizontal and vertical components. From the perspective of Rayleigh surface wave extraction and separation, the current popular algorithm is based on the continuous wavelet transform method, which extracts the part with stronger vertical component energy. However, the key to the continuous wavelet algorithm is to select a suitable wavelet basis function. There are currently dozens of wavelet basis functions. Therefore, the data segments of Rayleigh surface waves obtained by decomposing the original ground pulsation waveform using different basis functions are also different, and the polarization curves finally obtained are also different. To this end, the second improved technology of the present invention is to design an extraction method that conforms to the propagation characteristics of Rayleigh surface waves. The propagation of Rayleigh surface waves presents a special elliptical shape, that is, the vertical vibration component has a 1 / 4 phase lead compared to the horizontal component (body waves and Love surface waves in ground pulsations do not have this feature). Therefore, a zero-crossing phase algorithm is designed to extract data segments that meet this phenomenon; at the same time, the azimuth and correlation calculation of the horizontal and vertical component waveforms are used to ensure that the three-component waveforms are all from Rayleigh surface waves, and finally a stable Rayleigh surface wave polarization curve is calculated.

[0106] The combined application of the two technical innovations of the present invention can extract stable Rayleigh surface wave components at multiple scales and levels from complex ground pulsation waveforms, complete the calculation of polarization curves, and provide technical support for subsequent multi-level and multi-perspective interpretation of the physical properties of geological bodies.

[0107] Although embodiments of the present invention have been shown and described, it will be understood by those skilled in the art that various equivalent changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for calculating Rayleigh wave polarization curve based on a joint algorithm, characterized in that , the method comprises: Step S1, pre-processing the original waveform records of the three components of the Rayleigh surface wave polarization characteristic curve by removing instrument response, removing mean, removing potential terms, and filtering; Step S2, using the population search-migration forward greedy selection mechanism module to find the optimal solution for parameters, to provide support for the subsequent decomposition of waveform data; Step S3, combining the population optimization algorithm in step S2 with the fitness function of minimum envelope entropy, optimizing the parameters in the variational mode decomposition algorithm, making the variational mode decomposition algorithm adaptively determine the signal decomposition scale, decomposing the original data waveform in step S1, and obtaining sub-waveform data of different orthogonality; Step S4, for the sub-waveform data separated in step S3, using an adaptive algorithm with fewer input parameters, calculate the phase relationship of each signal component, and find a signal segment that meets the polarization characteristics of Rayleigh surface waves; Step S5, obtaining the signal component satisfying the Rayleigh wave characteristics through the method in step S4, and calculating the polarizability of different frequency bands; Step S6: Combine steps S2-S4 to calculate the Rayleigh surface wave polarization curves recorded by the three-component ground pulsation waveforms at different scales.

2. A method for calculating a Rayleigh wave polarization curve based on a joint algorithm as claimed in claim 1, characterized in that: The step S2 comprises: Step S2.1, randomly generate N initial solutions, which represent the population size, set the variable state dimension dim of each solution, and the upper boundary L of the dimension b , the lower boundary of dimension U b , set the maximum number of iterations Max_inter, set the function for evaluating the solution fitness, and apply the fitness function of minimum envelope entropy; Step S2.2, through the loop, calculate the fitness value of each solution in turn, and make a greedy selection, the judgment basis is the fitness function of the minimum envelope entropy; Step S2.3, updating the state position of the population according to the current iteration state of the biological population; Step S2.4, loop to evaluate the fitness of each solution in the new group, and update the solution and fitness value through the forward greedy selection mechanism. If the fitness of the new solution is better than the original solution, replace the original solution with the new solution and update the group fitness. At the same time, if the new solution is the current optimal solution, update the global optimal solution; Step S2.5, through a large number of iterations, the optimal fitness value Best_R_rate, the optimal solution Best_R and the fitness value convergence of the optimal solution in each iteration are finally calculated through iteration.

3. The method for calculating the Rayleigh wave polarization curve based on the joint algorithm according to claim 1, characterized in that: The step S2.3 updates the status position of the population, and the expression is: M new,j =M best,j +random1 cosγ·δ·γ·(U bj -L bj )+L bj ,random2<T Among them, M new,j represents the new state position of the j-th dimension of the population; M best,j represents the optimal state position of the j-th dimension of the population; random1: a random number ranging from (-1, 1), controlling the growth direction of the population; γ represents the relationship used to simulate the population state position and the current iteration state; δ represents the environmental factor, simulating the influence of the external environment on the convergence of the algorithm. As the number of iterations increases, the value gradually decreases to ensure the algorithm converges step by step; U bj and L bj represent the upper and lower boundaries of the j-th dimension of the population, used to ensure that the updated state position does not exceed the permitted range; random2 represents a random number. If random2 < T, the particle position is updated; this condition ensures a certain degree of randomness and probability during position update; T represents the convergence factor dependent on the current number of iterations, controlling the possibility of update as the number of iterations increases; t represents the current calculation iteration number.

4. The method for calculating the Rayleigh wave polarization curve based on the joint algorithm according to claim 1, characterized in that: The step S4 comprises: Step S4.1, set initial parameters, including the total number of data points K0, the number of time windows N win , time interval tau maximum time window width DT max , Nyquist frequency f nqy , the starting point f1 and end point f2 of the frequency range and the frequency list constant constlog. In addition to the number of time windows N win With the maximum time window width DT max It is self-defined, and the other parameters are time series data attributes. The parameters follow the following calculation formula: tau=t2-t1 f nqy =1 / (2*year) f start =max(f1,1 / DT max ) <h2 style=";text-align:left;direction:ltr">f<h2 style=";text-align:left;direction:ltr"> end <h2 style=";text-align:left;direction:ltr"> =min(f2,1 / DT<h2 style=";text-align:left;direction:ltr"> max <h2 style=";text-align:left;direction:ltr"> ) Among them: f1, f2 represent the starting and ending points of the frequency range of the original signal, and f start 、f end The start and end frequencies of the time window function; Step S4.2: embed a Chebyshev bandpass filter Chv_fliter, whose cutoff frequency ω0 is determined by the current frequency. The filter expression is: ε represents the ripple factor, which determines the ripple size in the stop band; R s represents the stopband ripple, set to 0.5; T n represents the filter order, which is set to 2; ω0 represents the cutoff frequency, which is [f0 / 2 / f nyq f0*2 / f nyq ] loop; S indicates the default parameter; Step S4.3: Use a Chebyshev filter to perform bandpass filtering on the signal within a certain time window, and use a cone window function to reduce the edge effect. The expression is: S(t)2=Chv_filter*(taper*S(t)1) Among them, S(t)1 is the original time waveform sequence in a certain time window, S(t)2 is the time waveform sequence after filtering by Chebyshev filter, and taper represents the cone window function, which is used to reduce the edge energy leakage phenomenon when the signal is truncate; Step S4.4, calculate the zero crossing point of the vertical component Vert(t) in the three-component signal, the zero crossing point from negative to positive, and the expression is: Where derive[i] represents the zero crossing point of the vertical component signal detected by the sign change function at a certain time node; When Vert(t) changes from negative to positive, sign(Vert(t+1)) is 1, sign(Vert(t)) is -1, the difference is 2, divided by 2 is 1; when Vert(t) changes from positive to negative, sign(Vert(t+1)) is -1, sign(Vert(t)) is 1, the difference is -2, divided by 2 is -1; when Vert(t) does not change its sign, the difference is 0, divided by 2 is still 0; Step S4.5, loop through all possible zero crossing points, and mark the east-west and south-north components that meet the zero crossing points at the same time, and obtain the three components E_applied(t), N_applied(t), and V_applied(t) that meet the requirements, but V_applied(t) is 1 / (4*f*tau) data points ahead of E_applied(t) and N_applied(t); Step S4.6: Based on the three-component data in steps S4.4-S4.5, calculate the horizontal azimuth θ and correlation corr. The azimuth calculation expression is: integral1=∑V_applied(t).*E_applied(t) integral2=∑V_applied(t).*N_applied(t) Among them, V_applied(t), E applied (t), N_applied(t) represent the three-component waveform data that meets the zero-crossing point detection in step S4.5; integral1 represents V_applied(t), E applied (t); integral2 represents V_applied(t), N applied (t) summation and integration operation; azimuth angle conversion between 0 and 360°; The correlation value of the horizontal component and the vertical component is calculated by correlation, and E_applied(t) and N_applied(t) are synthesized into H_applied(t) for horizontal components. Then, the correlation coefficient corr between the horizontal components H_applied(t) and V_applied(t) is calculated, ranging from -1 to 0, and the expression is: H_applied(t)=sin(Θ)*E_applied(t)+cos(Θ)*N applied (t) (15) Among them, H_applied(t) represents the horizontal waveform component formed by horizontal synthesis of E_applied(t) and N_applied(t) in step S4.6; V_applied(t) is the vertical waveform component.

5. The method for calculating the Rayleigh wave polarization curve based on the joint algorithm according to claim 1, characterized in that: In step S5, the polarizability H / V of different frequency bands is calculated, and the expression is: Among them, H_applied(t) represents the horizontal waveform component formed by horizontal synthesis of E_applied(t) and N_applied(t) in step S4.6; V_applied(t) is the vertical waveform component.

Images