A method and system for earthquake source location based on surface wave dispersion back propagation and superposition

Abstract

The present invention proposes a source location method and system based on surface wave dispersion backpropagation and superposition, by obtaining the surface wave phase velocity dispersion curve; based on the surface wave phase velocity dispersion curve, performing surface wave dispersion backpropagation processing to extract the seismic waveform envelope curve; based on the seismic waveform envelope curve, performing grid search and multi-channel superposition processing to obtain the source epicenter position. The problem that the processing scheme based on the arrival time information of relatively high-frequency body wave signals in the prior art cannot be applied to the source location of surface waves is solved, and the source epicenter location is realized by using relatively low-frequency and higher signal-to-noise ratio surface wave signals, and the accuracy and stability of the source location method of surface wave dispersion backpropagation and superposition are guaranteed.

Description

Technical Field

[0001] The present application relates to the technical field of earthquake location, and particularly relates to a method and system for earthquake source location based on surface wave dispersion backpropagation and superposition. Background Art

[0002] Earthquake source location is one of the most basic tasks in seismological research, and is of great significance for many aspects such as natural disaster monitoring, prevention, prediction, post-disaster assessment and emergency rescue.

[0003] Natural earthquakes, as well as a large number of instantaneous slip events and surface processes that are effectively coupled to the ground, such as landslides, glacier calving, volcanic eruptions, hurricanes and explosions, etc., may act as earthquake sources to excite outward-propagating seismic waves, which are ultimately recorded by artificially deployed observation systems such as seismic stations or DAS (Distributed Acoustic Sensing). The process of using these observed seismic records to infer the location where the earthquake source occurred is called earthquake source location.

[0004] The seismic waves excited by the earthquake source include body waves (P waves and S waves) and surface waves. Traditional earthquake source location techniques mainly process based on the arrival time information of relatively high-frequency body wave signals. However, in some cases, for example, some earthquake sources are difficult to excite high-signal-to-noise ratio body wave signals, the relatively high-frequency body wave signals rapidly attenuate during propagation, and in DAS observation data, the signal-to-noise ratio of body waves is usually low, making it difficult to extract high-quality body wave arrival time information in the prior art, resulting in limited accuracy of traditional earthquake source location techniques and making it difficult to carry out.

[0005] In addition to body wave phases, the earthquake source excitation also generates relatively low-frequency surface wave signals. In practice, the surface waves excited by the earthquake source directly acting on the Earth's surface usually have a higher signal-to-noise ratio. However, due to the dispersion characteristics of surface waves, that is, different frequency components of surface waves have different propagation speeds, the surface waves will have a wider wave packet, that is, the waveform diverges, as the propagation distance increases, making it difficult to accurately extract the arrival time information of surface waves. Therefore, it is difficult to directly apply traditional earthquake source location techniques to the application of surface waves.

[0006] Therefore, the prior art still needs to be improved and enhanced. Summary of the Invention

[0007] The technical problem to be solved by the present application is that traditional earthquake source location techniques, that is, the processing method based on the arrival time information of relatively high-frequency body wave signals, are difficult to be applied to the processing of surface waves; and the problem that is difficult to extract due to the dispersion effect of surface waves. In view of the deficiencies of the prior art, the present invention provides a method and system for earthquake source location based on surface wave dispersion backpropagation and superposition.

[0008] To solve the above technical problems, a first aspect of the embodiments of the present application provides a source location method based on surface wave dispersion backpropagation and superposition. The source location method based on surface wave dispersion backpropagation and superposition includes:

[0009] Obtain the surface wave phase velocity dispersion curve;

[0010] Based on the surface wave phase velocity dispersion curve, perform surface wave dispersion backpropagation processing to extract the seismic waveform envelope curve;

[0011] Based on the seismic waveform envelope curve, perform grid search and multi-channel superposition processing to obtain the epicenter location of the source.

[0012] Optionally, the step of obtaining the surface wave phase velocity dispersion curve includes:

[0013] Based on the common velocity model and the dominant frequency band of the observed data, perform forward calculation to obtain the surface wave phase velocity dispersion curve.

[0014] Optionally, the step of obtaining the surface wave phase velocity dispersion curve includes:

[0015] Obtain the seismic trace with the minimum offset and determine it as the virtual source;

[0016] Based on the observed data of the remaining seismic traces, perform cross-correlation calculation with respect to the virtual source to obtain the array data of the cross-correlation function;

