Method, device, computer equipment, and storage medium for determining phase velocity of seismic waves

Abstract

The present invention relates to a method, device, computer equipment and storage medium for determining the phase velocity of seismic waves, the method comprising: obtaining seismic data collected by a monitoring point in multiple time windows; determining the reverse azimuth from the earthquake source to the monitoring point; for each time window, determining the lateral acceleration according to the reverse azimuth and the north-south translation component and the east-west translation component; for each time window, determining the preliminary phase velocity of the time window according to the lateral acceleration and the vertical rotation component; determining the final phase velocity of the seismic wave according to the preliminary phase velocities of multiple time windows. The above method can determine the phase velocity of the seismic wave more accurately.

Description

Technical Field

[0001] The present invention relates to the field of earthquake monitoring technology, and in particular to a method, device, computer equipment and storage medium for determining the phase velocity of earthquake waves. Background Art

[0002] Seismic wave velocity is the most important parameter in seismic exploration. The propagation velocity of seismic waves can be understood as the distance that seismic waves propagate along the wave line in a unit time. The wave velocity is related to the type of wave and the medium. In the process of seismic data processing and interpretation, there are many places where velocity data is needed. For example, the average velocity in the velocity data is used for time-depth conversion or structural mapping. The so-called average velocity can be understood as the ratio of the total thickness of seismic waves passing through these interfaces vertically in a set of horizontal layered media to the time taken. The root mean square velocity is a velocity concept introduced when the time-distance curve of the reflection wave of the horizontal layered medium is simplified to a hyperbolic relationship, which is mainly used in horizontal stacking and migration. For a certain phase axis of a set of common reflection point gathers, a series of different velocities are selected using the hyperbolic formula to calculate the dynamic correction of each channel to perform dynamic correction on each channel in the gather. When a certain velocity can correct the phase axis into a horizontal straight line, this velocity is the stacking velocity required for calculating the dynamic correction. For homogeneous media with inclined interfaces, the stacking velocity is the equivalent velocity; for horizontal layered media, the stacking velocity is the root mean square velocity. In seismic exploration, the wave velocity of a certain stratum is called the layer velocity of the stratum, and the layer velocity is needed in the geological interpretation of seismic data. Theoretical research and actual data have confirmed that the propagation velocity of seismic waves in rock strata is related to factors such as the properties of the rock strata, the composition of the rock, density, burial depth, geological age, porosity, and fluid properties (Lu Jimeng. Principles of Seismic Exploration [M]. Petroleum University Press, 2009.). Extracting and analyzing velocity data is an important task in geological interpretation. It is extremely important to be familiar with various velocities and their application scenarios as well as the solution of velocity data.

[0003] The regular characteristics of the elastic wave field are essentially the regular characteristics of the velocity field. Through the analysis and study of time and space, kinematics and dynamics in the velocity field, the distribution characteristics and laws of the seismic wave field can be studied, and then the structure, structure and properties of the medium can be understood. The speed of elastic wave propagation is mainly divided into phase velocity and group velocity. Among them, the speed at which the equal phase surface (such as the crest surface or the trough surface) of a simple harmonic wave of a single frequency propagates in the medium is called the phase velocity. The speed at which the composite signal of many waves of different frequencies propagates in the medium is called the group velocity. Waves of different frequencies have different amplitudes and phases. The group velocity is the propagation speed of any constant phase point on the envelope wave, and it is a propagation speed representing energy. In a uniform isotropic medium, the phase velocity does not change with the propagation direction, the group velocity is consistent with the phase velocity, and the phase velocity surface and the group velocity surface coincide; in anisotropic media, the phase velocity changes with the propagation direction, the phase velocity and the group velocity no longer coincide, and they are equal only in a specific direction. Phase velocity is one of the parameters describing the physical properties of the medium. It is reflected by the physical property matrix of the medium. At the same time, the elastic matrix of the medium is implicit in the wave equation, that is, the velocity defined by the wave equation is the phase velocity. Combining the constitutive equation of elastic dynamics, Newton's differential equation of motion and geometric equation, the elastic wave equation and the Christoffel equation are obtained. The Christoffel equation is solved by the eigenvalue method, and the analytical expression of the elastic wave phase velocity is derived (Wu Guochen, Liang Kai, Qi Yanping. Phase velocity and group velocity of three-dimensional TTI medium [J]. Progress in Geophysics, 2009, 24 (006): 2097-2105.). Since the exact formula for solving the Christoffel equation is often complicated, it is difficult to use. Considering the common weak anisotropy of underground media, the weak anisotropy approximation assumption can be used to simplify the precise phase velocity formula of elastic waves. However, the weak anisotropy approximation may have large errors when facing moderately strong anisotropic media. Therefore, other approximate methods can be used to approximately characterize the elastic waves of moderately strong anisotropic media (Sun Shangrao, Liang Kai. General approximation of elastic wave phase velocity in three-dimensional TTI media [C] / / 2019 China Earth Science Joint Academic Annual Conference. 0.).

