A method for extracting dispersion of high-order surface waves of background noise

Abstract

The application provides a background noise high-order surface wave dispersion extraction method, and belongs to the technical field of seismology. The method first calculates the cross-correlation function of continuous noise data and performs preprocessing to obtain an original dispersion energy spectrum; then a reference phase shift curve between the energy of the base order and the first high-order mode is extracted; the cross-correlation function is phase-shifted based on the curve, so that the base order and the high-order surface wave are separated in time, the high-order part at the negative time is retained and inverse phase shift is performed, and a cross-correlation function containing only the high-order surface wave is obtained; then the high-order surface wave dispersion energy spectrum is obtained through frequency-Bessel transformation, and is spliced with the base order part in the original spectrum to form an enhanced dispersion energy spectrum; finally, the base order and the high-order dispersion curves are extracted after quality control. The application can effectively suppress the interference of the base order sidelobe on the high-order mode, and significantly improve the continuity and signal-to-noise ratio of the high-order surface wave dispersion energy.

Description

Technical Field

[0001] This invention relates to the field of seismology, and more particularly to a method for extracting the dispersion of higher-order surface waves in background noise. Background Technology

[0002] The widespread application of background noise cross-correlation techniques has freed seismic imaging from dependence on natural earthquakes. A relatively complete surface wave dispersion imaging technique has been developed around background noise, becoming an important tool for studying the near-surface and lithospheric shear wave velocity structure. Traditionally, background noise surface wave imaging techniques have primarily focused on the fundamental mode surface waves. In recent years, higher-order surface wave modes and their dispersion characteristics have also received widespread attention because higher-order modes can provide higher vertical resolution constraint information for model inversion. Therefore, combining the use of fundamental and higher-order dispersion in imaging allows for more refined modeling of stratigraphic structures. Theoretically, higher-order surface wave modes, like the fundamental mode, are inherent characteristics of layered half-space models. Based on the dispersion energy spectrum calculated using methods such as the frequency-Bessel transform method or the phase shift method, the information of higher-order surface wave modes contained in background noise can be widely identified. However, the dispersion energy of current higher-order modes is often accompanied by narrow bandwidth and low signal-to-noise ratio (SNR). The core reason for this phenomenon is the weak amplitude of higher-order modes, which leads to them being frequently interfered with or even masked by various artifacts in the dispersion energy spectrum. Among these, the sidelobes of the fundamental surface wave are one of the important sources of artifacts. These artifacts cut off and distort the originally smooth and continuous higher-order energy in the dispersion energy spectrum, ultimately resulting in a low SNR of the picked-up higher-order mode dispersion curve.

[0003] To suppress artifacts, some scholars have proposed a multi-time-window frequency-Bessel transform method. [2] The core of this method lies in using a pre-designed window function to suppress the fundamental mode with a large amplitude, thereby reducing the amplitude of the fundamental sidelobes in the dispersive energy spectrum and mitigating the interference of the sidelobes on the energy of higher-order modes. The limitation of this method is that the waveforms of the fundamental mode and higher-order modes partially overlap in time; therefore, this method cannot effectively protect higher-order modes while suppressing the fundamental waveform. Summary of the Invention

[0004] To overcome the shortcomings of existing technologies, this invention provides a method for extracting the dispersion of higher-order surface waves from background noise. This invention proposes a separation technique between the fundamental and higher-order modes, using it as a pre-process for calculating the dispersion energy spectrum. Based on this mode separation technique, this method can effectively suppress the interference of the fundamental surface wave sidelobes on the dispersion energy of higher-order surface waves, thereby improving the signal-to-noise ratio of the higher-order surface wave portion in the dispersion energy spectrum.