[0017] Use the surface wave dispersion extraction algorithm to extract the surface wave phase velocity dispersion curve from the array data of the cross-correlation function.

[0018] Optionally, the step of performing surface wave dispersion backpropagation processing based on the surface wave phase velocity dispersion curve to extract the seismic waveform envelope curve includes:

[0019] Obtain N-channel array seismic data;

[0020] Based on the N-channel array seismic data, calculate the corresponding seismic spectrum.

[0021] Optionally, the step of performing surface wave dispersion backpropagation processing based on the surface wave phase velocity dispersion curve to extract the seismic waveform envelope curve includes: performing surface wave dispersion backpropagation processing based on the seismic spectrum to obtain the seismic spectrum after surface wave dispersion backpropagation.

[0022] Optionally, the step of performing surface wave dispersion backpropagation processing based on the surface wave phase velocity dispersion curve to extract the seismic waveform envelope curve includes:

[0023] Based on the time-domain seismic waveform data after surface wave dispersion backpropagation, calculate the envelope curve of the time-domain seismic waveform.

[0024] Optionally, the step of performing grid search and multi-channel stacking processing based on the seismic waveform envelope curve to obtain the source epicenter position includes:

[0025] Based on the envelope curve of the seismic waveform in the time domain, obtain the maximum amplitude of the envelope curve of the seismic waveform in the time domain;

[0026] Based on the array seismic data, stack the maximum amplitudes of the envelope curves of the seismic waveforms in the time domain, and calculate the stacked amplitude corresponding to each coordinate point.

[0027] Optionally, the step of performing grid search and multi-channel stacking processing based on the seismic waveform envelope curve to obtain the source epicenter position further includes:

[0028] Compare all the stacked amplitudes, determine the coordinate point corresponding to the maximum stacked amplitude, and determine the coordinate point corresponding to the maximum stacked amplitude as the source epicenter position.

[0029] In addition, to achieve the above object, the present invention also provides a source location system based on surface wave dispersion backpropagation and stacking, and the system includes:

[0030] A surface wave phase velocity dispersion curve extraction module for obtaining a surface wave phase velocity dispersion curve;

[0031] A seismic waveform envelope curve extraction module for performing surface wave dispersion backpropagation processing based on the surface wave phase velocity dispersion curve and extracting a seismic waveform envelope curve;

[0032] A source epicenter position acquisition module for performing grid search and multi-channel stacking processing based on the seismic waveform envelope curve to obtain the source epicenter position.

[0033] As can be seen from the above, in the solution of the present invention, a source location method and system based on surface wave dispersion backpropagation and stacking are proposed. By obtaining a surface wave phase velocity dispersion curve; performing surface wave dispersion backpropagation processing based on the surface wave phase velocity dispersion curve and extracting a seismic waveform envelope curve; performing grid search and multi-channel stacking processing based on the seismic waveform envelope curve to obtain the source epicenter position. It solves the problem that the processing scheme based on the arrival time information of body wave signals with relatively high frequencies in the prior art cannot be applied to the source location of surface waves, realizes the source epicenter location using surface wave signals with relatively low frequencies and higher signal-to-noise ratios, and ensures the accuracy and stability of the source location method of surface wave dispersion backpropagation and stacking. Description of the Drawings

[0034] To more clearly illustrate the technical solutions in the embodiments of the present application, the following briefly introduces the accompanying drawings required for the description of the embodiments. Obviously, the accompanying drawings in the following description are only some embodiments of the present application. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on the structures shown in these drawings.

[0035] Figure 1 It is a flowchart of the seismic source location method based on surface wave dispersion backpropagation and superposition provided by the present invention.

[0036] Figure 2 It is a schematic diagram of the observation system provided by the present invention.

[0037] Figure 3 It is a multi-channel seismic velocity field record observed by a traditional seismic station provided by the present invention.

[0038] Figure 4 It is a multi-channel seismic strain rate field record observed by DAS provided by the present invention.

[0039] Figure 5 It is an effect diagram of seismic source location applied to seismic station observation data provided by the present invention.

[0040] Figure 6 It is an effect diagram of seismic source location applied to DAS observation data provided by the present invention.

[0041] Figure 7 It is a schematic diagram of the seismic source location system based on surface wave dispersion backpropagation and superposition provided by the present invention. Detailed implementation manners

[0042] To make the purpose, technical solutions and advantages of the present invention clearer and more definite, the following further elaborates on the present invention with reference to the accompanying drawings and by way of examples. It should be understood that the specific examples described herein are only used to explain the present invention and are not used to limit the present invention.

