Phase velocity imaging method based on Rayleigh wave time-frequency characteristic similarity analysis

Abstract

The invention discloses a phase velocity imaging method based on Rayleigh wave time-frequency characteristic similarity analysis, and the method comprises the steps: firstly, designing a reasonable working frequency band, and collecting a seismic signal; secondly, de-noising the seismic original signal by using wavelet transform; then, intercepting Rayleigh waves in a time domain; thirdly, obtaining the time waveform of each frequency component of each Rayleigh wave through time-frequency analysis, calculating the time difference through similarity analysis, and obtaining a frequency dispersion curve between two channels; and finally, in combination with a plurality of Rayleigh wave frequency dispersion curves, obtaining a phase velocity imaging graph below the seismic arrangement. The method is good in effect when applied to engineering investigation.

Description

Technical Field

[0001] The present invention relates to the research field of a Rayleigh wave exploration data processing method, in particular to a phase velocity imaging method based on similarity analysis of Rayleigh wave time-frequency characteristics. Background Art

[0002] Rayleigh waves are a type of seismic wave generated by the interaction of the vertical components of longitudinal and transverse waves on a free interface. They are characterized by low frequency, low velocity, and high energy. Since the mid-19th century, seismic surface wave exploration technology has gradually become an important means of studying information about near-surface subsurface media. Rayleigh (1885) first discovered Rayleigh waves and discussed their mathematical mechanisms in detail. Subsequent large-scale research has demonstrated that Rayleigh waves have a dispersion characteristic (the propagation velocity of Rayleigh waves in multi-layer media varies with frequency). Rayleigh waves of different frequencies penetrate to different depths. Based on this characteristic, the Rayleigh wave dispersion curve can be obtained to calculate the Rayleigh wave phase velocity of strata at different depths, thereby achieving the purpose of stratigraphic division.

[0003] Rayleigh wave exploration is widely used in various engineering fields. Current Rayleigh wave data processing methods, including phase difference and time-frequency characteristic peak detection, all present several challenges. First, Rayleigh wave data processing methods based on phase difference have low utilization of Rayleigh wave signal information, utilizing only the frequency domain information of a single mode. Second, Rayleigh wave data processing methods based on time-frequency characteristic peak detection often use peak detection to calculate dispersion curves after performing time-frequency analysis of the Rayleigh waves, resulting in poor stability. Third, the lateral resolution is insufficient; that is, traditional methods can only provide the average change in Rayleigh wave phase velocity beneath the seismic array. The complex and highly heterogeneous near-surface stratigraphic structure, coupled with the impact of human engineering activities such as underground pipelines and tunnels, pose challenges to the accuracy of Rayleigh wave exploration. Summary of the Invention

[0004] To this end, the present invention provides a phase velocity imaging method based on similarity analysis of Rayleigh wave time-frequency characteristics. Time-frequency analysis can be used to determine the relationship between Rayleigh wave time, frequency, and amplitude, simultaneously acquiring Rayleigh wave time and frequency domain information, improving information utilization. Similarity analysis can be used to obtain the time difference between traces of the same frequency component, thereby obtaining the dispersion curve between the two traces. This improves the stability of dispersion curve calculations in the time-frequency characteristic peak detection method. Combining the calculated results of multiple Rayleigh wave dispersion curves, a Rayleigh wave phase velocity image can be obtained from a single seismic array, clarifying the three-dimensional spatial position of geological anomalies and significantly improving lateral resolution.

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

[0006] A phase velocity imaging method based on similarity analysis of Rayleigh wave time-frequency characteristics comprises the following steps:

[0007] Step 1: Rayleigh wave exploration data collection:

[0008] The Rayleigh wave exploration frequency band plays an important role in physical filtering during field data collection and digital filtering during data processing. The wavelength of a Rayleigh wave at a certain frequency is proportional to the detection depth. The relationship between the Rayleigh wave phase velocity, frequency, and detection depth can be expressed as:

[0009]

[0010] β is the depth coefficient, which is a constant related to the Poisson's ratio of the rock and soil (as shown in Table 1).

[0011] Table 1 Relationship between Poisson's ratio and depth coefficient

[0012]

[0013] A series of geophones are placed simultaneously on pre-selected survey lines and points. An 18-pound sledgehammer is used to strike the ground to stimulate seismic waves. The geophones receive the waves, and the host computer records the waveforms. The selection of the data acquisition frequency band depends mainly on the required depth of detection and shallow resolution of the project. The maximum depth of detection is h max Determines the lower limit of the working frequency band (lowest frequency) and the minimum detection depth (or detection resolution) h min Determine the upper limit (maximum frequency) of the operating frequency band. Then, the maximum operating frequency is The minimum operating frequency is

