A method for locating the direction of incoming waves based on a single-station seismograph

Abstract

This invention discloses a method for locating the direction of incoming waves based on a single-station seismograph. The steps are as follows: 1) Install the seismograph station in a tightly coupled manner and acquire signals; 2) Perform transform domain analysis on the acquired signals; 3) Calculate the correlation coefficients between each component to obtain the adaptive covariance matrix of the signal in the transform domain; 4) Calculate the polarizability based on the adaptive covariance matrix; 5) Filter data based on the distribution of polarizability in the transform domain to obtain the motion trajectory of the particles at the station; 6) Locate the source location by combining the motion trajectory of the particles and the wave type. This invention can achieve source direction estimation even with a single station in practical applications. For multiple stations forming a monitoring array, it can also effectively improve the positioning accuracy, achieving a good trade-off between positioning accuracy and hardware complexity.

Description

Technical Field

[0001] This invention belongs to the field of earthquake source location technology, specifically relating to a method for locating the direction of incoming waves based on a single-station seismograph. Background Technology

[0002] Seismic source location technology is a key technology widely used in the geological field. It can calculate the geographical location of an earthquake source based on signals emitted by the source and received by established seismic stations. Seismic source location technology not only helps pinpoint the exact location of an earthquake but also assists in locating underground resources such as water, minerals, oil, and natural gas. It also plays a crucial role in locating moving targets exhibiting vibrational characteristics. Currently, various methods can be used for seismic source location, including back azimuth estimation, vibration amplitude analysis, and waveform inversion.

[0003] Currently, the most common method for seismic source location is using dual-difference seismic location technology for back azimuth estimation. Specifically, this involves deploying a large array of seismograph stations near the epicenter and using the time difference between signals received by each station and known distances to infer the location of the signal source. This method of locating the source using the time difference of mechanical wave arrival is called dual-difference seismic location technology. This technology requires at least two stations, and the more stations there are, the higher the location accuracy. Motion amplitude analysis compares the signal amplitude detected by different stations and calculates the direction of the signal source based on a known signal propagation attenuation model. Waveform inversion, based on the signals acquired by the station array, further utilizes the amplitude, frequency, and other information carried by the signal to comprehensively invert and estimate the location of the signal source.

[0004] Generally, existing seismic source location techniques require at least two seismograph stations to determine the approximate direction of incoming waves on the horizontal plane. For more precise source location estimation, an array of even more stations is needed. Some technologies claim to require only one station to determine the direction of incoming waves on the horizontal plane, but these utilize six-component seismographs with strapdown connections of translational accelerometers and rotational angular velocity meters. Compared to ordinary seismic stations, these instruments have higher hardware complexity, and their location accuracy is limited by the weakest link effect, making further improvements difficult. Summary of the Invention

[0005] In view of the fact that traditional earthquake source location methods do not make full use of various signal characteristics and require at least two seismograph stations for location in both theory and practical application, this invention provides a method for directional location of signal source on the horizontal plane by making full use of the polarization characteristics of incoming wave signals.

[0006] Current waveform inversion methods primarily utilize signal amplitude and frequency information, with limited utilization of other signal characteristics, resulting in a waste of prior information. This invention, by comprehensively analyzing various aspects of the signal and fully leveraging its characteristics, performs a more in-depth analysis of the signal data received by the seismograph, enabling the determination of the incoming wave direction on the horizontal plane using only a single seismograph station. The source location method of this invention uses only a single seismograph station, and compared to the cross-correlation method requiring a six-component seismograph, this invention requires at least a single three-component seismograph to locate the signal source direction on the horizontal plane. Furthermore, this invention can utilize an array of multiple stations for higher-precision direction estimation, achieving a good trade-off between positioning accuracy and hardware complexity.

[0007] The technical solution adopted by this invention to solve its technical problem is:

[0008] Traditional methods for earthquake source location, such as estimating the azimuth at the time of the earthquake, are as follows: Figure 1 As shown, this invention improves upon the method by conducting in-depth analysis of the received signal and using polarization characteristics as an auxiliary basis for orientation, thereby ensuring that a single station can locate the signal source. The technical solution is primarily based on seismic waves generated by rock fracturing for location. Since earthquakes, plate tectonics, volcanic activity, or man-made underground blasting and mining all generate seismic waves, and different seismic waves have consistent polarization characteristics, this technical solution can also be applied to fields such as mining.