[0043] The seismic source location method based on surface wave dispersion backpropagation and superposition described in the present invention can be based on traditional seismic observation data or new DAS observation data. As Figure 2 shown, it is a schematic diagram of the observation system.

[0044] Traditional seismic observation uses seismic stations to establish an observation system, which can collect information on seismic displacement fields, velocity fields, or acceleration fields. In this embodiment, a four-layer half-space geological model is taken as an example, and a single-force point source is used as the seismic source excitation.

[0045] DAS is a new type of dense array observation technology that has developed rapidly in recent years. It has advantages such as low observation cost and strong repeatability, and is very suitable for the location and monitoring of seismic sources. It can be used to observe information on the seismic strain field or strain rate field. In this embodiment, a four-layer half-space geological model is taken as an example, a single-force point source is used as the seismic source excitation, and two linear and perpendicular DAS arrays are used for observation.

[0046] Among them, Figure 2 The black five-pointed star represents the true seismic source position, the black straight line represents the optical fiber layout position, and they form the DAS observation system. The black triangle represents the seismic station layout position, and they form a randomly distributed array observation system.

[0047] The seismic source location method based on surface wave dispersion backpropagation and superposition according to the preferred embodiment of the present invention, as Figure 1 shown, the seismic source location method based on surface wave dispersion backpropagation and superposition includes the following steps:

[0048] Step S100, obtain the surface wave phase velocity dispersion curve.

[0049] Specifically, the surface wave phase velocity dispersion curve c(ω) is a curve that describes the variation of the surface wave phase velocity with frequency. The surface wave propagating in the medium has different phase velocities at different frequencies, that is, the velocity of the wavefront propagation. The phase velocity dispersion curve can be obtained through experiments or calculations. In this embodiment, the method for obtaining the surface wave phase velocity dispersion curve c(ω) can be carried out by the following method:

[0050] (1) Forward calculation method: According to the common velocity model of the area where the observation system is located and the dominant frequency band of the observation data, use the forward calculation method to calculate the surface wave phase velocity dispersion curve. The forward calculation of the surface wave phase velocity dispersion curve is essentially to solve the roots of the surface wave dispersion equation. The surface wave phase velocity is a function of geological model parameters such as longitudinal and transverse wave velocities and density, and is particularly sensitive to the transverse wave velocity. During the forward calculation process, different surface wave phase velocity dispersion curves can be obtained by changing different transverse wave velocity structures. Therefore, to obtain the surface wave phase velocity dispersion curve related to the observation data, it is necessary to obtain as accurate an underground velocity model as possible in the area where the observation system is located in advance.

[0051] (2) Cross-correlation method: Use the seismic trace with the smallest offset as the virtual source, and perform cross-correlation calculations on the observation data of other seismic traces with respect to the virtual source. After obtaining the array data of the cross-correlation function through cross-correlation calculations, a surface wave dispersion extraction algorithm, such as the frequency-Bessel transform method, can be used to extract the surface wave phase velocity dispersion curve.

[0052] In seismology, virtual sources are usually used to represent the excitation positions of virtual seismic wave sources, that is, imaginary starting points in the seismic wave propagation path, for reconstructing the propagation process of seismic waves relative to the virtual source points. This reconstruction process is achieved through cross-correlation operations. Cross-correlation is a widely adopted mathematical method, and the cross-correlation function can reveal the similar or correlated characteristics between two seismic signals. In signal processing, the cross-correlation function is usually used to analyze information such as the correlation, time shift, or phase difference between two signals.

[0053] Seismic trace observation data refers to the seismic wave signals recorded by seismic instruments at observation points after a seismic event occurs. The observation data usually includes information such as the amplitudes and arrival times of seismic waves at different time points. Each observation point is a seismic trace, and multiple observation points can collect seismic data of multiple seismic traces (i.e., multi-traces).

[0054] Multi-trace data is usually also often referred to as array data, which refers to a data set stored in the form of a two-dimensional array. The array data of the cross-correlation function is similar, where each cross-correlation function can be regarded as the data of a seismic trace. The array data of the cross-correlation function usually contains two dimensions. The first dimension represents the time axis, that is, the change of the signal values reflected by each cross-correlation function over time, and the second dimension represents the distance axis, that is, the offset information of the cross-correlation functions of different observation points relative to the virtual source point. The array data of the cross-correlation function, that is, the set of cross-correlation functions of the seismic signals of all seismic traces relative to the seismic signal of the virtual source, describes the propagation process of the reconstructed seismic waves.

