A single seismic location method based on staggered grid search

Abstract

The present invention discloses a single-station seismic location method based on staggered grid search, which relates to the technical field of seismic location. The present invention optimizes the problems of complex calculation and difficult determination of the focal depth of shallow earthquakes in the existing single-station seismic location method. By designing multiple sets of grids to intersect with each other in space, the spatial sampling rate of seismic location is improved, so that the location result has better stability and higher accuracy, and the reliability of the location result is proved by calculating the standard deviation of multiple sets of results.

Description

Technical Field

[0001] The present invention relates to the technical field of earthquake location, and in particular to a single-station earthquake location method based on staggered grid search. Background Art

[0002] The problems to be solved in earthquake location are to determine three basic parameters of the earthquake source location, origin time and magnitude based on the observation data of seismic stations. Accurate source parameters are of great significance for understanding issues such as the internal structure of the earth, earthquake generation mechanisms and evolution processes, and characteristics of volcanic activities, and also play an important role in earthquake prediction, earthquake prevention and disaster reduction, and earthquake rescue work. The accuracy of source parameters is related to factors such as the number and distribution of observation stations, the seismic phases used for location, arrival time picking, and velocity models. To improve the location accuracy, previous researchers have conducted a large number of studies and innovations on location methods from different perspectives. For the earthquake location of single-station observation data, the basic principle is to determine the azimuth and incident angle based on the P-wave first motion, determine the distance between the earthquake source and the station based on the arrival times of P-waves and S-waves, and obtain the earthquake source location by combining the azimuth, incident angle, and epicentral distance. The principle of this location method is easy to understand and the calculation is relatively simple. However, only a single-layer velocity model is used, which has a large difference from the real internal structure of the earth, and the influence of the earth's curvature on the location result is not considered, so the accuracy is relatively poor. There is a grid search method for automatic location using the three-component seismic waveform data of a single station. A radial grid is established with the station as the center, and different step sizes are set according to the epicentral distance and azimuth for searching. For each epicentral distance, calculate the cross-correlation function value of the predicted waveform of the earthquake event and the characteristic function of the recorded waveform. The epicentral distance corresponding to the maximum cross-correlation function value is the epicentral distance of the event. For each azimuth, synthesize the waveform components along the radial and tangential directions respectively, calculate the absolute integral difference between the radial and tangential directions of the P-wave, select the azimuth corresponding to the maximum value, and then judge the azimuth angle by the sign of the product integral of the synthesized radial component of the P-wave and the recorded vertical component. If the earthquake source depth is considered, the earthquake source depth condition is added during the search process of the epicentral distance, and the epicentral distance and earthquake source depth are determined simultaneously. Finally, the earthquake source location and origin time are determined according to the epicentral distance, earthquake source depth, and epicentral azimuth. However, this method is highly dependent on the data quality. If the data signal-to-noise ratio is low and the quality is poor, it may affect the location accuracy. The judgment of the azimuth angle is relatively complex. After rotating the horizontal components to the radial and tangential directions, it is necessary to select the correct azimuth angle from two results with a difference of 180°. Different results will be obtained by searching with different characteristic functions, and it is impossible to determine which characteristic function can obtain the optimal solution. The calculation of the earthquake source depth is ignored, and only the search along the depth is added in the area where deep-source earthquakes may exist, and the earthquake source depth is not searched for shallow-source earthquakes. Summary of the Invention

[0003] In view of the above deficiencies in the prior art, a single-station earthquake location method based on staggered grid search provided by the present invention solves the problems of complex calculation, difficult determination of the focal depth, and low location accuracy.

[0004] In order to achieve the above-mentioned invention object, the technical solution adopted by the present invention is as follows:

[0005] A single-station earthquake location method based on staggered grid search is provided, which includes the following steps:

[0006] S1. Obtain the original continuous seismic waveform data of a certain area;

[0007] S2. Pick up earthquake events from the original continuous seismic waveform data, calculate the corresponding time difference, back azimuth, and incident angle, and use them as true values;