[0014] Step 2: Waveform denoising:

[0015] The Rayleigh wave exploration work involved in the present invention can be carried out in cities or in the wild. If the exploration work is located in the city, it will be seriously disturbed by vibrations caused by various urban construction projects, vehicle types, etc. Therefore, it is considered to use the wavelet threshold filtering method to process the original signal (if the exploration work is located in the wild, denoising is not required). The process of wavelet threshold denoising is divided into three parts: wavelet decomposition, threshold processing and wavelet reconstruction. When decomposing the wavelet, the choice of the number of layers is crucial: the larger the number of layers, the more obvious the different characteristics of the noise and signal performance, and the more conducive to denoising. But at the same time, the reconstructed signal distortion will also be greater, which will in turn affect the denoising effect. The frequency band range of wavelet decomposition is related to the sampling frequency. If N layers of decomposition are used, the size of each frequency band is Fs / 2 / 2^N.

[0016] Step 2: Rayleigh wave separation:

[0017] After filtering, the Rayleigh waves need to be separated from the seismic waveform. This step can be performed directly on the time-domain seismic waveform. Compared to longitudinal and shear waves, Rayleigh waves are slower and have larger amplitudes. Based on this, each Rayleigh wave can be truncated. After truncating, the first point of the first Rayleigh wave waveform is used as the initial time, and the last point of the last waveform is used as the ending time. Missing parts of each waveform are padded with zeros to ensure that each waveform has the same length, facilitating subsequent calculations.

[0018] Step 4: Rayleigh wave time-frequency analysis:

[0019] The time-frequency analysis expression can be written as:

[0020] D(t,f,A)=TFA[g(t)] (6)

[0021] Where TFA is a time-frequency analysis method, D(t,f,A) is the time-frequency analysis result, and g(t) is the Rayleigh wave signal.

[0022] The Rayleigh wave is processed using the time-frequency analysis method to obtain the relationship between the time, frequency and amplitude of the Rayleigh wave.

[0023] Step 5: Calculate the Rayleigh wave time difference:

[0024] Set the amplitude threshold A0 and only analyze the frequency components with amplitudes greater than A0. Assume that the peak value of the amplitude of the time-frequency distribution of a waveform is A max , then A0=0.02A max For the nth and n+1th Rayleigh wave time-frequency distributions, extract the component with frequency f and obtain the time domain waveform of this frequency component, which are recorded as time series y1(x,f) and y2(x,f), respectively. Assuming that the number of signal sampling points is N0, the component y1(x,f) with frequency f is translated m(f) times in the direction of increasing time, with a translation step of t0, and the root mean square error is used to evaluate the similarity.

[0025]

[0026] When the MSE takes the minimum value, it means that the similarity between y1(x,f) and y2(x,f) is the highest. If the corresponding number of translations is m0(f), it can be considered that the component with frequency f in the Rayleigh wave signal has an inter-channel time difference t between the nth channel and the n+1th channel. d =m0(f)t0.

[0027] Step 6: Calculate the dispersion curve:

[0028] After obtaining the inter-channel time difference of each frequency component, the Rayleigh wave phase velocity V of each frequency component is obtained. R It can be expressed as:

[0029]

[0030] Where ΔX is the track spacing.

[0031] Using formula (5), the Rayleigh wave phase velocity-depth curve of the rock and soil between two adjacent detectors can be obtained.

[0032] Step 7: Draw the phase velocity imaging diagram:

[0033] Combining the phase velocity-depth curves calculated by pairwise combination of multi-channel geophones, the Rayleigh wave phase velocity-depth profile beneath the earthquake array can be obtained.

[0034] Compared with the prior art, the present invention has the following advantages:

[0035] (1) Based on time-frequency analysis, the time-frequency characteristics of Rayleigh waves can be obtained. In the subsequent analysis process, on the one hand, the time domain and frequency domain information of Rayleigh waves can be used simultaneously, and on the other hand, the information of Rayleigh waves of different modes can be used simultaneously, thereby greatly improving the utilization rate of Rayleigh wave information.

[0036] (2) When using time-frequency analysis to obtain the time difference between two adjacent channels with different frequency components, the similarity analysis method calculates the time difference between channels by the similarity of the time waveforms of the frequency components of the Rayleigh waves in each channel. The calculation results have better accuracy and stability.