[0055] In this embodiment, the seismic trace with the earliest arrival time of the seismic wave among all seismic traces is selected as the virtual source, and the cross-correlation calculation is performed on the observation data of other seismic traces relative to the virtual source. Finally, a surface wave dispersion extraction algorithm is used to extract the surface wave phase velocity dispersion curve from the array data of the cross-correlation function.

[0056] This cross-correlation method is relatively simple and does not require an accurate underground velocity model, but sufficient seismic data and relatively good signal-to-noise ratio are needed.

[0057] As Figure 3 shown, the figure shows the multi-trace seismic velocity field records observed by traditional seismic stations. In this embodiment, multi-trace seismic velocity field data is collected through numerical simulation. For example, an array with randomly distributed receiving points is used for observation, and it is assumed that the array contains 30 stations. Each station can observe seismic data of vertical (Z), radial (X), and tangential (Y) components (i.e., three-component seismic data). In this embodiment, the X-component seismic data is taken as an example. The seismic data of all seismic traces are arranged in ascending order of offset to form array data. It can be found that as the offset increases, the waveform curves continuously diverge and become more complex, reflecting the dispersion effect of surface waves.

[0058] As shown Figure 4 in the figure, it is a record of the multi-channel seismic strain rate field observed by DAS. In this embodiment, through numerical simulation, multi-channel seismic strain rate field data is collected. Here, taking the optical fiber in the X-axis direction of the DAS observation system as an example, the channel spacing is 0.15 km, with a total of 41 channels, arranged in ascending order along the X-axis, and the position of 4 km is the closest to the seismic source. It can be found that as the distance from 4 km increases, the waveform curves continuously diverge and become more complex, reflecting the dispersion effect of surface waves.

[0059] The dispersion effect of surface waves means that different frequency components of surface waves have different propagation speeds, and surface waves will have a wider wave packet as the propagation distance increases, that is, the waveform diverges. This makes it difficult to accurately extract the arrival time information of surface waves from the original observation data, so the traditional source location technology based on arrival time information is difficult to be directly applied to surface waves. Therefore, the present invention further processes the surface wave signals in the seismic observation data by using the obtained surface wave phase velocity dispersion curve, so as to analyze and extract the effective, stable and accurate source epicenter position.

[0060] Please further refer to Figure 1 , and perform step S200: Based on the surface wave phase velocity dispersion curve, perform surface wave dispersion backpropagation processing to extract the seismic waveform envelope curve.

[0061] Specifically, the seismic waveform curve refers to the curve of the amplitude of seismic waves changing with time, which describes the change of the intensity of seismic waves with time. A wave packet refers to a set of waves formed during the propagation of waves, which consists of many wave crests and wave troughs, and its shape and propagation direction are the same as the propagation direction of the waves, and it maintains a relatively stable waveform during propagation. Seismic waves include body waves and surface waves, and their wave packets also include the wave packets of body waves and surface waves.

[0062] The seismic waveform envelope curve is the curve connecting the wave crests or wave troughs in the seismic waveform curve, representing the overall shape of the seismic wave packet. The distribution of wave crests and wave troughs determines the shape of the envelope curve. When the wave crests and wave troughs are relatively dense and evenly distributed, the envelope curve will show a relatively smooth shape; while when the wave crests and wave troughs are relatively sparse or unevenly distributed, the envelope curve will show a relatively complex shape.

[0063] Surface waves are a special wave phenomenon, which is a wave form propagating on the surface of the medium. Different from body waves (such as longitudinal waves and transverse waves), for example, Rayleigh surface waves are formed by the vertical components of longitudinal waves and transverse waves interfering with each other on the free surface (such as the interface between air and solid), and Love surface waves are formed by the horizontal components of multiple reflected and refracted transverse waves interfering with each other on the free surface, and they propagate along the free surface, and the amplitude exponentially decays with the increase of depth.

[0064] During the propagation of surface waves, due to velocity dispersion, that is, surface waves with different frequency components propagate at different phase velocities. As the offset increases, the waveform of the surface wave will gradually diverge, the width of the surface wave packet will gradually increase, and the intensity of the surface wave signal will gradually weaken due to energy dispersion. Based on the obtained surface wave phase velocity dispersion curve c(ω), the reverse propagation factor of phase matching can be used to achieve the reverse transmission of surface wave dispersion, which can make the divergent surface wave waveforms in seismic data converge again.

