A calculation method for Rayleigh surface wave polarization curves based on a joint algorithm

Abstract

The present invention discloses a method for calculating the Rayleigh surface wave polarization curve based on a combined algorithm, and the method includes: preprocessing the three-component original waveform record; using a forward greedy selection mechanism module of population search - migration to find the optimal solution of parameters; optimizing the parameters in the variational mode decomposition algorithm to enable the variational mode decomposition algorithm to adaptively determine the signal decomposition scale and obtain sub-waveform data with different orthogonality; using an adaptive algorithm with fewer input parameters to calculate the phase relationship of each component of the signal and find the signal segment that satisfies the Rayleigh surface wave polarization characteristics; calculating the polarization rate of different frequency bands; and calculating the Rayleigh surface wave polarization curve of the three-component ground pulsation waveform record at different scales. The present invention can extract stable Rayleigh surface wave components from complex ground pulsation waveforms at multiple scales and levels, complete the calculation of the polarization curve, and provide technical support for subsequent multi-level and multi-perspective interpretation of the physical properties of geological bodies.

Description

Technical Field

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

[0002] With the rapid development of urbanization, the demand for underground space resources is increasing day by day. The high-resolution and high-precision exploration of the near-surface geological structure has become the safety foundation for engineering construction and operation. As a relatively safe and convenient method for non-destructive geological exploration, ground motion monitoring has been increasingly widely used in the fields of urban underground space exploration, geological disaster body exploration, and seismic site classification. Ground motion, also known as ambient vibration (in the Japanese region) or background noise (in the European and American regions), is a wave field with rich environmental seismic noise as the seismic source. Its main seismic phase is the surface wave (Rayleigh surface wave and Love surface wave) propagating along the near surface. Among them, the Rayleigh surface wave is more widely used in multi-scale geological-geophysical exploration due to its dispersion characteristics and polarization characteristics. The polarization characteristics of the Rayleigh surface wave refer to the ratio of its horizontal vibration component (H) to the vertical vibration component (V) at different frequencies, so it is reflected in the form of a frequency (X-axis)-polarization ratio-H / V (Y-axis) curve. The morphological changes of the polarization curve are highly correlated with the formation elastic wave velocity, density, layer thickness, and Poisson's ratio. Therefore, how to extract the Rayleigh surface wave in the in-situ monitored ground motion at multiple levels and calculate its polarization characteristic curve (H / V curve) is extremely important for inverting the optimal physical parameter model of the geological body (elastic wave velocity, density, layer thickness, etc.) and finally interpreting the characteristics of the geological body.

[0003] At present, the core of calculating the polarization characteristic curve of Rayleigh surface waves is to transform the observed vibration waveforms (time domain) of the original three components (east-west / north-south / vertical) into the spectral amplitude curves (frequency domain) of the three components through numerical integral transformation. Then, the amplitude curves of the two horizontal components are synthesized into a horizontal amplitude curve, and finally, the ratio is taken with the vertical amplitude curve to obtain the polarization characteristic curve. The X-axis of this curve is the frequency, and the Y-axis is the H / V value (polarization rate). The peak of the polarization curve represents obvious changes in physical parameters in the geological body, such as inherent properties like wave velocity and stiffness. However, the deficiencies of the current algorithm mainly include the following two aspects: 1. In terms of enhancing the Rayleigh surface wave component, the current continuous wavelet transform algorithm is adopted. Although this algorithm can take into account the signal energy resolution in both the time domain and the frequency domain, there are a variety of wavelet basis functions used in the algorithm, and the parameters of the wavelet basis function need to be set artificially. Therefore, the semi-empirical setting of the algorithm parameters has a certain impact on the calculation results of the final polarization curve. 2. The natural geological body scale has different responses to the vibration wave fields of noise sources in different frequency bands. Therefore, the surface wave polarization curves presented in different regions of the geological body are of different morphologies. Currently, the extraction and calculation of polarization curves are all based on aliased waveform signals. Although the interference signals in specific frequency bands can be separated through the filtering module in the algorithm, the remaining signal energy still cannot exist orthogonally. Therefore, the polarization curve also shows a mixed state of various morphologies and cannot be separated into multi-scale polarization curves to reflect the physical property characteristics of multi-scale geological bodies. Summary of the Invention