[0005] This invention is achieved through the following technical solution: a method for extracting the dispersion of high-order surface waves in background noise, characterized by comprising the following steps: Step 1: Calculate the cross-correlation function of continuous noise data; Step S1-1: Divide the noise data continuously recorded by each station into multiple time periods. The length of the time period depends on the research scale. For arrays with a coverage of hundreds of kilometers or more, the data is usually divided into days, while for small arrays, it can be shortened to several hours. Step S1-2: Perform basic processing on each data segment, including mean removal, trend removal, filtering, and downsampling; Steps S1-3: After the above basic processing, given the sliding window length and the amount of time window movement each time, perform cross-correlation calculations on the noise data of all station pairs within each sliding window; utilize the properties of Fourier transform to convert the time-domain cross-correlation operation into a frequency-domain multiplication operation: (2) in, The cross-correlation function representing the frequency domain. Represents frequency domain noise signal, and Represents the coordinates of two stations. ω represents the distance between the two stations, and ω represents the frequency. Steps S1-4: Superimpose all the calculation results by station pair to obtain the cross-correlation function; symmetrically superimpose the cross-correlation function of the negative time axis onto the positive time axis, and simultaneously set the negative time axis signal to 0 to obtain the cross-correlation function for subsequent processing; Step S2: Perform multichannel dispersion analysis on the cross-correlation function obtained in step S1 to obtain the original dispersion energy spectrum; use the adjacent channel superposition and group velocity filtering method to suppress interference signals on the cross-correlation function; Step S3: Extraction of the reference phase shift curve; From the original dispersion energy spectrum, determine the energy range of the surface wave fundamental and higher-order modes. The fundamental signal appears as a single energy band with the strongest energy and lowest velocity, while the higher-order signal appears as multiple energy bands with higher velocity and weaker energy. Pick a curve with a gradually decreasing phase velocity from low to high frequency between the fundamental and higher-order energies as the reference phase shift curve for subsequent mode separation, and record its corresponding wavenumber as... The curve lies between the surface wave fundamental order and the first higher-order mode energy. Step S4: Perform mode separation; first, transform each cross-correlation function to the frequency domain, then apply a phase shift, the phase shift factor being given by the reference phase shift curve picked up in step S3; for the frequency domain cross-correlation function of the nth station pair... , Given the corresponding station pair spacing, the phase shift factor can be determined as follows: The new cross-correlation function after applying the phase shift can be calculated by equation (4): (4) superscript ~ This indicates that the signal has undergone a phase shift; Then the cross-correlation function of the phase shift Applying an inverse Fourier transform yields the time series of phase-shift cross-correlation functions. For each phase-shift cross-correlation function, only the higher-order surface wave portion of the signal at the negative time is retained, and it is transformed back to the frequency domain, denoted as . Subsequently, by applying an inverse phase shift to the retained negative-time signal using equation (5), a cross-correlation function containing only higher-order surface waves can be obtained: (5); Step S5: Calculation of enhanced dispersion energy spectrum; cross-correlation function for higher-order surface wave components. The frequency-Bessel transformation is performed again using equation (3) to obtain the higher-order surface wave dispersion energy spectrum. The part of the original dispersion energy spectrum below the reference phase shift curve and the part of the higher-order surface wave dispersion energy spectrum above the reference phase shift curve are spliced ​​together to obtain a new dispersion energy spectrum, which is called the enhanced dispersion energy spectrum. Step S6, Quality Control: Determine whether the reference phase shift curve is centered based on the positions of the fundamental and higher-order energies in the enhanced dispersion energy spectrum. If it is not centered, return to step S3 and re-execute based on the enhanced dispersion energy spectrum. Step S7: Extract the fundamental and higher order dispersion curves from the enhanced dispersion energy spectrum.

[0006] As a preferred approach, in steps S1-2, mean removal involves subtracting the arithmetic mean from the entire seismic record's time series data, making the signal mean zero; filtering and downsampling involve using bandpass or lowpass filters to filter the signal into the desired frequency band, followed by reducing the signal sampling rate; detrending removes the linear trend by performing a linear fit on the noisy time series to obtain an optimal fit line. Then subtract this fitted line from the original signal: (1) in and These represent the original signal and the signal after de-linearization, respectively. As a preferred embodiment, the adjacent channel superposition method in step S2 divides the cross-correlation function into several groups according to its station spacing, averages the cross-correlation functions in each group to obtain a new cross-correlation function, and uses the average station spacing within the group as the station spacing of the new cross-correlation function; the group velocity filtering method first estimates the fastest arrival time T of the surface wave based on the surface wave signal shown in the time-domain cross-correlation function, then suppresses signals earlier than time T, and finally uses the frequency-Bessel transform method to obtain the original dispersion energy spectrum of the surface wave; the frequency-Bessel transform method is implemented by the following integration: (3) in It is a dispersive energy spectrum. Represents the 0th order Bessel function of the first kind. Represents wave number.

[0007] By employing the above technical solutions, this invention offers the following advantages compared to existing technologies: It utilizes mode separation technology to separate the fundamental and higher-order modes in the background noise cross-correlation function, suppressing the interference of fundamental surface wave sidelobes on the dispersion energy of higher-order surface waves, thereby significantly improving the signal-to-noise ratio of the higher-order surface wave portion of the dispersion energy spectrum. Compared to existing technologies, this invention effectively protects higher-order waveforms while separating the fundamental mode, resulting in more complete higher-order mode dispersion energy spectrum information and a higher signal-to-noise ratio.