[0065] The reverse propagation factor corresponds to the propagation factor. Both have the same functional form and are functions of the wave number, but the wave numbers of the two differ by a negative sign, which reflects that the propagation directions they control are opposite. After the source is excited, the generated surface waves propagate outward based on the propagation factor and at a specific phase velocity. By performing reverse transmission processing at the same phase velocity (i.e., phase matching) and based on the reverse propagation factor, the surface waves before propagation can be restored.

[0066] The specific implementation process is as follows:

[0067] Step S210: Obtain N-channel array data; calculate the corresponding seismic spectrum based on the N-channel array data.

[0068] Specifically, the observation systems established by seismic stations and DAS can both be used to receive the seismic waves excited by the source. Suppose a total of N-channel array seismic data are collected. In seismology, N-channel array data refers to the data of multiple seismic channels recorded simultaneously, usually an array composed of multiple sensors or seismic instruments. In this embodiment, the following operations are performed on each channel of seismic data. Suppose the nth channel of seismic data is S n (t), where n = 1, 2, …, N. Through Fourier transform (i.e., operation FFT), the seismic spectrum S n (ω) can be obtained:

[0069] S n (ω) = FFT[S n (t)]

[0070] Step S220: Perform surface wave dispersion reverse transmission processing based on the seismic spectrum to obtain the seismic spectrum after surface wave dispersion reverse transmission.

[0071] Specifically, based on the extracted surface wave phase velocity dispersion curve c(ω), through the reverse propagation factor of phase matching, the seismic spectrum S n (ω) is subjected to surface wave dispersion reverse transmission processing, and the seismic spectrum after surface wave dispersion reverse transmission can be obtained

[0072]

[0073] According to the theory of seismic plane wave propagation, the propagation factor and the back-propagation factor of a seismic plane wave can be defined by the exponential functions exp(-ik(ω)r) and exp(ik(ω)r) respectively. Among them, ω represents the angular frequency, r represents the offset, k(ω) = ω / c(ω) represents the radial wave number of the seismic plane wave, and c(ω) represents the phase velocity. It can be seen that the functional forms of the two are the same and both are functions of the wave number, but their wave numbers differ by a minus sign, which reflects that they control opposite propagation directions. After the source is excited, the generated surface wave propagates outward based on the propagation factor and at a specific phase velocity. To recover the surface wave before propagation, inverse propagation processing can be performed at the same phase velocity (i.e., phase matching) and based on the reverse propagation factor.

[0074] In addition, in the formula represents the offset, which is a function of the possible source coordinate points (x i , y j ) in the x-y coordinate grid (i = 1, 2,..., N x and j = 1, 2,..., N j ). The coordinate point (x n , y n ) represents the position of the nth receiving point.

[0075] Step S230: Process the seismic spectrum after surface wave dispersion inverse propagation to obtain the time-domain seismic waveform data after surface wave dispersion inverse propagation.

[0076] Specifically, it can be achieved by performing an inverse Fourier transform (i.e., operation IFFT) on the seismic spectrum after surface wave dispersion inverse propagation and taking the real part of the result of the inverse Fourier transform (i.e., operation Re), and calculating the time-domain seismic waveform data after surface wave dispersion inverse propagation

[0077] Exemplarily, the time-domain seismic waveform data after surface wave dispersion inverse propagation can be calculated by the following formula:

[0078]

[0079] Step S240: Calculate the envelope curve of the time-domain seismic waveform based on the time-domain seismic waveform data after surface wave dispersion inverse propagation; extract the surface wave packet with high time resolution based on the envelope curve of the time-domain seismic waveform.

[0080] Specifically, the seismic waveform envelope can describe the change of the amplitude of the seismic wave packet over time. A wave packet refers to a collection of waves formed by the propagation of waves in space. It consists of many peaks and troughs. Its shape and propagation direction are the same as the propagation direction of the wave, and it maintains a relatively stable waveform during the propagation process. The seismic wave packet contains the wave packets of body waves and surface waves.

[0081] The seismic waveform envelope curve is a curve connecting the peaks or troughs in the seismic waveform curve, which represents the overall shape of the seismic wave packet. The distribution of the peaks and troughs in the wave packet determines the shape of the envelope curve. When the peaks and troughs are relatively dense and evenly distributed, the envelope curve will present a smoother shape; when the peaks and troughs are relatively sparse or unevenly distributed, the envelope curve will present a more complex shape. The seismic waveform envelope curve also includes the envelope curve corresponding to the surface wave packet.