[0037] (3) Combining the calculation results of multi-channel Rayleigh wave dispersion curves and using a seismic array, the phase velocity profile of the stratum below the seismic array can be obtained, which greatly improves the lateral resolution of Rayleigh wave detection. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 A flow chart of the steps required for the present invention;

[0039] Figure 2 This is a seismic exploration data collection set according to an embodiment of the present invention;

[0040] Figure 3 This is a waveform diagram of the original seismic wave according to an embodiment of the present invention;

[0041] Figure 4 This is a comparison diagram of the waveform before and after wavelet transform denoising according to an embodiment of the present invention;

[0042] Figure 5 Schematic diagram of Rayleigh wave separation according to an embodiment of the present invention;

[0043] Figure 6 This is the Choi-Williams time-frequency distribution diagram of the Rayleigh wave in the embodiment of the present invention

[0044] Figure 7The waveform diagram of the component with frequency f is extracted from the nth and n+1th Rayleigh wave Choi-Williams time-frequency distributions respectively in the embodiment of the present invention;

[0045] Figure 8 Schematic diagram of the translation of two adjacent Rayleigh wave components with a frequency of f according to an embodiment of the present invention;

[0046] Figure 9 This is a Rayleigh wave phase velocity imaging diagram below the earthquake arrangement according to an embodiment of the present invention;

[0047] Figure 10 The data processing results are compared between the phase velocity imaging method based on the similarity analysis of the time-frequency characteristics of Rayleigh waves described in the present invention and the phase difference method, one of the current Rayleigh wave data processing methods.

[0048] Figure 11 This is a comparison diagram of the stability of obtaining the Rayleigh wave time difference between the phase velocity imaging method based on the similarity analysis of the Rayleigh wave time-frequency characteristics of the present invention and the time-frequency characteristic spectrum peak detection method, which is one of the current Rayleigh wave data processing methods.

[0049] Figure 12 This is a phase velocity imaging diagram of Rayleigh wave exploration applied to landslide sliding surface detection in the present invention.

[0050] Figure 13 This is a phase velocity imaging diagram of Rayleigh wave exploration applied to the detection of permafrost active layers according to the present invention.

[0051] Figure 14 This is a phase velocity imaging diagram of Rayleigh wave exploration applied to obstacle detection in pipe jacking construction according to the present invention. DETAILED DESCRIPTION

[0052] To help those skilled in the art better understand the technical solutions of this specification, the following detailed and complete description of the embodiments of this specification is provided with reference to the accompanying drawings. Obviously, the embodiments described are only a portion of the embodiments in this specification, not all of them. All other embodiments derived by those skilled in the art based on one or more embodiments in this specification without inventive effort should fall within the scope of protection of the embodiments of this specification.

[0053] The present invention will be described in detail below with reference to the accompanying drawings:

[0054] The method comprises the following steps:

[0055] Step 1: Rayleigh wave data collection.

[0056] In a certain city, the early detection of civil air defense projects is carried out with a planned detection depth of 20m. Figure 2As shown, an 18-pound sledgehammer was used as the excitation source, and 24 4.5Hz receivers were used with a channel spacing of 0.3m. From the surface to a depth of 20m, the strata consisted of Quaternary overburden and heavily weathered sandstone. Considering the site's rock and soil conditions, the Poisson's ratio was approximately 0.35. As shown in Table 1, β was chosen as 0.75. Field testing revealed that the Rayleigh wave velocity in the shallow overburden was approximately 450m / s, while the average Rayleigh wave velocity in the heavily weathered sandstone within the detection depth range was approximately 1200m / s, resulting in an operating frequency band of 37.5Hz-337.5Hz for this detection.

[0057] Step 2: Waveform denoising.

[0058] Figure 3 This is the original waveform of a shallow earthquake survey conducted along a main road during a city civil air defense project. The waveform contains a significant amount of high-frequency, low-amplitude noise, which severely impacts the subsequent calculation of Rayleigh wave velocity. This paper considers processing the original signal using wavelet threshold filtering.

[0059] For the signal g(t), the expression of its wavelet transform can be written as:

[0060]

[0061] Where a is the scale parameter, and a>0; b is the translation parameter. ψ(t) is a given function, called the mother wavelet or basic wavelet.

[0062] set up

[0063] ψ a,b (t) is called the wavelet basis function. Then, the wavelet transform of g(t) can also be written as:

[0064]

[0065] The process of wavelet threshold denoising is divided into three parts: wavelet decomposition, threshold processing and wavelet reconstruction.