[0008] Additional aspects and advantages of the invention will become apparent in the following description or may be learned by practice of the invention. Attached Figure Description

[0009] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which: Figure 1 This is a flowchart illustrating the implementation of the present invention. Figure 2 The image shows a partial distribution map of USArray stations (left) and the noise cross-correlation function (right). Figure 3 This is the original dispersion energy spectrum; Figure 4 To enhance the dispersion energy spectrum. Detailed Implementation

[0010] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0011] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0012] The following is combined with Figures 1 to 4 The dispersion extraction method for higher-order surface waves of background noise according to embodiments of the present invention will be described in detail.

[0013] This invention proposes a method for extracting the dispersion of higher-order surface waves in background noise. The method flow is as follows: Figure 1 As shown, the specific steps include: Step 1: Calculate the cross-correlation function of continuous noise data; Step S1-1: Divide the noise data continuously recorded by each station into multiple time periods. The length of the time period depends on the research scale. For arrays with a coverage of hundreds of kilometers or more, the data is usually divided into days, while for small arrays, it can be shortened to several hours. Steps S1-2 involve performing basic processing on each data segment, including mean removal, trend removal, filtering, and downsampling. Mean removal involves subtracting the arithmetic mean from the entire seismic record's time series data, making the signal mean zero. Since the original data bandwidth recorded by the instrument is usually much larger than the required signal bandwidth, filtering and downsampling use bandpass or lowpass filters to filter the signal into the required bandwidth, and then reduce the signal sampling rate to reduce subsequent computation and storage requirements. Trend removal aims to remove long-term trends in the signal. Trend removal uses linear trend removal, which involves linearly fitting the noisy time series to obtain an optimal fit line. Then subtract this fitted line from the original signal: (1) in and These represent the original signal and the signal after de-linearization, respectively. Steps S1-3: After the above basic processing, given the sliding window length and the amount of time window movement each time, perform cross-correlation calculations on the noise data of all station pairs within each sliding window; utilize the properties of Fourier transform to convert the time-domain cross-correlation operation into a frequency-domain multiplication operation: (2) in, The cross-correlation function representing the frequency domain. Represents frequency domain noise signal, and Represents the coordinates of two stations. ω represents the distance between the two stations, and ω represents the frequency. Step S1-4: Superimpose all the calculation results by station pair to obtain the cross-correlation function; since the cross-correlation function has signals on both the positive and negative time axes, the cross-correlation function on the negative time axis is superimposed on the positive time axis in a time-symmetric manner, and the negative time axis signal is set to 0 at the same time to obtain the cross-correlation function for subsequent processing; Step S2: Perform multichannel dispersion analysis on the cross-correlation function obtained in Step S1 to obtain the original dispersion energy spectrum; suppress interference signals using the neighboring channel superposition and group velocity filtering methods on the cross-correlation function; the neighboring channel superposition method divides the cross-correlation function into several groups (e.g., 200 groups) according to their station spacing, averages the cross-correlation functions in each group to obtain a new cross-correlation function, and uses the average station spacing within the group as the station spacing of the new cross-correlation function. This method can suppress random interference in the dispersion energy spectrum to a certain extent; the group velocity filtering method aims to remove interference waves whose propagation speed deviates significantly from the surface wave velocity range. This method first estimates the fastest arrival time T of the surface wave based on the surface wave signal shown in the time-domain cross-correlation function, then suppresses signals earlier than time T, and finally uses the frequency-Bessel transform method to obtain the original dispersion energy spectrum of the surface wave; the frequency-Bessel transform method is implemented by the following integration: (3) in It is a dispersive energy spectrum. Represents the 0th order Bessel function of the first kind. Represents wave number.

[0014] Step S3: Extraction of the reference phase shift curve; From the original dispersion energy spectrum, determine the energy range of the surface wave fundamental and higher-order modes. The fundamental signal appears as a single energy band with the strongest energy and lowest velocity, while the higher-order signal appears as multiple energy bands with higher velocity and weaker energy. Pick a curve with a gradually decreasing phase velocity from low to high frequency between the fundamental and higher-order energies as the reference phase shift curve for subsequent mode separation, and record its corresponding wavenumber as... There are no strict requirements for picking this curve; it is sufficient to place it as close as possible to the energy between the surface wave fundamental order and the first higher-order mode. If the higher-order mode energy is difficult to identify, it can be picked first along the upper bound of the fundamental surface wave energy.