[0082] Hilbert transform (i.e., operation H) can be used to calculate the envelope curve of the seismic waveform in the time domain. Hilbert transform is widely used in signal processing. The physical meaning of this transform is to delay the phase of all frequency components of the signal by 90 degrees. It can be used as a demodulator to decode both amplitude modulation and frequency modulation. Hilbert transform is also often used to construct the complex analytical signal of the original signal. The construction form is the original signal plus the imaginary unit multiplied by the Hilbert transform of the original signal. Finally, by taking the modulus of the complex analytical signal, the envelope curve corresponding to the original signal can be obtained.

[0083] Through the above method, the envelope curve of the time domain seismic waveform can be obtained, which includes the envelope curve corresponding to the surface wave packet. Due to the surface wave dispersion backpropagation, the surface wave packet appears in the time window with the largest amplitude in the envelope curve. Compared with the original surface wave signal, its width will be greatly compressed and its amplitude will be significantly enhanced.

[0084] Exemplarily, the envelope curve of the time domain seismic waveform It can be calculated by the following formula:

[0085]

[0086] Of course, other ways of calculating waveform envelopes are also suitable for the purpose of this step, such as using other envelope detection algorithms, such as envelope detectors, to calculate the envelope curve of the seismic waveform.

[0087] The envelope curve is further filtered to highlight the surface wave signal, wherein a low-pass filter can be used to retain the lower frequency components. Detection of surface wave packets in the envelope curve can also be performed by identifying other characteristics of the wave packet, such as the amplitude, duration and frequency range of the wave packet.

[0088] Finally, an automatic detection algorithm, such as threshold detection or correlation detection, can be used to automatically determine the presence and location of the surface wave packet. It should be noted that the effect of extracting the surface wave packet may be affected by the quality of the seismic waveform and the noise level. Therefore, before extraction, it may be necessary to preprocess and optimize the waveform to obtain more optimized results.

[0089] Step S300: Based on the seismic waveform envelope curve, perform grid search and multi-channel stacking processing to obtain the source epicenter location.

[0090] Specifically, due to the conservation of surface wave energy, the width of the surface wave packet after dispersion backpropagation of the surface wave is greatly compressed, and the amplitude of the surface wave packet is significantly enhanced. By stacking the envelope curves corresponding to multi-channel seismic data and obtaining the maximum amplitude (i.e., the stacked amplitude) of the stacked envelope curve, and performing grid search to determine the grid coordinates that maximize the stacked amplitude, the source epicenter location can be finally determined.

[0091] Grid search is an exhaustive search method that finds the parameter combination that optimizes the mathematical model by traversing all possible combinations of the parameters to be determined. The present invention also adopts the idea of grid search for source location, and the parameter to be determined is the horizontal coordinate point of the source. The target area is grid-divided (i = 1, 2,..., N x and j = 1, 2,..., N y ), and all grid coordinate points are all possible combinations of the parameters to be determined. For each coordinate point (x i , y j ) as the possible source epicenter location, the following operations are performed in sequence.

[0092] The specific implementation process is as follows:

[0093] Step S310: Based on the envelope curve of the time-domain seismic waveform corresponding to the array seismic data, stack the envelope curves of the time-domain seismic waveforms of all channels to obtain the stacked envelope curve; based on the stacked envelope curve, obtain its maximum amplitude as the stacked amplitude corresponding to each coordinate point.

[0094] Specifically, through the steps of S200, the envelope curve of the time-domain seismic waveform corresponding to each grid coordinate point (x i , y j ) has been calculated To improve the reliability and stability of source location, all envelope curves corresponding to the array seismic data containing N channels can be stacked, and the maximum amplitude of the stacked envelope curve can be obtained to obtain the stacked amplitude E(x i , y j ) corresponding to each coordinate point (x i , y j)。

[0095] Exemplarily, for each coordinate point (x i , y j ), the superimposed amplitude E(x i , y j ) can be calculated by the following formula:

[0096]

[0097] When the coordinate point (x i , y j ) is closer to the true source position, the quality of the surface wave dispersion backpropagation is higher, the divergent surface wave waveforms converge better, the width of the surface wave packet is compressed narrower, the amplitude of the surface wave packet is stronger, and finally the superimposed amplitude E(x i , y j ) is larger; conversely, when the coordinate point (x i , y j ) is farther from the true source position, the superimposed amplitude E(x i , y j ) is smaller.