[0008] S3. Obtain the epicentral distance and focal depth range corresponding to the earthquake event, and establish a three-dimensional grid centered on the observation station location based on the epicentral distance range and focal depth range;

[0009] S4. Construct a grid travel time table;

[0010] S5. Use a single set of three-dimensional grids to search and locate earthquake events, screen out grid points that meet the screening conditions, and set an objective function using the back azimuth, incident angle, and time difference in the grid points and true values. Take the grid point coordinates with the minimum objective function as the focal source location located by the single set of three-dimensional grids;

[0011] S6. Translate and combine a single set of three-dimensional grids N times to obtain a staggered grid model;

[0012] S7. Use the staggered grid model to search and locate earthquake events to obtain the final focal source location and error evaluation;

[0013] S8. Obtain the origin time of the earthquake event by forward modeling seismic rays according to the velocity model to complete the location of the earthquake; among them, the location result includes the final focal source location and the origin time.

[0014] Further, step S2 includes the following steps:

[0015] S2-1. Select earthquake events that meet the signal-to-noise ratio threshold in the original continuous seismic waveform data, perform arrival time picking and measurement to obtain the corresponding P-wave arrival time, S-wave arrival time, and amplitude data; among them, the amplitude data is the maximum amplitude in the first half cycle before the P-wave arrival of the three-component, including the vertical record u U , the east-west record u E and the north-south record u N ;

[0016] S2-2. Calculate the time difference, back azimuth, and incident angle.

[0017] Further, the vertical recording u in step S2-1 U is positive upward, the east-west recording u E is positive eastward, and the north-south recording u N is positive northward; when u U > 0, the seismic ray moves away from the epicenter. At this time, the east-west recording u E and the north-south recording u N take the opposite values.

[0018] Further, the calculation formula for the arrival time difference in step S2-2:

[0019] t S-P = t S - t P

[0020] where t S-P represents the arrival time difference, t S represents the arrival time of the S-wave, and t P represents the arrival time of the P-wave;

[0021] The calculation formula for the back azimuth:

[0022]

[0023] where baz represents the back azimuth, and tan -1 represents the arctangent function, u E represents the east-west recording, that is, the maximum amplitude of the P-wave first arrival in the east-west recording, and u N represents the north-south recording, that is, the maximum amplitude of the P-wave first arrival in the north-south recording;

[0024] The calculation formula for the incident angle:

[0025]

[0026] where i p represents the incident angle, represents the initial incident angle, u U represents the vertical recording, that is, the maximum amplitude of the P-wave first arrival in the vertical recording, |·| represents the absolute value, and sin represents the sine function. V S represents the propagation speed of the S-wave, and V P represents the propagation speed of the P-wave.

[0027] Further, the three-dimensional grid in step S3 includes a horizontal grid and a depth grid;

[0028] Step S3 includes the following steps:

[0029] S3-1. Determine the source depth range of the measured earthquake based on geological and geophysical data, and based on the formula:

[0030]

[0031] Obtain the epicentral distance Δ of the events monitored by the station; determine the epicentral distance range based on the arrival time differences of each earthquake event; where, t S-P represents the maximum arrival time difference of the locatable earthquakes monitored by the station, V P represents the propagation speed of the P wave, V S represents the propagation speed of the S wave;

[0032] S3-2. Based on the epicentral distance range and the location of the observation station, and set up a longitude-latitude grid centered on the location of the observation station to obtain a horizontal grid; where, the unit of the horizontal grid is degree;

[0033] S3-3. Based on the epicentral distance range and the source depth range, and set up an epicentral distance-depth grid with the location of the observation station as the vertex to obtain a depth grid; where, the unit of the depth grid is kilometer.

[0034] Furthermore, the specific process of step S4:

[0035] Calculate the distance Δ' i and back azimuth baz i ' of each horizontal grid point (lon i , lat i ') relative to the station, and use the taup function of the obspy library to calculate the depth grid point (depth i , Δ' i ) to obtain the corresponding theoretical travel time difference and incident angle That is, the travel time table.

[0036] Furthermore, step S5 includes the following steps:

[0037] S5-1. Take the back azimuth, incident angle, and arrival time difference in step S2 as the observed values (baz, i p , t S-P );

[0038] S5-2. Based on the screening conditions, and use each distance to screen the grid points of a single set of three-dimensional grids to obtain the screened grid points;

[0039] S5-3. Establish a three-dimensional coordinate system with the back azimuth as the x-axis, the incident angle as the y-axis, and the arrival time differences of the P wave and S wave as the z-axis;

[0040] S5-4. According to the formula:

[0041]

[0042] Obtain the distance Δd between the filtered grid points in the three-dimensional coordinate system and the observed values (baz, i p , t S-P ), and use it as the objective function; where baz i ”, represents the three-dimensional coordinates of the filtered grid points;

[0043] S5-5. Select the minimum objective function and obtain the corresponding grid point coordinates (lon i ', lat i ', depth i ) as the source location of the single set of three-dimensional grid positioning.

[0044] Furthermore, the screening conditions include a first screening condition and a second screening condition; where:

[0045] The first screening condition is:

[0046]

[0047]

[0048] where |·| represents the absolute value;

[0049] The second screening condition is:

[0050] Δs i ' ≤ Δ_max + 2

[0051] Δs i ' ≥ Δ_min - 2

[0052] |baz i ' - baz| ≤ 1

[0053] where Δs i ' represents the distance between the i-th horizontal grid point and (baz, i p , t S-P ), Δ_max represents the maximum epicentral distance of the depth grid points obtained by the first screening, and Δ_min represents the minimum epicentral distance of the depth grid points obtained by the first screening;

[0054] Step S5-4 includes the following steps:

[0055] S5-4-1. Screen the grid points of the depth grid based on the first screening condition;

[0056] S5-4-2. Screen the grid points of the horizontal grid based on the second screening condition;

[0057] The screened grid points include screened depth grid points and screened horizontal grid points.

[0058] Further, the specific process of step S6 is as follows:

[0059] The horizontal grid is translated along longitude and latitude with a step size of 1 / 4 of the small grid size, and the depth grid is translated along the horizontal direction and vertically downward with a step size of 1 / 4 of the small grid size, and is moved N times in total to form N + 1 mutually staggered grids; the N + 1 horizontal grids and N + 1 depth grids are combined in pairs to obtain corresponding grid cells; where the value of N is 15;

[0060] The staggered grid model includes (N + 1) 2 sets of grid cells.

[0061] Further, step S7 includes the following steps:

[0062] S7-1. Using (N + 1) 2 sets of grid cells and adopting the same method as steps S4 to S5 to search and locate the earthquake to obtain the corresponding earthquake source location;

[0063] S7-2. Calculate the average values of the longitudes, latitudes, and depths corresponding to all the earthquake source locations in step S7-1, and use them as the final earthquake location;

[0064] S7-3. According to the (N + 1) 2 sets of grid cells, obtain (N + 1) 2 earthquake source results (lone i , late i , depthe i ) and calculate the standard deviation as the error analysis of the location result. The formula is as follows:

[0065]

[0066] Among them, σ lon represents the longitude standard deviation, σ lat represents the latitude standard deviation, σ depth represents the depth standard deviation, lone i represents the longitude of the i-th earthquake source location result, ∑(·) represents the summation function, represents the average value of the longitudes among the (N + 1) 2 earthquake source location results, late i represents the latitude of the i-th earthquake source location result, represents the average value of the latitudes among the (N + 1) 2 earthquake source location results, depthe iRepresents the depth of the i-th seismic source location result, Indicates (N + 1) 2 The average value of the depths in the seismic source location results.

[0067] The beneficial effects of the present invention are as follows: This method solves the problems of complex calculation, difficult determination of the seismic source depth, and low positioning accuracy in the existing single-station seismic positioning method. By designing multiple sets of grids that intersect with each other in space, the spatial sampling rate of seismic positioning is improved, making the positioning result have better stability and higher accuracy. And the reliability of the positioning result is improved by calculating the standard deviation of multiple sets of results. Description of the Drawings