[0009] The flowchart of the wave direction localization technology based on a single-station three-component seismograph is as follows: Figure 2 As shown, a single station is tightly coupled to a solid propagation medium of the seismic source, such as rock strata in the region where the seismic source is located, to collect and process the signals emitted by the seismic source. The specific steps include:

[0010] S1: During installation, ensure that the z-axis points perpendicular to the solid propagation medium. A single seismograph station tightly coupled to the medium acquires seismic wave signals emitted by the source at a fixed sampling rate. Select a signal segment from the acquired seismic wave signals as the signal to be analyzed and processed. The selected signal segment (i.e., the signal to be analyzed) needs to contain a complete signal generated by the source. The z-axis of the seismograph station is a vertical axis perpendicular to the horizontal plane. The z-axis is a straight line, but it has a positive direction (i.e., direction). Generally, the direction perpendicular to the fixed surface and pointing upwards is taken as the positive direction.

[0011] S2: Perform time-frequency domain analysis on the signal to be analyzed, obtain the time-frequency domain parameter matrix, and save the time window width information used in the analysis process;

[0012] S3: For each element in the time-frequency domain parameter matrix, calculate the correlation coefficient between the signals on the two components (i.e., the corresponding components of the x-axis and y-axis) in the horizontal plane within a given time window width, with the time corresponding to the element as the center time, at the frequency point corresponding to the element, and obtain the time-frequency domain adaptive covariance matrix of the corresponding element; the given time window width can be selected according to actual needs, and this invention selects the time window width saved in the previous step;

[0013] S4: After calculating the time-frequency domain adaptive covariance matrix for each frequency point and each time step, calculate the eigenvalues ​​for each time-frequency domain adaptive covariance matrix and save them in descending order; two eigenvalues ​​will be calculated for the time-frequency domain adaptive covariance matrix corresponding to each element. The two eigenvalues ​​corresponding to the same time-frequency domain adaptive covariance matrix are saved in descending order. The eigenvalues ​​corresponding to different time-frequency domain adaptive covariance matrices are not compared.

[0014] S5: Divide the smaller eigenvalue by the larger eigenvalue and take the square root of the quotient. The result is the polarizability value at that time-frequency point. Arrange the result into a new matrix according to the corresponding time-frequency point index coordinates, which is the polarizability matrix.

[0015] S6: Based on the magnitude of the polarization at each time frequency point, the time frequency points corresponding to each polarization value are divided into two polarization types: linear polarization and elliptic polarization.

[0016] S7: Filter the previously acquired signal according to the filtering method described in S6, extract the corresponding linearly polarized signals from the two components on the horizontal plane (i.e., the corresponding components of the x-axis and y-axis), and calculate the back azimuth angle; the direction of the seismic source can be determined based on the obtained back azimuth angle;

[0017] S8: Since the determined back azimuth angle still has a 180° angular ambiguity, the elliptic polarization signals in the two components on the horizontal plane (i.e., the corresponding components on the x and y axes) and the vertical component (i.e., the corresponding component on the z axis) are extracted to determine the actual direction of arrival of the wave. That is, two sets of components are constructed: the first set consists of the x and z components, and the second set consists of the y and z components. For both sets of components, other types of wave signals are extracted following the operations in S3 to S7, and the characteristics of these signals are used to eliminate angular ambiguity. For example, Rayleigh wave signals are extracted, and then the uncertainty is eliminated based on the rotation direction of the Rayleigh wave.

[0018] Preferably, when performing time-frequency domain analysis on the acquired signal in step S2, wavelet transform can also be used for analysis, and subsequent steps can be performed accordingly in the time-scale domain.

[0019] Preferably, the method for calculating the azimuth angle in step S7 is to extract the Love wave signal collected by the station, determine the trajectory direction of the particle motion at the station caused by the Love wave on the horizontal plane, and then the direction of the incoming wave is the direction perpendicular to the trajectory of the particle motion on the horizontal plane.

[0020] Preferably, the elimination of 180° angular blur in step S8 mainly includes the following steps:

[0021] S81: If the approximate direction of the seismic source is known during installation (such as in man-made underground blasting, mining, and other application scenarios), the direction can be determined directly based on the installation location, eliminating angular ambiguity.