[0098] Step S320: Compare all the superimposed amplitudes, determine the coordinate point corresponding to the maximum value of the superimposed amplitude, and determine the coordinate point corresponding to the maximum value of the superimposed amplitude as the epicenter position of the earthquake source.

[0099] Furthermore, in order to further accurately determine the obtained epicenter position of the earthquake source, this method may further include:

[0100] Step S400: Iteratively update the epicenter position of the earthquake source.

[0101] Specifically, obtain the coordinates of all seismic receiving points, and based on the obtained epicenter position of the earthquake source, calculate the offset distances of the coordinates of all the seismic receiving points relative to the epicenter position of the earthquake source.

[0102] Based on the offset distances, directly extract the surface wave phase velocity dispersion curve from the array seismic data based on the N-channel array seismic data. Replace the surface wave phase velocity dispersion curve obtained in the original step S100 with the updated surface wave phase velocity dispersion curve.

[0103] Based on the updated surface wave phase velocity dispersion curve, repeat the operations of steps S200 - S300, and a grid search can be performed in the adjacent area of the obtained epicenter coordinate position of the earthquake source, and finally the updated epicenter coordinate position of the earthquake source can be obtained. The above steps can be iterated several times until, under the preset accuracy of the current grid division interval, the offset distance of the currently obtained epicenter position of the earthquake source compared to the epicenter position of the earthquake source obtained in the previous cycle is within the preset error range, then determine the currently obtained epicenter position of the earthquake source as the final epicenter position of the earthquake source.

[0104] Among them, the thresholds of the preset precision and the preset error range can be set to different values according to different events due to the different sizes of the research areas. In actual operation, they can be set artificially according to the precision requirements, such as meter level or kilometer level.

[0105] As Figure 5 shown, it is the source location effect diagram of the application of the present invention to the seismic station observation data. Based on the surface wave signals in the multi-channel seismic data, the source location technology proposed by the present invention is adopted to perform processing such as surface wave dispersion backpropagation and superposition. Finally, the epicenter position of the source is successfully determined. The black five-pointed star in the figure represents the real source position; the black triangles represent the layout positions of the seismic stations, forming a randomly distributed array observation system. The contour map represents the positioning result obtained by applying the present invention: among them, the position enclosed by the contour line with the largest value, that is, the pure white area enclosed by contour line 14, represents the inferred source position. The black five-pointed star represents the real source position, and the fact that the black five-pointed star falls into the position enclosed by the contour line with the largest value verifies the accuracy and effectiveness of the present invention.

[0106] As Figure 6 shown, it is the source location effect diagram of the application of the present invention to DAS observation data. According to the source location technical solution proposed by the present invention, operations such as surface wave dispersion backpropagation and superposition are performed, and finally the epicenter position of the source is also successfully determined. Among them, the black five-pointed star represents the real source position; the black straight line represents the fiber layout position, forming a DAS observation system. The contour map represents the positioning result obtained by applying this method. Among them, the position enclosed by the contour line with the largest value, that is, the pure white area enclosed by contour line 100, represents the inferred source position. The black five-pointed star, that is, the real source position, falls into the inferred source position, verifying the accuracy and effectiveness of the present invention.

[0107] Based on the above source location method based on surface wave dispersion backpropagation and superposition, this embodiment provides a source location system based on surface wave dispersion backpropagation and superposition, as Figure 7 shown, the system includes:

[0108] A surface wave phase velocity dispersion curve extraction module 51, configured to obtain a surface wave phase velocity dispersion curve;

[0109] An earthquake waveform envelope curve extraction module 52, configured to perform surface wave dispersion backpropagation processing based on the surface wave phase velocity dispersion curve and extract an earthquake waveform envelope curve;

[0110] A source epicenter position acquisition module 53, configured to perform grid search and multi-channel superposition processing based on the earthquake waveform envelope curve to obtain the source epicenter position.

[0111] In the solution of the present invention, a method and system for source location based on surface wave dispersion backpropagation and stacking are proposed. By obtaining the surface wave phase velocity dispersion curve, based on the surface wave phase velocity spectrum curve, surface wave dispersion backpropagation processing is carried out to extract the seismic waveform envelope curve. Based on the seismic waveform envelope curve, grid search and multi-channel stacking processing are carried out to obtain the source epicenter position. It solves the problem that the processing solution based on the arrival time information of body wave signals with relatively high frequencies in the prior art cannot be applied to the source location of surface waves, realizes the source epicenter location using surface wave signals with relatively low frequencies and higher signal-to-noise ratios, and ensures the accuracy and stability of the source location method of surface wave dispersion backpropagation and stacking. At the same time, it is applicable to traditional seismic observation data and new DAS observation data, expanding the application scope of source location technology.