[0068] Figure 1 Is the flowchart of the method of the present invention;

[0069] Figure 2 Is the schematic diagram of the three-dimensional grid;

[0070] Figure 3 Is the schematic diagram of the three-dimensional grid intersection;

[0071] Figure 4 Is the waveform schematic diagram of a specific embodiment. Detailed Embodiment

[0072] The following describes the detailed embodiment of the present invention to facilitate those skilled in the art of this technology to understand the present invention. However, it should be clear that the present invention is not limited to the scope of the detailed embodiment. For those of ordinary skill in the art of this technology, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are within the scope of protection.

[0073] As Figure 1 shown, a single-station seismic positioning method based on cross-grid search includes the following steps:

[0074] S1. Obtain the original continuous seismic waveform data of a certain area;

[0075] S2. Pick up seismic events from the original continuous seismic waveform data, calculate the corresponding time difference, back azimuth, and incident angle, and use them as true values;

[0076] Step S2 includes the following steps:

[0077] S2-1. Select seismic events that meet the signal-to-noise ratio threshold in the original continuous seismic waveform data and perform arrival time picking and measurement to obtain the corresponding P-wave arrival time, S-wave arrival time, and amplitude data; among them, the amplitude data is the maximum amplitude in the first half cycle before the P-wave arrival of the three-component, including the vertical recording u U 、east-west recording u Eand the north-south record u N ;

[0078] the vertical record u in step S2-1 U Upward is positive, the east-west record u E Eastward is positive, the north-south record u N Northward is positive; when u U > 0, the seismic ray moves away from the epicenter. At this time, the east-west record u E and the north-south record u N take the opposite value.

[0079] S2-2. Calculate the time difference, back azimuth, and incident angle.

[0080] The calculation formula for the time difference in step S2-2:

[0081] t S-P = t S - t P

[0082] where t S-P represents the time difference, t S represents the arrival time of the S wave, t P represents the arrival time of the P wave;

[0083] The calculation formula for the back azimuth:

[0084]

[0085] where baz represents the back azimuth, tan -1 represents the arctangent function, u E represents the east-west record, that is, the maximum amplitude of the P-wave first arrival in the east-west record, u N represents the north-south record, that is, the maximum amplitude of the P-wave first arrival in the north-south record;

[0086] The calculation formula for the incident angle:

[0087]

[0088] where i p represents the incident angle, represents the initial incident angle, u U represents the vertical record, that is, the maximum amplitude of the P-wave first arrival in the vertical record, |·| represents the absolute value, sin represents the sine function, V S represents the propagation speed of the S wave, V P represents the propagation speed of the P wave.

[0089] Step S2 ensures the quality of data preprocessing by manually measuring the amplitude and arrival time. To a certain extent, it can filter out event waveforms with low signal-to-noise ratio and poor quality, and only measure and locate event waveforms with high signal-to-noise ratio and direct P-wave phases. Moreover, it directly uses the P-wave first arrival amplitude measured from the three-component seismic waveforms to calculate the back azimuth and incident angle of the event for searching, which is highly reliable, simple to calculate, and not prone to errors.

[0090] S3. Obtain the epicentral distance and focal depth range corresponding to the seismic event, and establish a three-dimensional grid centered on the location of the observation station based on the epicentral distance range and focal depth range;

[0091] The three-dimensional grid in step S3 includes a horizontal grid and a depth grid;

[0092] As Figure 2 shown, step S3 includes the following steps:

[0093] S3-1. Determine the focal depth range of the measured earthquake according to geological and geophysical data, and based on the formula:

[0094]

[0095] Obtain the epicentral distance Δ of the event monitored by the station; determine the epicentral distance range based on the arrival time differences of each seismic event; where, t S-P represents the maximum arrival time difference of the locatable earthquake monitored by this station, V P represents the P-wave propagation velocity, and V S represents the S-wave propagation velocity;

[0096] S3-2. Based on the epicentral distance range and the location of the observation station, and with the location of the observation station as the center, set the longitude-latitude grid to obtain the horizontal grid; where, the unit of the horizontal grid is degree;