[0004] Phase velocity and group velocity are important parameters for studying the laws of seismic wave propagation and describing medium characteristics. They are the core content of elastic wave propagation theory and play an important role in theoretical research and practical applications. Summary of the invention

[0005] The present application provides a method, apparatus, computer equipment and storage medium for determining the phase velocity of a seismic wave.

[0006] A first aspect provides a method for determining the phase velocity of a seismic wave, the method comprising:

[0007] Acquire seismic data collected by the monitoring point in multiple time windows, the seismic data including: a north-south translation component in the north-south direction, an east-west translation component in the east-west direction, and a vertical rotation component perpendicular to a plane formed by the north-south direction and the east-west direction;

[0008] Determine the back azimuth from the earthquake source to the monitoring point;

[0009] For each time window, determining a lateral acceleration according to the back azimuth and the north-south translational component and the east-west translational component;

[0010] For each time window, determining a preliminary phase velocity of the time window according to the lateral acceleration and the vertical rotation component;

[0011] The final phase velocity of the seismic wave is determined according to the preliminary phase velocities of a plurality of time windows.

[0012] In some embodiments, determining the back azimuth from the earthquake source to the monitoring point includes:

[0013] Determine a plurality of candidate back azimuths from the earthquake source to the monitoring point;

[0014] For each time window, determining the lateral acceleration corresponding to each candidate back azimuth angle according to the north-south translation component, the east-west translation component and each candidate back azimuth angle;

[0015] determining, for each time window, a zero-lag cross-correlation coefficient between the lateral acceleration and the vertical rotation component;

[0016] For each time window, determining a preliminary back azimuth according to the candidate back azimuth and the zero-lag cross-correlation coefficient;

[0017] The back azimuth from the earthquake source to the monitoring point is determined based on the preliminary back azimuths of multiple time windows.

[0018] In some embodiments, determining a plurality of candidate back azimuth angles from the earthquake source to the monitoring point comprises:

[0019] The back azimuth angles from 0° to 360° are divided into multiple parts at preset degree intervals to obtain multiple candidate back azimuth angles.

[0020] In some embodiments, determining a preliminary back azimuth angle according to the candidate back azimuth angle and the zero-lag cross-correlation coefficient comprises:

[0021] The candidate back azimuth angle corresponding to the largest zero-lag cross-correlation coefficient is selected as the preliminary back azimuth angle of the time window.

[0022] In some embodiments, determining the preliminary phase velocity of the time window according to the lateral acceleration and the vertical rotation component comprises:

[0023] Determining the maximum-amplitude lateral acceleration and the maximum-amplitude vertical rotation component in the time window;

[0024] The preliminary phase velocity of the time window is determined according to the lateral acceleration with the largest amplitude and the vertical rotation component with the largest amplitude.

[0025] In some embodiments, after determining the lateral acceleration of each time window and before determining the preliminary phase velocity of the time window according to the phase velocities of the multiple time windows, the method further includes:

[0026] determining, for each time window, a zero-lag cross-correlation coefficient between the lateral acceleration and the vertical rotation component;

[0027] The phase velocity of the time window that meets the preset condition is filtered out, and the preset condition is: the zero-lag mutual correlation coefficient of the time window is less than a preset threshold.

[0028] In some embodiments, determining a zero-lag correlation coefficient between the lateral acceleration and the vertical rotation component comprises:

[0029] Determine the mean value of the lateral acceleration and the mean value of the vertical rotation component;

[0030] A zero-lag correlation coefficient between the lateral acceleration and the vertical rotation component is determined based on the lateral acceleration, the vertical rotation component, the average value of the lateral acceleration, and the average value of the vertical rotation component.

[0031] A second aspect provides a device for determining the phase velocity of a seismic wave, comprising:

[0032] A data acquisition unit is used to acquire seismic data collected by the monitoring point in multiple time windows, wherein the seismic data includes: a north-south translation component in the north-south direction, an east-west translation component in the east-west direction, and a vertical rotation component perpendicular to a plane formed by the north-south direction and the east-west direction;

[0033] A back azimuth angle determination unit, used to determine the back azimuth angle from the earthquake source to the monitoring point;

[0034] A lateral acceleration calculation unit, configured to determine the lateral acceleration according to the back azimuth angle and the north-south translational component and the east-west translational component for each time window;

[0035] a phase velocity calculation unit, configured to determine, for each time window, a preliminary phase velocity of the time window according to the lateral acceleration and the vertical rotation component;