[0112] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application, not to limit them. Although the present application has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some of the technical features. And these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for earthquake source location based on surface wave dispersion back propagation and superposition, characterized in that: The earthquake source location method based on surface wave dispersion back propagation and superposition includes: Obtain surface wave phase velocity dispersion curve; Based on the surface wave phase velocity dispersion curve, performing surface wave dispersion back-propagation processing to extract the seismic waveform envelope curve; Based on the seismic waveform envelope curve, grid search and multi-channel stacking processing are performed to obtain the epicenter position of the earthquake source; The step of performing surface wave dispersion back propagation processing based on the surface wave phase velocity dispersion curve to extract the seismic waveform envelope curve comprises: Acquire N-channel array seismic data; Calculate and obtain a corresponding seismic spectrum based on the N-channel array seismic data; Assume that the nth seismic data is S n (t), where n = 1, 2, ..., N, the earthquake spectrum S is obtained by Fourier transform n (ω):S n (ω) = FFT[S n (t)]; Performing surface wave dispersion backpropagation processing based on the seismic spectrum to obtain a seismic spectrum after surface wave dispersion backpropagation; Based on the extracted surface wave phase velocity dispersion curve c(ω), the seismic spectrum S n (ω) Performing surface wave dispersion back propagation processing to obtain the seismic spectrum after the surface wave dispersion back propagation Where ω represents the angular frequency, r represents the offset, k(ω)=ω / c(ω) represents the radial wave number of the seismic plane wave, and c(ω) represents the phase velocity; in the formula represents the offset, which is the possible source coordinate point (x i ,y j ), where i=1,2,…,N x and j = 1, 2, ..., N y ; Coordinate point (x n ,y n ) represents the position of the nth receiving point; Processing the seismic spectrum after the surface wave dispersion backpropagation to obtain time domain seismic waveform data after the surface wave dispersion backpropagation; Based on the time-domain seismic waveform data after the surface wave dispersion back-propagation, an envelope curve of the time-domain seismic waveform is calculated; The step of performing grid search and multi-channel stacking processing based on the seismic waveform envelope curve to obtain the epicenter position of the earthquake source also includes: Based on the envelope curve of the time-domain seismic waveform, obtaining the maximum amplitude of the envelope curve of the time-domain seismic waveform; Based on the array seismic data, the maximum amplitude of the envelope curve of the time-domain seismic waveform is superimposed to calculate the superimposed amplitude corresponding to each coordinate point; All the superimposed amplitudes are compared to determine the coordinate point corresponding to the maximum value of the superimposed amplitude, and the coordinate point corresponding to the maximum value of the superimposed amplitude is determined as the epicenter position of the earthquake source.

2. The earthquake source location method based on surface wave dispersion back propagation and superposition according to claim 1 is characterized in that: The step of obtaining the surface wave phase velocity dispersion curve comprises: The surface wave phase velocity dispersion curve is obtained by forward modeling based on the public velocity model and the dominant frequency band of the observation data.

3. The earthquake source location method based on surface wave dispersion back propagation and superposition according to claim 1, characterized in that: The step of obtaining the surface wave phase velocity dispersion curve comprises: Get the seismic trace with the smallest offset and determine it as the virtual source; Performing cross-correlation calculations relative to the virtual source based on the remaining seismic channel observation data to obtain array data of the cross-correlation function; The surface wave phase velocity dispersion curve is extracted from the array data of the cross-correlation function using a surface wave dispersion extraction algorithm.

4. A source location system based on surface wave dispersion back propagation and superposition using the method described in any one of claims 1 to 3, characterized in that: The system comprises: A surface wave phase velocity dispersion curve extraction module is used to obtain a surface wave phase velocity dispersion curve; A seismic waveform envelope curve extraction module, used to perform surface wave dispersion back-propagation processing based on the surface wave phase velocity dispersion curve to extract the seismic waveform envelope curve; The module for obtaining the epicenter position of the earthquake source is used to perform grid search and multi-channel stacking processing based on the earthquake waveform envelope curve to obtain the epicenter position of the earthquake source.

Images