[0097] S3-3. Based on the epicentral distance range and the focal depth range, and with the location of the observation station as the apex, set the epicentral distance-depth grid to obtain the depth grid; where, the unit of the depth grid is kilometer.

[0098] In Figure 2 , the right grid is a quarter of the horizontal grid, and the left grid is a schematic diagram of the depth grid. The triangle represents the station, and the pentagram represents the earthquake.

[0099] The scale of the horizontal grid is set according to the locatable earthquake epicentral distance range estimated from the arrival time difference and virtual wave velocity, and the depth grid is set according to the estimated epicentral distance and the focal depth range of the earthquakes in the region. The internal small grid size can also be set as required, greatly shortening the time for generating the grid model.

[0100] S4. Construct a grid travel time table;

[0101] Specific process of step S4:

[0102] Calculate the distance Δ' i ' and back azimuth baz i ' of each horizontal grid point (lon i , lat i ') relative to the station, and use the taup function of the obspy library to calculate the depth grid point (depth i ', Δ' i ), obtaining the corresponding theoretical travel time difference and incident angle That is, the travel time table.

[0103] S5. Use a single set of three-dimensional grids to search and locate seismic events, screen out the grid points that meet the screening conditions, and set the objective function using the back azimuth, incident angle, and arrival time difference in the grid points and the true values. Take the grid point coordinates with the minimum objective function as the source location of the single set of three-dimensional grid positioning;

[0104] Step S5 includes the following steps:

[0105] S5-1. Use the back azimuth, incident angle, and arrival time difference in step S2 as the observed values (baz,i p , t S-P );

[0106] S5-2. Based on the screening conditions, and use each distance to screen the grid points of the single set of three-dimensional grids, obtaining the screened grid points;

[0107] S5-3. Establish a three-dimensional coordinate system with the back azimuth as the x-axis, the incident angle as the y-axis, and the arrival time difference of P-wave and S-wave as the z-axis;

[0108] S5-4. According to the formula:

[0109]

[0110] Obtain the distance Δd between the screened grid points in the three-dimensional coordinate system and the observed values (baz,i p , t S-P ), and use it as the objective function; where baz i ” represents the three-dimensional coordinates of the screened grid points;