[0036] The result output unit is used to determine the final phase velocity of the seismic wave according to the preliminary phase velocities of multiple time windows.

[0037] The third aspect provides a computer device, including a memory and a processor, wherein the memory stores computer-readable instructions, and when the computer-readable instructions are executed by the processor, the processor executes the steps of the above-mentioned method for determining the phase velocity of seismic waves.

[0038] A fourth aspect provides a storage medium storing computer-readable instructions, which, when executed by one or more processors, causes the one or more processors to execute the steps of the above-mentioned method for determining the phase velocity of seismic waves.

[0039] The above-mentioned method, device, computer equipment and storage medium for determining the phase velocity of seismic waves determine the reverse azimuth from the earthquake source to the monitoring point; obtain the north-south translational acceleration, the east-west translational acceleration and the vertical rotational motion component collected by the monitoring point, and determine the lateral acceleration according to the north-south translational acceleration and the east-west translational acceleration; divide the seismic data collected in the historical time period according to the preset time window; determine the phase velocity of the time window according to the lateral acceleration and the vertical rotational component of each time window; fit the phase velocity line of multiple time windows to obtain the phase velocity of the seismic wave. Therefore, the phase velocity solution is realized based on the six-component seismic data actually collected by a single station. Compared with the traditional array method for solving the rotation component, there is no site restriction and the cost of manual layout is greatly reduced. Compared with the traditional approximate derivation based on the wave equation, the method of solving the phase velocity based on the combination of the single rotation component and the translation component is more convenient in practical application, and the similar information obtained is quite stable. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 A flow chart of a method for determining the phase velocity of seismic waves in one embodiment of the present application;

[0041] Figure 2 This is a flow chart of an algorithm for solving phase velocity using six-component data of a single station in a method for determining phase velocity of seismic waves in one embodiment of the present application;

[0042] Figure 3 The waveform of the lateral acceleration and the actual measured rotation component of the method for determining the phase velocity of seismic waves in one embodiment of the present application is shown in the case where the back azimuth angle is 280°; the abscissa corresponds to time (unit, seconds), and the ordinate corresponds to the normalized amplitude; the solid line represents the lateral acceleration, and the dotted line represents the vertical rotation component;

[0043] Figure 4The lateral acceleration and the actual measured rotation component waveform of the phase velocity determination method of seismic waves in one embodiment of the present application; the abscissa corresponds to time (unit, seconds), and the ordinate corresponds to the normalized amplitude; the solid line represents the lateral acceleration, the dotted line represents the vertical rotation component, and the horizontal dotted line is the dividing line when the correlation threshold is 0.75;

[0044] Figure 5 This is a structural block diagram of a device for determining the phase velocity of seismic waves in one embodiment of the present application;

[0045] Figure 6 FIG. 4 is a block diagram of the internal structure of a computer device in one embodiment. DETAILED DESCRIPTION

[0046] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0047] It can be understood that the terms "first", "second", etc. used in the present application can be used in this article to describe various elements, but these elements are not limited by these terms. These terms are only used to distinguish the first element from another element.

[0048] An implementation environment diagram of an earthquake source monitoring method provided in one embodiment, in which the implementation environment may include computer equipment and earthquake monitoring instruments.

[0049] The computer device has an interface, for example, the interface may be an API (Application Programming Interface). The earthquake monitor is a device set at a detection point to collect earthquake data, and the computer device performs earthquake source monitoring (such as determining the reverse azimuth of the earthquake source) based on the data collected by the earthquake monitor.

[0050] It should be noted that the terminal and the computer device may be a smart phone, a tablet computer, a laptop computer, a desktop computer, etc., but are not limited thereto. The computer device and the terminal may be connected via Bluetooth, USB (Universal Serial Bus) or other communication connection methods, and the present invention does not limit this.

[0051] In one embodiment, a method for determining the phase velocity of seismic waves is proposed. The method for determining the phase velocity of seismic waves can be applied to the above-mentioned computer device. The method is implemented based on a python program; Figure 1 As shown, the following steps may be specifically included:

[0052] Step 101, obtaining seismic data collected by monitoring points in multiple time windows, the seismic data including: a north-south translation component in the north-south direction, an east-west translation component in the east-west direction, and a vertical rotation component perpendicular to a plane formed by the north-south direction and the east-west direction.

[0053] It can be understood that the earthquake data collected by the monitoring point within a preset historical time period is obtained, and the historical time end is divided into multiple time windows. For example, the earthquake data within 12 seconds before the current moment is obtained, and the above earthquake data can be divided into 3 time windows (each time window is 4 seconds long).

[0054] In this step, the north-south direction and the east-west direction are two directions perpendicular to each other in the three-dimensional space, and the translational components and rotational components of the three coordinate axes of the three-dimensional space collected by the monitoring point are obtained, including:

[0055] A six-component seismograph is used to collect translational components and rotational components in three directions (mutually orthogonal X-axis, Y-axis and Z-axis) of a three-dimensional space; the six-component seismograph may include: three fiber optic gyroscopes, three accelerometers and a six-component signal processing module, wherein the sensitive axes of the three fiber optic gyroscopes are mutually orthogonal, and the sensitive axes of the three accelerometers are mutually orthogonal; the sensitive axis of each fiber optic gyroscope is parallel to or coincides with the three sensitive axes of an accelerometer respectively; the input end of the six-component signal processing module is connected to the detection signal output end of each fiber optic gyroscope and the detection signal output end of each accelerometer, the six-component signal processing module is used to generate a modulation signal required by the fiber optic gyroscope, and obtain a detection angular velocity according to a detection signal output from the detection signal output end of the fiber optic gyroscope, and perform error compensation on the detection angular velocity; the six-component signal processing module is also used to obtain a detection translation acceleration according to a detection signal output from the detection signal output end of the three-axis accelerometer, and perform error compensation on the detection translation acceleration. A spatial rectangular coordinate system is established, which includes mutually orthogonal X-axis, Y-axis and Z-axis, and the sensitive axes of the fiber optic gyroscope and the accelerometer are on the X-axis, Y-axis and Z-axis respectively.

[0056] Furthermore, the six-component seismograph includes three orthogonal fiber optic gyroscopes and three equally orthogonal accelerometers. The fiber optic gyroscopes are used to measure the angular velocity information of the object, and the three-axis accelerometers are used to measure the linear motion information of the carrier. The six-component signal processing module can provide complete seismic wave field information by processing the outputs of the fiber optic gyroscopes and accelerometers.

[0057] Of course, seismic data can include: the north-south translation component and the north-south rotation component in the north-south direction; the east-west translation component and the east-west rotation component in the east-west direction, and the vertical translation component and the vertical rotation component perpendicular to the plane formed by the north-south and east-west directions.

[0058] Step 102, determining the back azimuth from the earthquake source to the monitoring point;

[0059] In this step, the back azimuth from the earthquake source to the monitoring point may be known in advance. If it is unknown, the above step 102 may include:

[0060] Step 1021, determining multiple candidate back azimuths from the earthquake source to the monitoring point;

[0061] In this step, the candidate back azimuth angle may be determined empirically or based on seismic data.

[0062] In some embodiments, the cross-correlation method is used to determine the final back azimuth of the earthquake source. In the cross-correlation method, in step 102, setting the candidate back azimuth based on experience may include:

[0063] Step 102A1: Divide the back azimuth angle from 0° to 360° into multiple parts at preset degree intervals to obtain multiple candidate back azimuth angles.

[0064] In some embodiments, a polarization analysis method is used to determine the final back azimuth of the earthquake source. In the polarization analysis method, determining the candidate back azimuth according to the seismic data in step 102 may include:

[0065] Step 102B1, obtaining an estimated back azimuth angle according to the first rotation component with the largest amplitude and the second rotation component with the largest amplitude;

[0066] In this step, the first rotation component with the largest amplitude and the second rotation component with the largest amplitude are substituted into the calculation formula for estimating the back azimuth angle to obtain the estimated back azimuth angle, wherein the calculation formula for estimating the back azimuth angle is:

[0067]

[0068] In the formula, θ BAz To estimate the back azimuth; R e is the first rotation component with the largest amplitude; R n is the second rotation component with the largest amplitude.

[0069] Step 102B2, determine whether the estimated back azimuth angle is less than 0. If the estimated back azimuth angle is not less than 0, use the estimated back azimuth angle as the candidate back azimuth angle; if the estimated back azimuth angle is less than 0, add 180° to the estimated back azimuth angle to obtain the candidate back azimuth angle.

[0070] Step 1022: for each time window, determine the lateral acceleration corresponding to each candidate back azimuth angle according to the north-south translation component, the east-west translation component and each candidate back azimuth angle;

[0071] In this step, the north-south translational component, the east-west translational component and each candidate back azimuth angle are substituted into the lateral acceleration formula to obtain the lateral acceleration of each candidate back azimuth angle.

[0072] The lateral acceleration calculation formula is:

[0073]

[0074] In the formula, A e A is the north-south translational component; n is the east-west translational component; θ Baz is the candidate back azimuth; A t is the lateral acceleration.

[0075] It can be understood that the six-component seismograph collects seismic data at preset time intervals (such as once per second), and the lateral acceleration (A) corresponding to each collection moment is t ) is obtained by substituting the north-south translational component, east-west translational component, and back azimuth (candidate back azimuth or final back azimuth) collected at each collection moment into the lateral acceleration calculation formula.