[0022] S82: If the approximate direction of the earthquake source cannot be clearly determined (e.g., in applications such as earthquake monitoring), then select the signal on the vertical component (i.e., the z-axis) and the two signals on the horizontal components (i.e., the x-axis and y-axis), and perform the operations described in S3 to S7 in sequence to extract other types of wave signals to eliminate angular ambiguity, such as extracting Rayleigh wave signals, and then eliminating angular ambiguity based on the rotation direction of Rayleigh waves. Specifically, this invention divides the three components of the signal to be analyzed into two groups. The first group includes the vertical component in the z-axis direction and the horizontal component in the x-axis direction, and the second group includes the vertical component in the z-axis direction and the horizontal component in the y-axis direction. Then, for each element in the time-frequency domain parameter matrix, the correlation coefficient between the two component signals in the first group is calculated at the corresponding frequency point, with the time corresponding to the element as the center time, within a given time window width, to obtain the time-frequency domain adaptive covariance matrix corresponding to the element. Then, the method of steps S4 to S7 is used to extract the first set type signal. And the correlation coefficient between the two component signals in the second group is calculated at the corresponding frequency point, with the time corresponding to the element as the center time, within a given time window width, to obtain the time-frequency domain adaptive covariance matrix corresponding to the element. Then, the method of steps S4 to S7 is used to extract the second set type signal. Then, the angular ambiguity of the source direction is eliminated according to the first set type signal and the second set type signal. The set type signal is a signal other than the linearly polarized signal.

[0023] Furthermore, the set type signal is an elliptical polarization signal or a Rayleigh wave signal.

[0024] The features of this invention are:

[0025] By utilizing the different polarization characteristics corresponding to different types of seismic wave signals, different types of seismic wave signals are separated. The motion trajectory of particles at the seismograph station under the influence of different seismic wave signals is obtained through numerical calculation. Finally, the location of the earthquake source is determined by combining the motion trajectory of the particles with the characteristics of the corresponding seismic waves themselves to invert the incoming wave direction.

[0026] Compared with the prior art, the positive effects of the present invention are as follows:

[0027] The wave direction localization technology based on a single-station seismograph proposed in this invention can estimate the direction of the seismic source with only a single station in practical applications. For the case of multiple stations forming a monitoring array, it can also effectively improve the positioning accuracy, achieving a good trade-off between positioning accuracy and hardware complexity. Attached Figure Description

[0028] Figure 1 This is an example diagram illustrating the method of estimating the azimuth angle based on the arrival time in earthquake source location.

[0029] Figure 2 This is a flowchart of the wave direction localization technology based on a single-station three-component seismograph proposed in this invention.

[0030] Figure 3 This is a diagram illustrating the specific implementation effect of the present invention applied to the motion trajectory of a mass point at a seismic source location station in a coal mine in China. Detailed Implementation

[0031] The flowchart of the wave direction localization technology based on a single-station three-component seismograph is as follows: Figure 2 As shown, a single station is tightly coupled to a solid propagation medium to acquire and process signals emitted by the seismic source. The specific steps include:

[0032] S1: During installation, ensure the z-axis points in the opposite direction to gravitational acceleration. A single seismograph station tightly coupled to the medium acquires signals emitted by the seismic source at a fixed sampling rate, ensuring that the complete propagation process of the signal to be measured from a given source is included within the start and end time period of the acquisition. After sampling is completed, the data is saved and denoted as s. x (t), s y (t),s z (t); where s x (t), s y (t) represents the horizontal component data of the seismic wave, s z (t) represents the vertical component data of the seismic wave;

[0033] S2: Perform time-frequency domain analysis on the acquired signal to obtain the time-frequency domain parameter matrix S, and save the time window width information T used in the analysis process;

[0034] S3: For each element S(t,f) in the time-frequency domain parameter matrix S, calculate the correlation coefficient I between the signal m component and the signal n component located in the horizontal plane (i.e., the corresponding components of the x-axis and y-axis) at the frequency point f corresponding to that element, with the time t corresponding to that element as the center time, within a given time window width T. mn(t,f), to obtain the time-frequency domain (or time-scale domain) adaptive covariance matrix MST, where m, n∈{x,y}; I mn This refers to including I xx I xy I yx and I yy Four types.

[0035] S4: After calculating the adaptive covariance matrix in the time-frequency domain for each frequency point and each time step, calculate the eigenvalues ​​for each matrix MST(t,f) and denote them as λ1, λ2 in descending order;