[0004] In order to overcome the disadvantages and deficiencies of the existing technology, the present invention provides a method for calculating the Rayleigh surface wave polarization curve based on a joint algorithm. This invention mainly improves the extraction technology of the Rayleigh surface wave polarization curve for the current ground motion source. Aiming at the traditional extraction method of performing overall numerical integral transformation on in-situ observation data, the objectives of this invention are as follows: 1. By optimizing the population algorithm, the in-situ observation waveform data is separated into multi-scale independent orthogonal time series components to minimize signal correlation; 2. Through the improved zero-crossing phase analysis method, the waveform time series of different scales are calculated to obtain the 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 the physical characteristic inversion of the geological structure with multi-resolution, reduce the inversion error rate, and thus more comprehensively and accurately interpret the structural characteristics of the geological body.

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

[0007] Step S1, perform pre-processing on the three-component original waveform records of the Rayleigh surface wave polarization characteristic curve, including removing instrument response, removing the mean, removing the trend term, and filtering.

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

[0009] Step S3: Combine the population optimization algorithm in Step S2 with the fitness function of minimum envelope entropy to optimize the parameters K and α in the variational mode decomposition algorithm, enabling the variational mode decomposition algorithm to adaptively determine the signal decomposition scale, decompose the original data waveform in Step S1, and obtain sub - waveform data with different orthogonality;

[0010] Step S4: For the sub - waveform data separated in Step S3, use an adaptive algorithm with fewer input parameters to calculate the phase relationship of each signal component and find the signal segment that satisfies the polarization characteristics of Rayleigh surface waves;

[0011] Step S5: Through the method in Step S4, obtain the signal components that satisfy the Rayleigh surface wave characteristics and calculate the polarization rates in different frequency bands;

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

[0013] Furthermore, Step S2 includes:

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

[0015] Step S2.2: Through a loop, calculate the fitness value of each solution in turn and perform greedy selection, with the judgment basis being the fitness function of minimum envelope entropy;

[0016] Step S2.3: Update 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 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;

[0018] 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.

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

[0020] M new,j = M best,j + random1cosγ·δ·γ·(U bj - L bj ) + L bj , random2 < 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; random1: a random number ranging from (-1, 1), controlling 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, 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 gradual convergence of the algorithm; 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 allowed range; random2 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, controlling the possibility of update as the number of iterations increases; t represents the current calculation iteration number.

[0023] Further, step S4 includes:

[0024] Step S4.1: Set the initial parameters required for this calculation, mainly including the total number of data points K0, 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 f1 and the ending point f2 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 the attributes of the time series data, and the parameters follow the following calculation formulas:

[0025] tau = t2 - t1

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

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

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

[0029] Where: f1 and f2 represent the starting and ending points of the frequency range of the original signal, and f start , f end are the starting and cut-off frequencies of the time window function.

[0030] Step S4.2: Embed a Chebyshev band-pass filter Chv_fliter, whose cut-off frequency ω0 is determined by the current frequency. The filter is as follows:

[0031]

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

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

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

[0035] Where, 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 taper window function, which is used to reduce the edge energy leakage phenomenon when the signal is truncated.