[0076] Step 1023, for each time window, determine the zero-lag correlation coefficient between the lateral acceleration and the vertical rotation component;

[0077] In this step, each lateral acceleration and vertical rotation component is substituted into the relevant calculation formula to obtain the zero-lag mutual correlation coefficient between the lateral acceleration and the vertical rotation component. The relevant calculation formula is:

[0078]

[0079] In the formula, A t The lateral acceleration corresponding to each acquisition moment; R z It is the vertical rotation component collected at each collection moment; is the average value of lateral acceleration; is the average value of the vertical rotation component; xcorr(A t , R z ) is the zero-lag correlation coefficient between the lateral acceleration and the vertical rotation components.

[0080] It can be understood that the average value of the lateral acceleration is the average value of the lateral acceleration in a time window (such as 3 seconds); the average value of the vertical rotation component is the average value of the vertical rotation component in a time window.

[0081] Step 1024: for each time window, determine a preliminary back azimuth according to the candidate back azimuth and the zero-lag cross-correlation coefficient;

[0082] In this step, there is a certain relationship between the zero-lag mutual correlation coefficient between the lateral acceleration and the vertical rotation component and the back-square angle of the earthquake source. Therefore, the final back-azimuth angle from the earthquake source to the monitoring point can be determined based on the candidate back-azimuth angle and the zero-lag mutual correlation coefficient.

[0083] In some embodiments, the cross-correlation method is used to determine the final back azimuth of the earthquake source. In the cross-correlation method, in step 1024, determining the preliminary back azimuth according to the candidate back azimuth and the zero-lag cross-correlation coefficient may include:

[0084] Step 1024A1: Select the candidate back azimuth angle corresponding to the largest zero-lag mutual correlation coefficient as the final back azimuth angle.

[0085] In some embodiments, a polarization analysis method is used to determine the final back azimuth of the earthquake source. In the polarization analysis method, in step 1024, determining the preliminary back azimuth according to the candidate back azimuth and the zero-lag cross-correlation coefficient may include:

[0086] Step 1024B1, determine whether the zero-lag mutual correlation coefficient is greater than 0. If it is greater than 0, add 180° to the candidate back azimuth angle to obtain the final back azimuth angle; if it is less than or equal to 0, the candidate back azimuth angle is the final back azimuth angle.

[0087] Step 1025: Determine the back azimuth from the earthquake source to the monitoring point based on the preliminary back azimuths of multiple time windows.

[0088] It is understandable that the back azimuth from the earthquake source to the monitoring point is determined based on the preliminary back azimuths of multiple time windows. That is, the preliminary back azimuths are mostly distributed near the theoretical back azimuths, and the results of individual time windows have large errors, which may be affected by other external factors. Overall, the prediction results are consistent with the theoretical values. Therefore, the preliminary back azimuths of multiple time windows are linearly fitted to obtain the back azimuth from the earthquake source to the monitoring point.

[0089] like Figures 2 to 4 In one application scenario, the feasibility of the cross-correlation method was tested using data from a magnitude 6.4 earthquake in Dali, Yunnan, collected by Wuhan instruments on May 21, 2021. The theoretical back azimuth is 280°, and the source (epicenter) is about 1583km away from the test station.

[0090] like Figure 2 As shown, the above method may include:

[0091] Step A1, obtain three translational components and rotational component seismic data from the rotational seismograph, where the length of the recorded data is 390s, starting from the time of the earthquake. Then perform preprocessing on the data such as removing instrument response;

[0092] Step A2: Divide the back azimuth angle from 0° to 360° into 72 parts at intervals of 5°. For each given back azimuth angle, calculate the first horizontal component and the second horizontal component according to the lateral acceleration formula to obtain the lateral acceleration.

[0093] Step A3: Divide the recorded 390s seismic data into 13 parts evenly with 30s as a time window. In each time window, loop through all candidate back azimuths and each time window. t and the vertical rotation component R z Substituting into the relevant calculation formula, the zero-lag mutual correlation coefficient between the lateral acceleration and the vertical rotation component.

[0094] Step A4: Select the back azimuth angle corresponding to the maximum position of the mutual correlation coefficient as the final result of the time window.

[0095] When all time window cycles are completed, the reverse azimuth prediction result of the entire time channel can be obtained.

[0096] In one application scenario, the feasibility of the polarization analysis method was tested using a methane explosion experiment conducted in Guangzhou on January 12, 2021. The theoretical back azimuth angle was 190°, and the distance between the methane source and the test station was about 305m. The above method may include:

[0097] Step B1, download three translational components and rotational component seismic data from the seismograph, where the length of the recorded data is 15s, starting from the explosion of the methane source. Then pre-process the downloaded data by removing the instrument response;