[0036] S5: Divide the smaller eigenvalue by the larger eigenvalue and take the square root of the quotient. The result is the polarizability value at that time-frequency point, denoted as polarizability ρ(t,f). Arrange ρ(t,f) into a new matrix according to the corresponding time-frequency point index coordinates, which is the polarizability matrix P.

[0037] S6: Based on the magnitude of the polarizability ρ(t,f) at each time frequency point, divide the time frequency point ST(t,f) corresponding to each polarizability value ρ(t,f) into two polarization types: linear polarization and elliptic polarization.

[0038] S7: Filter the previously acquired signal according to the filtering method described in S6, extract the corresponding linearly polarized signals from the two components on the horizontal plane (i.e., the corresponding components of the x-axis and y-axis), and calculate the azimuth angle α.

[0039] S8: Since there is still 180° angular ambiguity, the elliptic polarization signals in the two components on the horizontal plane (i.e., the components corresponding to the x-axis and y-axis) and the vertical component (i.e., the component corresponding to the z-axis) are extracted to determine the actual direction of the incoming wave.

[0040] Preferably, the transform domain analysis method used in step S2 includes, but is not limited to, short-time Fourier transform, S-transform, wavelet transform, etc. Taking the S-transform as an example, the specific formula is as follows:

[0041]

[0042] If wavelet transform is used to analyze the acquired signal in step S2, then subsequent steps are performed accordingly in the time-scale domain.

[0043] Preferably, the formula for calculating the elements of the time-frequency domain adaptive covariance matrix used in step S3 is as follows:

[0044]

[0045] Among them, I mnThe method for calculating (t,f) is to take the time t corresponding to the element as the center time at the corresponding frequency point f, and calculate the correlation coefficient between the m component and the n component within a given time window width.

[0046] Preferably, the method for calculating the polarizability ρ(t,f) in step S5 is to first calculate the eigenvalues ​​λ1 and λ2 of the matrix MST, where λ1>λ2, and then calculate the polarizability ρ(t,f). The specific formula is as follows:

[0047]

[0048] Preferably, the polarization type classification method in step S6 involves pre-setting a threshold. In this specific embodiment, considering the presence of noise, the threshold is set to 0.2. Filtering is performed based on the relationship between the polarizability ρ(t,f) and the threshold, and the time-frequency domain parameter matrix ST corresponding to the linearly polarized and elliptically polarized signals is calculated. re (t,f) and ST ep The specific formula for (t,f) is:

[0049]

[0050]

[0051] The inverse S-transform of the time-frequency domain parameter matrix is ​​then performed to obtain the corresponding linearly polarized and ellipticly polarized signals.

[0052] Preferably, the method for calculating the azimuth angle in step S7 is to extract the Love wave signal collected by the station, determine the trajectory direction of the particle motion at the station caused by the Love wave on the horizontal plane, and then the incoming wave direction is the direction perpendicular to the particle motion trajectory on the horizontal plane. The angle between the incoming wave direction and the x-component direction is denoted as α.

[0053] Preferably, in step S7, the trajectory direction of the particle's motion can be calculated by performing least-squares fitting on the discrete points of the particle's force on the two orthogonal components of the horizontal plane over a period of time, and using the arctangent trigonometric function to obtain the angle between the force direction and the horizontal orthogonal component. Due to the properties of the arctangent trigonometric function, the angle between the calculated incoming wave direction and the x-component direction is α, but the actual angle between the incoming wave direction and the x-component direction may be α or α+π.

[0054] Preferably, the elimination of 180° angular blur in step S8 mainly includes the following steps:

[0055] S81: If the approximate direction of the seismic source is known during installation (such as in man-made underground blasting, mining, and other application scenarios), the direction can be determined directly based on the installation location, eliminating angular ambiguity.

[0056] S82: If the approximate direction of the earthquake source cannot be clearly determined (e.g., in earthquake monitoring applications), select the signal on the vertical component (i.e., the z-axis) and the two signals on the horizontal components (i.e., the x-axis and y-axis), and perform the operations described in S3 to S7 in sequence to extract other types of wave signals to eliminate angular ambiguity. For example, extract the Rayleigh wave signal, and then determine the azimuth angle calculated in S7 as α or α+π based on whether the rotation direction of the Rayleigh wave is clockwise or counterclockwise, thereby eliminating the angular ambiguity caused by the properties of the arctangent trigonometric function.