[0111] S5-5. Select the minimum objective function, and obtain the corresponding grid point coordinates (lon i ', lat i ', depth i') as the source location for single-set 3D grid positioning.

[0112] The screening conditions include a first screening condition and a second screening condition; where:

[0113] The first screening condition is:

[0114]

[0115] where, |·| represents the absolute value;

[0116] The second screening condition is:

[0117] Δs i ' ≤ Δ_max + 2

[0118] Δs i ' ≥ Δ_min - 2

[0119] |baz i ' - baz| ≤ 1

[0120] where, Δs i ' represents the distance between the i-th horizontal grid point and (baz, i p , t S-P ), Δ_max represents the maximum epicentral distance of the depth grid points obtained from the first screening, and Δ_min represents the minimum epicentral distance of the depth grid points obtained from the first screening;

[0121] Step S5-4 includes the following steps:

[0122] S5-4-1. Screen the grid points of the depth grid based on the first screening condition;

[0123] S5-4-2. Screen the grid points of the horizontal grid based on the second screening condition;

[0124] The screened grid points include the screened depth grid points and the screened horizontal grid points.

[0125] S6. Perform N translations and combinations on the single-set 3D grid to obtain an interleaved grid model;

[0126] As Figure 3 shown, the specific process of step S6 is:

[0127] The horizontal grid is translated along longitude and latitude in steps of 1 / 4 of the small grid size, and the depth grid is translated along the horizontal direction and vertically downward in steps of 1 / 4 of the small grid size. This is done N times to form N + 1 interleaved grids. The N + 1 horizontal grids and N + 1 depth grids are combined pairwise to obtain the corresponding grid cells. Here, N takes the value of 15. This step improves the accuracy and search stability of single - station earthquake location, and the spatial sampling rate reaches 16 times that of a single grid.

[0128] The staggered - grid model includes (N + 1) 2 sets of grid cells.

[0129] S7. Use the staggered - grid model to search for and locate earthquake events to obtain the final source location and error evaluation;

[0130] Step S7 includes the following steps:

[0131] S7 - 1. Use (N + 1) 2 sets of grid cells and use the same method as in steps S4 to S5 to search for and locate the earthquake to obtain the corresponding source location;

[0132] S7 - 2. Calculate the average values of the longitudes, latitudes, and depths corresponding to all the source locations in step S7 - 1 and use them as the final earthquake location;

[0133] S7 - 3. Calculate the standard deviation based on the (N + 1) 2 sets of grid cells to obtain the (N + 1) 2 source results (lone i , late i , depthe i ) and use it as the error analysis of the location result. The formula is as follows:

[0134]

[0135] where σ lon represents the longitude standard deviation, σ lat represents the latitude standard deviation, σ depth represents the depth standard deviation, lone i represents the longitude of the i - th source location result, ∑(·) represents the summation function, represents the average value of the longitudes in the (N + 1) 2 source location results, late i represents the latitude of the i - th source location result, represents the average value of the latitudes in the (N + 1) 2 source location results, depthe i represents the depth of the i - th source location result, denotes (N + 1) 2 the average value of the depths in the source location results

[0136] In the location of the next earthquake event, referring to the error analysis of the previous earthquake event, the location result of the next earthquake event is adjusted.

[0137] S8. Forward model seismic rays according to the velocity model to obtain the origin time of the earthquake event, and complete the location of the earthquake; among them, the location result includes the final source location and the origin time.

[0138] The present invention also includes the reliability evaluation of the location result. The standard deviation can reflect the dispersion degree of a set of data. Therefore, the present invention calculates the standard deviation of the epicenter positions matched in different grid combinations for an event to evaluate the stability of the location result. The smaller the standard deviation, the smaller the influence of the grid model difference on the result, and the more stable the single-station earthquake location result.

[0139] In an embodiment of the present invention, a near-earthquake event recorded by a seismic station in a certain area at UTC 2019-11-01 22:56:18.964 is selected. As Figure 4 shown, the arrival time of the P wave is 22:56:24.594, and the arrival time of the S wave is 22:56:29.454; the vertical component record u U , the east-west component record u E and the north-south component record u N are -231 counts, 163 counts and -1200 counts respectively; calculate the time difference t S-P of this event to be 4.86 s, the back azimuth baz to be 305.21°, and the incident angle i p to be 11.99°; search in the travel time table of each set of grid models with the measurement parameters as the standard, and take the average value of the 256 source positions with the minimum objective function to obtain the final result (lon i ', lat i ', depth i ') to be (-59.286°, -62.114°, 20.36 km), the corresponding back azimuth to be 304.44°, and the incident angle to be 12.00°; waveform forward modeling gives the P wave travel time of 5.133 s, the origin time of 22:56:19.460, the longitude standard deviation of 0.01°, the latitude standard deviation of 0.00°, and the depth standard deviation of 0.31 km. In Figure 4 , a is the schematic diagram of the three-component waveform of the locatable near-earthquake event, b is the schematic diagram of the maximum amplitude waveform corresponding to the three-component P wave first arrival of the locatable event, c is the three-component original waveform record of the non-locatable earthquake, and d is the three-component record of the P wave first motion amplification of the earthquake event in c.

[0140] In summary, the present invention solves the problems of complex calculation, difficult determination of focal depth and low positioning accuracy in the existing single-station earthquake positioning method. By designing multiple sets of grids to intersect with each other in space, the spatial sampling rate of earthquake positioning is improved, so that the positioning result has good stability and high accuracy. And the reliability of the positioning result is proved by calculating the standard deviation of multiple sets of results, and good constraint effects on the focal depths of shallow, intermediate and deep earthquakes are achieved. By calculating the time differences between P-waves and S-waves at grid points through the TauP function, this function directly solves the equation in the spherical coordinate system without flattening transformation, taking into account the influence of the earth's curvature and having strong adaptability. Especially in the Antarctic and Arctic regions where the station distribution is sparse, this method can weaken the influence of curvature on the positioning accuracy to the greatest extent. By combining the staggered grids in pairs to form 256 sets of three-dimensional grids, 256 positioning results are correspondingly obtained, and the standard deviations of the longitude, latitude and depth of the seismic event are calculated as error evaluation to verify the stability of the positioning result.

Claims

1. A single-seismic positioning method based on staggered grid search, characterized in that: It includes the following steps: S1. Obtain the original continuous seismic waveform data of a certain area; S2. Pick up seismic events from the original continuous seismic waveform data, calculate the corresponding time difference, back azimuth, and incident angle, and use them as true values; S3. Obtain the epicentral distance and focal depth range corresponding to the seismic events in this area, and establish a single set of three-dimensional grids centered on the location of the observation station based on the epicentral distance range and focal depth range; S4. Construct a grid travel time table; S5. Use the single set of three-dimensional grids to search and locate seismic events, screen out the grid points that meet the screening conditions, and set the objective function using the back azimuth, incident angle, and time difference in the grid points and true values. Take the grid point coordinates with the minimum objective function as the focal source position located by the single set of three-dimensional grids; S6. Translate and combine the single set of three-dimensional grids N times to obtain a staggered grid model; S7. Use the staggered grid model to search and locate seismic events to obtain the final focal source position and error evaluation; S8. Forward model seismic rays according to the velocity model to obtain the origin time of the seismic events, and complete the location of the earthquake; among them, the location result includes the final focal source position and the origin time; Step S5 includes the following steps: S5-1. Take the back azimuth, incident angle, and arrival time difference in step S2 as the observed values (baz, i p , t S-P ); S5-2. Based on the screening conditions and using each distance, screen the grid points of a single set of three-dimensional grids to obtain the screened grid points; S5-3. Establish a three-dimensional coordinate system with the back azimuth as the x-axis, the incident angle as the y-axis, and the time difference between P-wave and S-wave as the z-axis; S5-4. According to the formula: Obtain the distance Δd between the filtered grid points in the three-dimensional coordinate system and the observed values (baz, i p , t S-P ), and use it as the objective function; where baz” i 、 represents the three-dimensional coordinates of the filtered grid points; S5-5. Select the minimum objective function and obtain the corresponding grid point coordinates (lon’ i , lat’ i , depth’ i ) as the source location for single-set three-dimensional grid positioning; The screening conditions include the first screening condition and the second screening condition; where: The first screening condition is: where, |·| represents the absolute value; The second screening condition is: Δs’ i ≤Δ_max + 2 Δs’ i ≥Δ_min - 2 |baz’ i -baz| ≤ 1 where, Δs’ i represents the distance between the i-th horizontal grid point and (baz,i p ,t S-P ), Δ_max represents the maximum epicentral distance of the depth grid points obtained by the first screening, and Δ_min represents the minimum epicentral distance of the depth grid points obtained by the first screening; Step S5-4 includes the following steps: S5-4-1. Based on the first screening condition, screen the grid points of the depth grid; S5-4-2. Based on the second screening condition, screen the grid points of the horizontal grid; The screened grid points include the screened depth grid points and the screened horizontal grid points.

2. The single seismic location method based on staggered grid search according to claim 1, characterized in that: Step S2 includes the following steps: S2-1. Select seismic events that meet the signal-to-noise ratio threshold from the original continuous seismic waveform data, and perform arrival time picking and measurement to obtain the corresponding P-wave arrival time, S-wave arrival time, and amplitude data; among them, the amplitude data is the maximum amplitude in the first half cycle before the P-wave first arrival of the three-component, including the vertical component record u U , the east-west component record u E and the north-south component record u N ; S2-2. Calculate the time difference, back azimuth, and incident angle.

3. The single seismic location method based on staggered grid search according to claim 2, wherein: The vertical recording u in step S2-1 U is positive upward, and the east-west recording u E is positive eastward, and the north-south recording u N is positive northward; when u U > 0, the seismic ray moves away from the epicenter. At this time, the east-west recording u E and the north-south recording u N take the opposite value.

4. The single-seismic event location method based on staggered grid search according to claim 2, wherein: The calculation formula for the time difference in step S2-2: t S-P = t S -t P where t S-P represents the time difference, t S represents the S-wave arrival time, and t P represents the P-wave arrival time; The calculation formula for the back azimuth: where baz represents the back azimuth, and tan -1 represents the arctangent function, and u E represents the east-west record, that is, the maximum amplitude of the P-wave first arrival in the east-west record, and u N represents the north-south record, that is, the maximum amplitude of the P-wave first arrival in the north-south record; The calculation formula for the incident angle: where i p represents the incident angle, represents the apparent incident angle, u U represents the vertical recording, that is, the maximum amplitude of the P-wave first arrival in the vertical recording, |·| represents the absolute value, and sin represents the sine function, V S represents the propagation velocity of the S-wave, V P represents the propagation velocity of the P-wave.

5. The single-seismic event location method based on staggered grid search according to claim 2, wherein: The three-dimensional grids in step S3 include horizontal grids and depth grids; Step S3 includes the following steps: S3-1. Determine the focal depth range of the measured earthquake according to geological and geophysical data, and based on the formula: Obtain the epicentral distance Δ of the events monitored by the station; determine the epicentral distance range based on the arrival-time differences of each earthquake event; where, t S-P represents the maximum arrival-time difference of the locatable earthquakes monitored by the station, V P represents the propagation speed of the P wave, V S represents the propagation speed of the S wave; S3-2. Based on the epicentral distance range and the location of the observation station, and set a longitude-latitude grid centered on the location of the observation station to obtain a horizontal grid; where, the unit of the horizontal grid is degree; S3-3. Based on the epicentral distance range and the focal depth range, and set an epicentral distance-depth grid with the location of the observation station as the vertex to obtain a depth grid; where, the unit of the depth grid is kilometer.

6. The single-seismic event location method based on staggered grid search according to claim 5, wherein: The specific process of step S4: Calculate the distance Δ' i and back azimuth baz' i of each horizontal grid point (lon’ i , lat’ i ) relative to the station, and use the taup function of the obspy library to calculate the depth grid point (depth’ i , Δ' i ) to obtain the corresponding theoretical travel time difference and the incident angle i.e., the travel time table.

7. The single-seismic event location method based on staggered grid search according to claim 5, wherein: The specific process of step S6 is: The horizontal grid is translated along the longitude and latitude in steps of 1 / 4 of the small grid size, and the depth grid is translated along the horizontal direction and vertically downward in steps of 1 / 4 of the small grid size. This is done N times to form N + 1 mutually staggered grids. The N + 1 horizontal grids and N + 1 depth grids are combined pairwise to obtain the corresponding grid cells. Here, the value of N is 15; The staggered grid model includes (N + 1) 2 sets of grid cells.

8. The single seismic positioning method based on staggered grid search according to claim 7, wherein: Step S7 includes the following steps: S7-1. Use (N + 1) 2 sets of grid cells and uses the same method as in steps S4 to S5 to search for and locate the earthquake, and obtains the corresponding epicenter location; S7-2. Calculate the average values of the longitudes, latitudes, and depths corresponding to all the seismic source positions in step S7-1, and use them as the final earthquake position; S7-3. According to (N+1) 2 sets of grid cells to obtain (N+1) 2 source results (lone i , late i , depth i ) Calculate the standard deviation and use it as the error analysis of the positioning result. The formula is as follows: Among them, σ lon represents the longitude standard deviation, σ lat represents the latitude standard deviation, σ depth represents the depth standard deviation, lone i represents the longitude of the i-th earthquake source location result, and ∑(·) represents the summation function, represents (N + 1) 2 the average value of longitudes in i the (N + 1) earthquake source location results, late represents (N + 1) 2 the average value of latitudes in i the (N + 1) earthquake source location results, depthe represents (N + 1) 2 the average value of depths in the (N + 1) earthquake source location results.

Images