[0098] Step B2: The recorded 15s seismic data is evenly divided into 15 parts with 1s as a time window. In each time window, the first rotation component and the second rotation component with the largest amplitude are selected to calculate the candidate back azimuth.

[0099] Step B3: In each 1 s time window, the horizontal component is converted into lateral acceleration according to the calculated candidate back azimuth angle.

[0100] Step B4: Substitute the north-south translational component, the east-west translational component and the candidate back azimuth angles into the lateral acceleration formula to obtain the lateral acceleration of each candidate back azimuth angle.

[0101] Step B5: If the lateral acceleration is positively correlated with the vertical rotation component, the obtained reverse azimuth angle is increased by 180°, otherwise it remains unchanged.

[0102] Step B6, looping each time window to obtain the back azimuth prediction result of the entire time channel. The predicted back azimuth is relatively evenly distributed near the theoretical azimuth, but some time windows are affected by other external factors, resulting in large errors.

[0103] Step 103: for each time window, determine the lateral acceleration according to the back azimuth and the north-south translational component and the east-west translational component;

[0104] In this step, the north-south translational component, the east-west translational component and the back azimuth angle are substituted into the lateral acceleration formula to obtain the back azimuth angle lateral acceleration.

[0105] It can be understood that the north-south translational component and the east-west translational component collected at each collection moment of the detection point are substituted into the lateral acceleration formula together with the back azimuth angle to obtain the lateral acceleration corresponding to each collection moment.

[0106] Step 104: for each time window, determine the preliminary phase velocity of the time window according to the lateral acceleration and the vertical rotation component;

[0107] In this step, the phase velocity is obtained based on the lateral acceleration and vertical rotation components, so that the phase velocity can be solved based on the six-component seismic data actually collected by a single station. Compared with the traditional array method for solving the rotation component, there is no site restriction and the cost of manual layout is greatly reduced. Compared with the traditional approximate derivation based on the wave equation, the method of solving the phase velocity based on the combination of the single rotation component and the translation component is more convenient in practical applications, and the similar information obtained is quite stable.

[0108] In some embodiments, the above step 104 may include:

[0109] Step 1041, determining the lateral acceleration with the largest amplitude and the vertical rotation component with the largest amplitude in the time window;

[0110] Among them, in each time window, the maximum lateral acceleration and the maximum vertical rotation component are found. The two values ​​may not be collected at the same collection time, as long as they are the maximum values.

[0111] Step 1042: Determine the preliminary phase velocity of the time window according to the lateral acceleration with the largest amplitude and the vertical rotation component with the largest amplitude.

[0112] Specifically, the lateral acceleration and vertical rotation components with the largest amplitude are substituted into the phase velocity formula to obtain the phase velocity, which is:

[0113]

[0114] In the formula, A t is the lateral acceleration; Rz is the vertical rotation component; c is the phase velocity (initial acceleration).

[0115] Step 105: Determine the final phase velocity of the seismic wave according to the preliminary phase velocities of multiple time windows.

[0116] Among them, the phase velocity lines of multiple time windows are fitted to obtain the final phase velocity of the seismic wave.

[0117] Generally, the preliminary phase velocity of each time window is evenly distributed near the theoretical phase velocity. However, some time windows may have large errors due to the influence of other external factors. Therefore, the preliminary phase velocity lines of multiple time windows are fitted, and the preliminary phase velocities with large errors are eliminated to obtain the final phase velocity of the seismic wave.

[0118] In an application scenario, such as Figure 3 , 4 As shown, the feasibility of this application is tested using the data of the 6.4 magnitude earthquake in Dali, Yunnan collected by Wuhan instruments on May 21, 2021. Its theoretical back azimuth is 280°, and the distance between the epicenter and the test station is about 1583km. (Angle between the connecting line and the true north direction)

[0119] A. Based on the back azimuth angle of 280°, convert the north-south horizontal component and the east-west horizontal component into lateral acceleration;

[0120] The recorded 550 s seismic data were divided into 22 parts with a time window of 25 s.

[0121] For each time window, the zero-lag correlation coefficient between the lateral acceleration and the vertical rotation component is calculated, as Figure 4 shown.

[0122] The time window was selected based on a threshold of 0.75 correlation between the lateral acceleration component and the vertical rotation rate.

[0123] In each filtered time window, the phase velocity of the time window is calculated.

[0124] In the present application, the method for solving the phase velocity based on the combination of a single rotational component and a translational component is more convenient in practical applications, and the similar information obtained is quite stable.

[0125] In some embodiments, after step 103 and before step 104, the following may be included:

[0126] Step 103a, for each time window, determining the zero-lag correlation coefficient between the lateral acceleration and the vertical rotation component;

[0127] The method for determining the zero-lag mutual correlation coefficient in this step is the same as described above and will not be repeated here.