[0066] Specifically, the number of layers in wavelet decomposition is crucial. The greater the number of layers, the more pronounced the differences between noise and signal, facilitating denoising. However, this also increases the distortion of the reconstructed signal, which in turn affects the denoising effect. The frequency band of wavelet decomposition is related to the sampling frequency. For an N-layer decomposition, the size of each frequency band is Fs / 2 / 2^N. The raw seismic exploration signal collected in this paper has a sampling interval of 0.0001s, a sampling frequency of 10kHz, 3400 sampling points, and a total signal duration of 0.34s. The sampling theorem indicates that the maximum frequency of this signal is 5kHz. Therefore, a three-layer wavelet decomposition of this signal yields a frequency band of 2.5-5kHz for first-order details and a frequency band of less than 2.5kHz for first-order approximation. A frequency band of 1.25-2.5kHz for second-order details and a frequency band of less than 1.25kHz for approximation. A frequency band of approximately 0.625-1.25kHz for third-order details and a frequency band of less than 0.625kHz for approximation.

[0067] In the wavelet domain, the effective signal coefficient is large, while the noise coefficient is small. Set the coefficient threshold to:

[0068]

[0069] Where N is the signal length and σ is the noise variance. When the absolute value of the wavelet coefficient is less than the threshold, it is set to 0; when it is greater than the threshold, it remains unchanged. After thresholding, wavelet reconstruction is performed to restore the signal. Figure 4 Comparison of seismic waveforms before and after denoising.

[0070] Step 3: Rayleigh wave separation.

[0071] This step can be done directly on the time domain seismic waveform (e.g. Figure 5 (as shown). Compared to longitudinal and transverse waves, Rayleigh waves have slower speeds and larger amplitudes. Based on this, each Rayleigh wave can be truncated. After truncating, the first point of the first Rayleigh wave waveform is used as the initial time, and the last point of the last waveform is used as the ending time. Missing parts of each waveform are padded with zeros to ensure that each waveform has the same length, facilitating subsequent calculations.

[0072] Step 4: Rayleigh wave time-frequency analysis.

[0073] The Choi-Williams distribution, one of the time-frequency analysis methods, is used to process Rayleigh waves. The Choi-Williams distribution expression can be written as:

[0074]

[0075] Where, CWD g (t,f) is the Choi-Williams distribution, and σ is the kernel function parameter.

[0076] The relationship between the time, frequency and amplitude of Rayleigh waves can be obtained by using the Choi-Williams distribution of Rayleigh waves (such as Figure 6 shown).

[0077] Step 5: Obtain the Rayleigh wave time difference.

[0078] by Figure 6 Take the nth and n+1th Rayleigh waves in the example. Set the amplitude threshold A0 and only analyze the frequency components with amplitudes greater than A0. Assume that the peak value of the time-frequency distribution amplitude of a waveform is A max , then A0=0.02A max For the nth and n+1th Rayleigh wave time-frequency distributions, extract the component with frequency f (e.g. Figure 7 As shown in the figure), the time domain waveform of the frequency component is obtained, which is recorded as time series y1(x,f) and y2(x,f). Assuming that the number of signal sampling points is N0, the component y1(x,f) with frequency component f is shifted m(f) times along the direction of increasing time (as shown in the figure). Figure 8 As shown), the translation step is t0, and the root mean square error is used to evaluate the similarity.

[0079]

[0080] When the MSE takes the minimum value, it means that the similarity between y1(x,f) and y2(x,f) is the highest. If the corresponding number of translations is m0(f), it can be considered that the component with frequency f in the Rayleigh wave signal has an inter-channel time difference t between the nth channel and the n+1th channel. d =m0(f)t0.

[0081] Step 6: Dispersion Curve Calculation

[0082] After obtaining the inter-channel time difference of each frequency component, the Rayleigh wave phase velocity V of each frequency component is obtained. R It can be expressed as:

[0083]

[0084] Where ΔX is the track spacing.

[0085] Using formula (12), the Rayleigh wave phase velocity-depth curve or frequency-depth curve of the rock and soil between two adjacent detectors can be obtained.

[0086] Step 7: Draw the phase velocity imaging diagram.

[0087] Combining the phase velocity-depth curve or frequency-depth curve calculated by the two-way combination of multi-channel geophones, the Rayleigh wave phase velocity and frequency profile below the earthquake array can be obtained (such as Figure 9 ). Using Figure 9 , the spatial position of the civil air defense project and the position of the bedrock interface can be clearly determined.