[0015] Step S4: Perform mode separation; first, transform each cross-correlation function to the frequency domain, then apply a phase shift, the phase shift factor being given by the reference phase shift curve picked up in step S3; for example, for the frequency domain cross-correlation function of the nth station pair... , Given the corresponding station pair spacing, the phase shift factor can be determined as follows: The new cross-correlation function after applying the phase shift can be calculated by equation (4): (4) superscript ~ This indicates that the signal has undergone a phase shift; Then the cross-correlation function of the phase shift Applying an inverse Fourier transform yields the time series of the phase-shift cross-correlation function. This time series exhibits the following characteristics: the waveforms of the fundamental surface wave are primarily concentrated in the positive time intervals, while higher-order surface waves primarily appear in the negative time intervals; both the fundamental and higher-order surface waves retain dispersion characteristics. For each phase-shift cross-correlation function, only the higher-order surface wave portion of the signal at the negative time intervals is retained, and it is transformed back to the frequency domain, denoted as […]. Subsequently, by applying an inverse phase shift to the retained negative-time signal using equation (5), a cross-correlation function containing only higher-order surface waves can be obtained: (5); Step S5: Calculation of enhanced dispersion energy spectrum; cross-correlation function for higher-order surface wave components. The frequency-Bessel transformation is performed again using equation (3) to obtain the higher-order surface wave dispersion energy spectrum. The part of the original dispersion energy spectrum below the reference phase shift curve and the part of the higher-order surface wave dispersion energy spectrum above the reference phase shift curve are spliced ​​together to obtain a new dispersion energy spectrum, which is called the enhanced dispersion energy spectrum. Step S6, Quality Control: Determine whether the reference phase shift curve is centered based on the positions of the fundamental and higher-order energies in the enhanced dispersion energy spectrum. If it is not centered, return to step S3 and re-execute based on the enhanced dispersion energy spectrum. Step S7: Extract the fundamental and higher order dispersion curves from the enhanced dispersion energy spectrum.

[0016] To verify the advantages of this mode separation method in high-order surface wave extraction, experiments were conducted using background noise data recorded by some USArray stations in the central United States. Figure 2 The experiment used continuous background noise data recorded by 96 stations between June 1 and December 1, 2011. The raw data from each station was split into multiple data files stored on a daily basis. After basic processing such as mean removal, detrending, and downsampling, pairwise cross-correlation was performed between stations to obtain the cross-correlation function. Figure 2 Using the method described in step 2, dispersion analysis is performed on the above noise cross-correlation function to obtain the original dispersion energy spectrum. Figure 3 Then, using the methods described in steps 3 to 6, the enhanced dispersion energy spectrum is obtained. Figure 4Compared to the original dispersive energy spectrum, the enhanced dispersive energy spectrum shows more continuous and wider bandwidth higher-order surface wave dispersive energy. Furthermore, while only the third higher order can be identified in the original dispersive energy spectrum, the enhanced dispersive energy spectrum clearly identifies the fifth higher order. This indicates that the proposed method can significantly improve the pickup performance of higher-order surface waves in background noise.

[0017] In the description of this invention, the term "a plurality of" refers to two or more. Unless otherwise explicitly defined, the terms "upper," "lower," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, and are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. The terms "connection," "installation," "fixing," etc., should be interpreted broadly. For example, "connection" can be a fixed connection, a detachable connection, or an integral connection; it can be a direct connection or an indirect connection through an intermediate medium. For those skilled in the art, the specific meaning of the above terms in this invention can be understood according to the specific circumstances.