[0036] Step S4.4: Calculate the zero crossing points (the zero crossings 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 sign change function at a certain time node to detect the zero crossing points of the vertical component signal;

[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 there is no sign change in Vert(t), the difference is 0, and dividing by 2 is still 0.

[0040] Step S4.5: Traverse all possible zero-crossing points in a loop, and simultaneously mark the east-west and south-north components that meet the zero-crossing points to obtain the three components E_applied(t), N_applied(t), V_applied(t) that meet the requirements. However, 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] integral1 = ΣV_applied(t).*E_applied(t)

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

[0045] where V_applied(t), E applied (t), N_applied(t) represent the three-component waveform data that meet the zero-crossing detection in Step S4.5; integral1 represents the summation integral operation of the above V_applied(t), E applied (t); integral2 represents the summation integral operation of the above V_applied(t), N applied (t).

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

[0047] To calculate the correlation value between the horizontal component and the vertical component through correlation, first, the horizontal components E_applied(t) and N_applied(t) need to be combined into H_applied(t), and then the correlation coefficient corr between the horizontal component H_applied(t) and V_applied(t) is calculated, which ranges from -1 to 0. The calculation is as follows:

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

[0049]

[0050] Among them, H_applied(t) represents the horizontal waveform component formed by horizontally synthesizing 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 horizontally synthesizing 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 Rayleigh surface wave polarization curve based on a joint algorithm. The present invention designs an optimized population algorithm, and uses the forward greedy selection mechanism module of search - migration in the algorithm to continuously iterate to obtain the required optimal parameters. This algorithm is applied to the VMD - variational mode decomposition method (superior to the filter effect) to find the best decomposition parameters. Finally, the original ground motion records are multi - scale decomposed into sub - waveform records with high orthogonality. The waveforms of each scale carry the physical property information of geological bodies that do not interfere with each other, ensuring multi - level and multi - perspective exploration of geological bodies. 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 (body waves and Love surface waves in ground motion do not have this characteristic). By designing a zero - crossing phase algorithm, the data segments that meet this phenomenon are extracted; at the same time, by calculating the azimuth angle and correlation of the horizontal and vertical component waveforms, it is ensured that the three - component waveforms all come from Rayleigh surface waves, and finally a stable Rayleigh surface wave polarization curve is calculated. The present invention can multi - scale and multi - level extract stable Rayleigh surface wave components from complex ground motion waveforms, complete the calculation of the polarization curve, and provide technical support for subsequent multi - level and multi - perspective interpretation of the physical properties of geological bodies. Description of the drawings

[0056] Figure 1 is the flowchart of the method steps of the present invention;

[0057] Figure 2Optimization example diagram of the decomposition parameters of the signal waveform by the optimized population algorithm of the present invention;

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

[0059] Figure 4 Schematic diagram of the polarization curve of Rayleigh surface waves of the present invention. Detailed implementation manners

[0060] It should be noted that, without conflict, the embodiments in this application and the features in the embodiments can be combined with each other. The following further describes this application in detail with reference to the drawings and specific embodiments.

[0061] As Figure 1 shown, a method for calculating the polarization curve of Rayleigh surface waves based on a joint algorithm. The core module of this method mainly consists of two parts. First, through the optimized population algorithm, the original three-component waveform data is decomposed into mutually orthogonal sub-waveforms, and the orthogonality ensures that each waveform carries the physical properties of geological bodies at 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 sub-waveform, extract the phase relationship satisfied by the Rayleigh surface waves, and retain all signal segments that satisfy the relationship. Finally, the polarization curves of Rayleigh surface waves at different scales are calculated. This method includes:

[0062] Step S1: Perform pre-processing on the three-component original waveform record of the Rayleigh surface wave polarization characteristic curve, including removing instrument response, removing the mean value, removing the trend term, and filtering.

[0063] Step S2: Use the forward greedy selection mechanism module of population search - migration to find the optimal solution of parameters and provide support for subsequent decomposition of waveform data;

[0064] Step S2.2: Through 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 basis here is through 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 rule:

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

[0067]

[0068] 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. random1 is a random number in the range (-1, 1), which controls 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, 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 allowable range. random2 represents a random number. If random2 < T, the particle position is updated. This condition ensures a certain degree of randomness and probability in position updating. T represents the convergence factor depending on the current number of iterations, which controls the possibility of updating as the number of iterations increases. t represents the current calculation iteration number.

[0069] Specifically, for each dimension of each solution, it is judged whether to execute the population evolution search strategy according to the probability random1. If the condition is met, the corresponding dimension of the solution is updated, and the upper and lower bounds of the solution and the current optimal solution are considered. It is judged whether to execute the population migration mechanism according to the probability random2. If the condition is met, a certain dimension of the solution is directly set to the corresponding dimension of the current optimal solution, and boundary processing is performed on each new solution to ensure that the value of the solution does not exceed the upper and lower boundaries.

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

[0071] Step S2.5: Through a large number of iterations, finally 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 characteristics).

[0073] The process of Step S2 is relatively abstract, so Figure 2 and Figure 3Shows an example of signal decomposition using step S2 (step S3). Figure 2 a is a time series signal with a sampling interval of 0.01 s. Considering that the optimization 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 lower and upper bounds L b and U b are usually taken as 100 and 2500, the number of decompositions K can be given a range from 1 to 20, that is, the original signal can be decomposed into 1 to 20 intrinsic signals; the maximum number of iterations Max_iter is taken as 10; the number of populations is taken as 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 best number of decomposed sub-signals is 4; and the best α value is 174, that is, the smoothness of the signal is 174.

[0074] Figure 3 a shows the waveform of the sub-signals obtained by decomposing this signal. The waveforms of each sub-signal are different, indicating that they carry different characteristics of the original total signal, while Figure 3 b shows the frequency characteristics of each sub-signal. The energies carried by each sub-signal are different, presenting orthogonal and 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 of each component of the signal, so as to find the signal segment that satisfies the polarization characteristics of Rayleigh surface waves.

[0076] Step S4.1: Set the initial parameters required for this calculation, mainly including the total number of data points (K0), 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 (f1) and ending point (f2) 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 the attributes of time series data, and the parameters follow the following calculation formulas:

[0077] tau = t2 - t1

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

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

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

[0081] Step S4.2: Embed a Chebyshev band-pass filter Chv_fliter, whose cut-off frequency ω0 is determined by the current frequency. The filter is as follows:

[0082]

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

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

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

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

[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) has no sign change, the difference is 0, and dividing by 2 is still 0.