[0057] Figure 3 A schematic diagram of the particle trajectory at the station mentioned in step S7 on the horizontal plane was drawn, and the direction perpendicular to the direction of motion is the direction of the earthquake source.

[0058] The present invention has been described in detail above, but it is obvious that the specific implementation of the present invention is not limited thereto. For those skilled in the art, various obvious modifications made to the method without departing from the spirit and scope of the claims are within the protection scope of the present invention.

Claims

1. A method for locating the direction of incoming waves based on a single-station seismograph, comprising the following steps: S1) A single seismograph station is installed on the source propagation medium, with its z-axis pointing opposite to the direction of gravitational acceleration in the vertical horizontal plane; the seismograph station is used to acquire seismic wave signals emitted from the source; each seismic wave signal includes three components s x (t), s y (t), s z (t), s x (t) represents the horizontal component signal along the x-axis at time t, s y (t) represents the horizontal component signal along the y-axis at time t, s z (t) represents the vertical component signal in the z-axis direction at time t; S2) Select a complete signal generated by the source from the collected seismic wave signals as the signal to be analyzed, and then perform time-frequency domain analysis on the signal to be analyzed to obtain the time-frequency domain parameter matrix, and save the time window width information used in the time-frequency domain analysis process; S3) For each element in the time-frequency domain parameter matrix, calculate the correlation coefficient between the two horizontal component signals within a given time window, taking the time corresponding to the element as the center time, and obtain the time-frequency domain adaptive covariance matrix corresponding to the element at the frequency point corresponding to the element. S4) Calculate the two eigenvalues ​​of each of the time-frequency domain adaptive covariance matrices; S5) Divide the smaller eigenvalue by the larger eigenvalue of the two eigenvalues ​​of each time-frequency domain adaptive covariance matrix, and take the square root of the quotient as the polarizability value at the corresponding time-frequency point. Then, form a polarizability matrix according to the index coordinates of the corresponding time-frequency points of each polarizability value. S6) Determine whether the seismic wave signal at each time and frequency point is linearly polarized or elliptically polarized based on the magnitude of the polarizability at each time and frequency point. S7) Filter the corresponding seismic wave signals according to the signal type of each time and frequency point, obtain the linear polarization signals of the two horizontal components in each seismic wave signal, and calculate the back azimuth angle of the corresponding seismic wave signal based on the linear polarization signals of the two horizontal components in the seismic wave signal. S8) Determine the source direction based on the back azimuth angle and perform angle ambiguity elimination on it; wherein the angle elimination method is as follows: if the approximate direction of the source has been determined when installing the seismograph station in step S1), then the source direction is ambiguity eliminated according to the installation position; otherwise, the three components in the signal to be analyzed are divided into two groups, the first group includes the vertical component in the z-axis direction and the horizontal component in the x-axis direction, and the second group includes the vertical component in the z-axis direction and the horizontal component in the y-axis direction; Then proceed to step S9); S9) For each element in the time-frequency domain parameter matrix, calculate the correlation coefficient between the two component signals in the first group within a given time window, with the time corresponding to the element as the center time, at the frequency point corresponding to the element, to obtain the time-frequency domain adaptive covariance matrix corresponding to the element; then extract the first set type signal using the method of steps S4) to S7); and calculate the correlation coefficient between the two component signals in the second group within a given time window, with the time corresponding to the element as the center time, at the frequency point corresponding to the element, to obtain the time-frequency domain adaptive covariance matrix corresponding to the element; then extract the second set type signal using the method of steps S4) to S7); then perform angular ambiguity elimination on the source direction based on the first set type signal and the second set type signal; the set type signal is a signal of other types besides linearly polarized signals.

2. The method according to claim 1, characterized in that, In step S2), wavelet transform is performed on the acquired signal to obtain a time-scale domain matrix that replaces the time-frequency domain parameter matrix.

3. The method according to claim 1 or 2, characterized in that, The time-frequency domain analysis methods include short-time Fourier transform, S-transform, and wavelet transform.

4. The method according to claim 1 or 2, characterized in that, The specified signal type is either an elliptic polarization signal or a Rayleigh wave signal.

5. The method according to claim 1 or 2, characterized in that, The polarizability values ​​at each time and frequency point are compared with a set threshold. If the value is greater than the set threshold, the seismic wave signal at the corresponding time and frequency point is determined to be linearly polarized; otherwise, it is determined to be elliptically polarized.

6. The method according to claim 5, characterized in that, The set threshold is 0.2.

Images