[0018] In the description of this specification, the terms "one embodiment," "some embodiments," "specific embodiment," etc., refer to a specific feature, structure, material, or characteristic described in connection with that embodiment or example, which is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0019] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for extracting the dispersion of higher-order surface waves in background noise, characterized in that, Specifically, the following steps are included: Step 1: Calculate the cross-correlation function of continuous noise data; Step S1-1: Divide the noise data continuously recorded by each station into multiple time periods. The length of the time period depends on the research scale. For arrays with a coverage of hundreds of kilometers or more, the data is usually divided into days, while for small arrays, it can be shortened to several hours. Step S1-2: Perform basic processing on each data segment, including mean removal, trend removal, filtering, and downsampling; Steps S1-3: After the above basic processing, given the sliding window length and the amount of time window movement each time, perform cross-correlation calculations on the noise data of all station pairs within each sliding window; utilize the properties of Fourier transform to convert the time-domain cross-correlation operation into a frequency-domain multiplication operation: (2) Where C represents the cross-correlation function in the frequency domain, U represents the frequency domain noise signal, x1 and x2 represent the coordinates of the two stations, r represents the distance between the two stations, and ω represents the frequency. Steps S1-4: Superimpose all the calculation results by station pair to obtain the cross-correlation function; symmetrically superimpose the cross-correlation function of the negative time axis onto the positive time axis, and simultaneously set the negative time axis signal to 0 to obtain the cross-correlation function for subsequent processing; Step S2: Perform multichannel dispersion analysis on the cross-correlation function obtained in step S1 to obtain the original dispersion energy spectrum; use the adjacent channel superposition and group velocity filtering method to suppress interference signals on the cross-correlation function; Step S3: Extraction of the reference phase shift curve; From the original dispersion energy spectrum, determine the energy range of the surface wave fundamental and higher-order modes. The fundamental signal appears as a single energy band with the strongest energy and lowest velocity, while the higher-order signal appears as multiple energy bands with higher velocity and weaker energy. Pick a curve with a gradually decreasing phase velocity from low to high frequency between the fundamental and higher-order energies as the reference phase shift curve for subsequent mode separation, and record its corresponding wavenumber as... The curve lies between the surface wave fundamental order and the first higher-order mode energy. Step S4: Perform mode separation; first, transform each cross-correlation function to the frequency domain, then apply a phase shift, the phase shift factor being given by the reference phase shift curve picked up in step S3; for the frequency domain cross-correlation function of the nth station pair... r n Given the corresponding station pair spacing, the phase shift factor can be determined as follows: The new cross-correlation function after applying the phase shift can be calculated by equation (4): (4) superscript ~ This indicates that the signal has undergone a phase shift; Then the cross-correlation function of the phase shift Applying an inverse Fourier transform yields the time series of phase-shift cross-correlation functions. For each phase-shift cross-correlation function, only the higher-order surface wave portion of the signal at the negative time is retained, and it is transformed back to the frequency domain, denoted as . Subsequently, by applying an inverse phase shift to the retained negative-time signal using equation (5), a cross-correlation function containing only higher-order surface waves can be obtained: (5)； Step S5: Calculation of enhanced dispersion energy spectrum; cross-correlation function for higher-order surface wave components. The frequency-Bessel transformation is performed again using equation (3) to obtain the higher-order surface wave dispersion energy spectrum. The part of the original dispersion energy spectrum below the reference phase shift curve and the part of the higher-order surface wave dispersion energy spectrum above the reference phase shift curve are spliced ​​together to obtain a new dispersion energy spectrum, which is called the enhanced dispersion energy spectrum. Step S6, Quality Control: Determine whether the reference phase shift curve is centered based on the positions of the fundamental and higher-order energies in the enhanced dispersion energy spectrum. If it is not centered, return to step S3 and re-execute based on the enhanced dispersion energy spectrum. Step S7: Extract the fundamental and higher order dispersion curves from the enhanced dispersion energy spectrum.

2. The method for extracting the dispersion of higher-order surface waves in background noise according to claim 1, characterized in that... In steps S1-2, the mean removal process involves subtracting the arithmetic mean from the time series data of the entire earthquake record, making the mean of the signal zero; filtering and downsampling involve using a bandpass or lowpass filter to filter the signal into the frequency band to be studied, and then reducing the signal sampling rate. Detrending removes linear trends by performing a linear fit on the noisy time series to obtain an optimal fit line. Then subtract this fitted line from the original signal: (1) in and These represent the original signal and the signal after de-linearization, respectively.

3. The method for extracting the dispersion of higher-order surface waves in background noise according to claim 1, characterized in that... In step S2, the adjacent channel superposition method is to divide the cross-correlation function into several groups according to its station spacing, take the average of the cross-correlation functions in each group to obtain a new cross-correlation function, and take the average of the station spacing in the group as the station spacing of the new cross-correlation function. The group velocity filtering method first estimates the fastest arrival time T of the surface wave based on the surface wave signal displayed in the time-domain cross-correlation function. Then, signals arriving earlier than T are suppressed. Finally, the original dispersion energy spectrum of the surface wave is obtained using the frequency-Bessel transform method. The frequency-Bessel transform method is implemented through the following integration: (3) in It is a dispersive energy spectrum. represents a 0th-order Bessel function of the first kind, and k represents the wave number.

Images