Claims

1. A phase velocity imaging method based on similarity analysis of Rayleigh wave time-frequency characteristics, comprising the following steps: Step 1: Rayleigh wave exploration data collection: Basic data acquisition method: A series of geophones are deployed simultaneously along pre-selected survey lines and points. An 18-pound sledgehammer is used to strike the ground, generating seismic waves. The geophones receive the waves, and the main unit records the waveforms. The operating frequency band of Rayleigh wave exploration plays an important role in physical filtering during field data collection and digital filtering during data processing. The wavelength of a Rayleigh wave at a certain frequency is proportional to the detection depth. The relationship between the Rayleigh wave phase velocity, frequency, and detection depth can be expressed as: β is the depth coefficient, which is a constant related to the Poisson's ratio of the rock mass (as shown in Table 1); Table 1 Relationship between Poisson's ratio and depth coefficient The selection of the working frequency band for Rayleigh wave data acquisition mainly depends on the detection depth and shallow resolution required by the project. The maximum detection depth h max Determines the lower limit of the working frequency band (lowest frequency) and the minimum detection depth (or detection resolution) h min Determine the upper limit (maximum frequency) of the operating frequency band, then the maximum operating frequency is The minimum operating frequency is Step 2: Denoising of seismic waveform: The Rayleigh wave exploration work involved in the present invention can be carried out in cities or in the wild. If the exploration work is carried out in cities, it is severely disturbed by vibrations caused by various urban construction projects and vehicle types. Therefore, the original signal can be processed using a wavelet threshold filtering method (if the exploration work is carried out in the wild, denoising is not required). The wavelet threshold denoising process is divided into three parts: wavelet decomposition, threshold processing, and wavelet reconstruction. When wavelet decomposition is carried out, the selection of the number of layers is crucial: the larger the number of layers, the more obvious the different characteristics of the noise and signal performance, which is more conducive to denoising. However, at the same time, the reconstructed signal distortion will also be greater, which in turn affects the denoising effect. The frequency band range of wavelet decomposition is related to the sampling frequency. If the decomposition is carried out in N layers, the size of each frequency band is Fs / 2 / 2^N. Step 3: Rayleigh wave separation: After filtering, the Rayleigh wave needs to be separated from the seismic wave waveform. This step can be done directly on the time-domain seismic wave waveform. After the truncation is completed, the first point of the first Rayleigh wave waveform is used as the initial time, and the last point of the last waveform is used as the end time. The missing parts of each waveform are padded with 0 to ensure that the length of each waveform is the same, which is convenient for subsequent calculations. Step 4: Rayleigh wave time-frequency analysis: The time-frequency analysis expression can be written as: D(t,f,A)=TFA[g(t)] (2) Where TFA is a time-frequency analysis method, D(t,f,A) is the time-frequency analysis result, and g(t) is the Rayleigh wave signal; The Rayleigh wave is processed using the time-frequency analysis method to obtain the relationship between the time, frequency and amplitude of the Rayleigh wave; Step 5: Calculate the Rayleigh wave time difference: Set the amplitude threshold A0 and only analyze the frequency components with amplitudes greater than A0. Assume that the peak value of the time-frequency distribution amplitude of a waveform is A max , then A0=0.02A max For the nth and n+1th Rayleigh wave time-frequency distributions, extract the component with frequency f from them and obtain the time domain waveform of this frequency component, which are recorded as time series y1(x,f) and y2(x,f) respectively; Assuming that the number of signal sampling points is N0, the component y1(x,f) with frequency f is translated m(f) times in the direction of increasing time with a translation step of t0, and the root mean square error is used to evaluate the similarity; When the MSE takes the minimum value, it means that the similarity between y1(x,f) and y2(x,f) is the highest. If the corresponding number of translations is m0(f), it can be considered that the component with frequency f in the Rayleigh wave signal has an inter-channel time difference t between the nth channel and the n+1th channel. d =m0(f)t0; Step 6: Calculate the dispersion curve: After obtaining the inter-channel time difference of each frequency component, the Rayleigh wave phase velocity V of each frequency component is obtained. R It can be expressed as: Where ΔX is the track spacing; Using formula (1), the Rayleigh wave phase velocity-depth curve of the rock and soil between two adjacent detectors can be obtained; Step 7: Draw the phase velocity imaging diagram: Combining the phase velocity-depth curves calculated by pairwise combination of multi-channel geophones, the Rayleigh wave phase velocity-depth profile beneath the earthquake array can be obtained.

Images