[0128] Step 103b: filtering out the phase velocities of the time windows that meet a preset condition, wherein the preset condition is that the zero-lag mutual correlation coefficient of the time window is less than a preset threshold.

[0129] In this step, for each time window, the zero-lag mutual correlation coefficient between the lateral acceleration and the vertical rotation component is calculated, and the time window is screened according to the zero-lag mutual correlation coefficient. When the zero-lag mutual correlation coefficient is less than the preset correlation coefficient threshold, the seismic data of the time window will be eliminated. For example, when the zero-lag mutual correlation coefficient of a time window is less than 0.75, the seismic data of the time window will be eliminated, that is, the preliminary phase velocity of the time window will not be calculated.

[0130] Of course, the zero-lag mutual correlation coefficient of this step may be the zero-lag mutual correlation coefficient of each time window obtained when determining the back azimuth from the earthquake source to the monitoring point in the above step 102.

[0131] In one embodiment, a seismic wave phase velocity determination device is provided, and the seismic wave phase velocity determination device can be used to run the seismic wave phase velocity determination method in the above embodiment, such as Figure 5 As shown, the device may specifically include:

[0132] The data acquisition unit 511 is used to acquire seismic data collected by the monitoring point in multiple time windows, the seismic data including: a north-south translation component in the north-south direction, an east-west translation component in the east-west direction, and a vertical rotation component perpendicular to the plane formed by the north-south direction and the east-west direction;

[0133] A back azimuth angle determination unit 512, used to determine the back azimuth angle from the earthquake source to the monitoring point;

[0134] A lateral acceleration calculation unit 513, for determining the lateral acceleration according to the back azimuth and the north-south translational component and the east-west translational component for each time window;

[0135] A phase velocity calculation unit 514 is used to determine the preliminary phase velocity of each time window according to the lateral acceleration and the vertical rotation component;

[0136] The result output unit 515 is used to determine the final phase velocity of the seismic wave according to the preliminary phase velocities of multiple time windows.

[0137] like Figure 6As shown, in one embodiment, a computer device is proposed, which may include a processor, a storage medium, a memory, and a network API interface connected via a system bus. Among them, the storage medium of the computer device stores an operating system, a database, and computer-readable instructions, and the database may store a control information sequence. When the computer-readable instructions are executed by the processor, the processor can implement an earthquake source monitoring method. The processor of the computer device is used to provide computing and control capabilities to support the operation of the entire computer device. The memory of the computer device may store computer-readable instructions. When the computer-readable instructions are executed by the processor, the processor can execute an earthquake source monitoring method. The network API interface of the computer device is used to connect and communicate with a terminal. Those skilled in the art can understand that Figure 6 The structure shown in the figure is only a block diagram of a part of the structure related to the solution of the present application, and does not constitute a limitation on the computer device to which the solution of the present application is applied. The specific computer device may include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components.

[0138] The processor implements the following steps when executing the computer program: obtaining seismic data collected by the monitoring point in multiple time windows, the seismic data including: a north-south translation component in the north-south direction, an east-west translation component in the east-west direction, and a vertical rotation component perpendicular to the plane formed by the north-south direction and the east-west direction; determining the back azimuth from the earthquake source to the monitoring point; for each time window, determining the lateral acceleration based on the back azimuth and the north-south translation component and the east-west translation component; for each time window, determining the preliminary phase velocity of the time window based on the lateral acceleration and the vertical rotation component; and determining the final phase velocity of the seismic wave based on the preliminary phase velocities of multiple time windows.

[0139] In one embodiment, a storage medium storing computer-readable instructions is proposed. When the computer-readable instructions are executed by one or more processors, the one or more processors perform the following steps: obtain seismic data collected by a monitoring point within multiple time windows, the seismic data including: a north-south translation component in the north-south direction, an east-west translation component in the east-west direction, and a vertical rotation component perpendicular to a plane formed by the north-south direction and the east-west direction; determine a back azimuth from the earthquake source to the monitoring point; for each time window, determine a lateral acceleration based on the back azimuth and the north-south translation component and the east-west translation component; for each time window, determine a preliminary phase velocity of the time window based on the lateral acceleration and the vertical rotation component; and determine a final phase velocity of the seismic wave based on the preliminary phase velocities of multiple time windows.

[0140] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiments can be implemented by instructing the relevant hardware through a computer program, and the computer program can be stored in a computer-readable storage medium. When the program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, the aforementioned storage medium can be a non-volatile storage medium such as a disk, an optical disk, a read-only memory (ROM), or a random access memory (RAM).

[0141] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0142] The above embodiments only express several implementation methods of the present invention, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the present invention. It should be pointed out that, for a person of ordinary skill in the art, several modifications and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention patent shall be subject to the attached claims.

Claims