[0089] Step S4.5: Loop through all possible zero-crossing points, and at the same time mark the east-west and south-north components that meet the zero-crossing points to obtain the three-component E_applied(t), N_applied(t), V_applied(t) that meet the requirements. However, 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 calculation is as follows:

[0091]

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

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

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

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

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

[0097]

[0098] Step S5: Through the method in Step S4, obtain the signal components that meet the characteristics of Rayleigh surface waves, and calculate the polarizability (H / V) of different frequency bands. The calculation is as follows:

[0099]

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

[0101] Step S6: Combine the two algorithms in Steps S2 to S4 to calculate the Rayleigh surface wave polarization curve of the three-component ground motion waveforms recorded at different scales.

[0102] As Figure 4 shown, it shows the extraction and calculation of the Rayleigh surface wave polarization curve using the zero-crossing algorithm in Steps S4 and S5. Figure 4 a is the original waveform record of typical three-component ground motion (band-pass filtered at 0.5 - 20 Hz). Figure 4b is the three-component waveform after filtering by the designed Chebyshev band-pass filter; calculate the zero-crossing points of the vertical component V, loop through all possible zero-crossing points, and find the data segment where the vertical component is 1 / (4*f*tau) ahead of the two horizontal components; as Figure 4 c, Figure 4 d shows the horizontal azimuth (Θ) and correlation (corr) of all data segments; as Figure 4 e shows that for the data segments meeting the correlation requirements, the polarizability (H / V) of different frequency bands is calculated to form the polarization curve of a single data segment; as Figure 4 f shows that the polarization curves of all data segments are averaged and smoothed using a smoothing function to obtain the final Rayleigh surface wave polarization curve; as Figure 4 g shows that multiple peaks in the polarization curve represent that the geological body structure is relatively complex and there are physical stratifications at different depths.

[0103] Therefore, the forward greedy selection mechanism module of population search - migration designed in step S2 is used to optimize the parameters of the decomposition function, and the original three-component ground pulsation waveform is decomposed at 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; finally, the polarization curves of the multi-scale Rayleigh surface waves are calculated using steps S5 and S6.

[0104] In terms of the analysis scale level, the current algorithms usually filter the original ground pulsation records and then perform the integral transformation from the time domain to the frequency domain to obtain the polarization curve. The source of the ground pulsation is background noise, covering a complex combination of noise sources. Simple filters cannot completely separate the vibrations of different frequency bands, so the final obtained polarization curve is a situation where the responses of various geological bodies are aliased. For this reason, the first improved technique of the present invention is to design an optimized population algorithm. Using the search - migration forward greedy selection mechanism module in the algorithm, the required optimal parameters are continuously iteratively obtained. This algorithm is applied to the VMD - variational mode decomposition method (superior to the filter effect) to find the best decomposition parameters. Finally, the original ground pulsation record is decomposed into sub-waveform records with high orthogonality at multiple scales. The waveforms of each scale carry the physical property information of geological bodies that do not interfere with each other, ensuring a multi-level and multi-perspective exploration of geological bodies.

[0105] Different from the Love surface wave which only has horizontal vibrations, the vibrations of Rayleigh surface wave exist in both horizontal and vertical components. In terms of the extraction and separation of Rayleigh surface wave, the current popular algorithm is based on the continuous wavelet transform method, and 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 dozens of current wavelet basis functions, so the data segments of Rayleigh surface wave finally decomposed from the original ground pulsation waveform are different for different basis functions, and the final polarization curve shapes obtained are also different. Therefore, the second improved technology of the present invention is to design an extraction method that conforms to the propagation characteristics of Rayleigh surface wave. The propagation of Rayleigh surface wave presents a special elliptical shape, that is, the vertical vibration component has a 1 / 4 phase lead compared with the horizontal component (body wave and Love surface wave in ground pulsation do not have this characteristic). Therefore, a zero-crossing phase algorithm is designed to extract the data segment that satisfies this phenomenon; at the same time, the azimuth angle and correlation of the horizontal and vertical component waveforms are calculated to ensure that the three-component waveforms all come from Rayleigh surface wave, and finally a stable Rayleigh surface wave polarization curve is calculated.

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

[0107] Although the embodiments of the present invention have been shown and described, for those of ordinary skill in the art, it can be understood that various equivalent changes, modifications, substitutions and variations can be made to these embodiments without departing from the principle and spirit of the present invention. The scope of the present invention is defined by the appended claims and their equivalent scope.

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. A method for calculating Rayleigh wave polarization curve based on a joint algorithm as claimed in claim 2, 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 in the range (-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 allowed 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 depending 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. A method for calculating Rayleigh wave polarization curve based on a joint algorithm as claimed in claim 4, 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