1. A method for determining the phase velocity of seismic waves, characterized in that: The method comprises: Acquire seismic data collected by the monitoring point in multiple time windows, the seismic data including: a north-south translation component in the north-south direction, an east-west translation component in the east-west direction, and a vertical rotation component perpendicular to a plane formed by the north-south direction and the east-west direction; Determine the back azimuth from the earthquake source to the monitoring point; For each time window, determining a lateral acceleration according to the back azimuth and the north-south translational component and the east-west translational component; For each time window, determining a preliminary phase velocity of the time window according to the lateral acceleration and the vertical rotation component; The final phase velocity of the seismic wave is determined according to the preliminary phase velocities of a plurality of time windows.

2. The method for determining the phase velocity of seismic waves according to claim 1, characterized in that: The determining of the back azimuth from the earthquake source to the monitoring point comprises: Determine a plurality of candidate back azimuths from the earthquake source to the monitoring point; For each time window, determining the lateral acceleration corresponding to each candidate back azimuth angle according to the north-south translation component, the east-west translation component and each candidate back azimuth angle; For each time window, determining a zero-lag correlation coefficient between each of the lateral acceleration and the vertical rotation component; For each time window, determining a preliminary back azimuth of the time window according to the candidate back azimuth and the zero-lag cross-correlation coefficient; The back azimuth from the earthquake source to the monitoring point is determined based on the preliminary back azimuths of multiple time windows.

3. The method for determining the phase velocity of seismic waves according to claim 2, characterized in that: The step of determining a plurality of candidate back azimuth angles from the earthquake source to the monitoring point comprises: The back azimuth angles from 0° to 360° are divided into multiple parts at preset degree intervals to obtain multiple candidate back azimuth angles.

4. The method for determining the phase velocity of seismic waves according to claim 2, characterized in that: The step of determining the preliminary back azimuth of the time window according to the candidate back azimuth and the zero-lag cross-correlation coefficient comprises: The maximum zero-lag cross-correlation coefficient of the time window is determined, and the candidate back azimuth angle corresponding to the maximum zero-lag cross-correlation coefficient is determined as the preliminary back azimuth angle of the time window.

5. The method for determining the phase velocity of seismic waves according to claim 1, characterized in that: Determining the preliminary phase velocity of the time window according to the lateral acceleration and the vertical rotation component includes: Determine the maximum-magnitude lateral acceleration and the maximum-magnitude vertical rotation component in the time window; The preliminary phase velocity of the time window is determined according to the lateral acceleration with the largest amplitude and the vertical rotation component with the largest amplitude.

6. The method for determining the phase velocity of seismic waves according to claim 1, characterized in that: After determining the lateral acceleration of each time window and before determining the preliminary phase velocity of the time window, the method further includes: determining, for each time window, a zero-lag cross-correlation coefficient between the lateral acceleration and the vertical rotation component; The phase velocity of the time window that meets the preset condition is filtered out, and the preset condition is: the zero-lag mutual correlation coefficient of the time window is less than a preset threshold.

7. The method for determining the phase velocity of seismic waves according to claim 2 or 6, characterized in that: Determining a zero-lag correlation coefficient between the lateral acceleration and the vertical rotation component comprises: Determine the mean value of the lateral acceleration and the mean value of the vertical rotation component; A zero-lag correlation coefficient between the lateral acceleration and the vertical rotation component is determined based on the lateral acceleration, the vertical rotation component, the average value of the lateral acceleration, and the average value of the vertical rotation component.

8. A device for determining the phase velocity of seismic waves, characterized in that: The device comprises: A data acquisition unit is used to acquire seismic data collected by the monitoring point in multiple time windows, wherein the seismic data includes: a north-south translation component in the north-south direction, an east-west translation component in the east-west direction, and a vertical rotation component perpendicular to a plane formed by the north-south direction and the east-west direction; A back azimuth angle determination unit, used to determine the back azimuth angle from the earthquake source to the monitoring point; A lateral acceleration calculation unit, configured to determine the lateral acceleration according to the back azimuth angle and the north-south translational component and the east-west translational component for each time window; a phase velocity calculation unit, configured to determine, for each time window, a preliminary phase velocity of the time window according to the lateral acceleration and the vertical rotation component; The result output unit is used to determine the final phase velocity of the seismic wave according to the preliminary phase velocities of multiple time windows.

9. A computer device comprising a memory and a processor, wherein the memory stores computer-readable instructions, and when the computer-readable instructions are executed by the processor, the processor executes the steps of the method for determining the phase velocity of seismic waves as claimed in any one of claims 1 to 7.

10. A storage medium storing computer-readable instructions, which, when executed by one or more processors, causes the one or more processors to perform the steps of the method for determining the phase velocity of seismic waves as claimed in any one of claims 1 to 7.

Images