A method and apparatus for earthquake association and earthquake location

Abstract

The embodiment of the application provides a method and device for earthquake association and earthquake positioning, comprising: generating a plurality of initial sampling points in a preset search area of a seismic source parameter; calculating a target value of each initial sampling point according to a preset target function; wherein the target function is constructed according to the number of seismic phases matched by the seismic source parameter and the positioning reliability of the seismic source parameter; determining a local area in which the initial sampling points with a target value greater than a preset target threshold are located according to the distribution of the target values of the initial sampling points; generating a plurality of local sampling points in the local area; calculating the target value of each local sampling point; finding an optimal sampling point in which the target function reaches a maximum target value based on the target values of the local sampling points; and determining the seismic source parameter corresponding to the optimal sampling point and the matched seismic phase. The application can quickly and efficiently obtain the earthquake association and positioning result, improve the estimation accuracy of the seismic source parameter, and enhance the accuracy and reliability of the earthquake association result.

Description

Technical Field

[0001] This application relates to the field of earthquake observation technology, and in particular to a seismic correlation and earthquake location method and apparatus. Background Technology

[0002] Earthquake catalogs are crucial for characterizing regional fault structures, studying the activity characteristics of large earthquakes, and conducting seismic statistical parameter studies. After phase picking, seismic correlation is a key step in constructing the seismic catalog. Its purpose is to identify phase groups belonging to the same seismic event from a large number of phase picking results, effectively distinguish them from other seismic events, and eliminate erroneous or unreliable phases, thereby improving the data quality for subsequent earthquake location and analysis. Relevant seismic correlation methods include grid search-based methods, deep learning neural network-based methods, and clustering-based methods. These methods suffer from drawbacks such as long processing times, poor correlation performance under high noise levels, and complex and diverse dependent parameters. Summary of the Invention

[0003] In view of this, the purpose of this application is to provide a seismic correlation and earthquake location method and apparatus.

[0004] To achieve the above objectives, embodiments of this application provide a seismic correlation and earthquake location method, including:

[0005] Multiple initial sampling points are generated within the preset source parameter search area;

[0006] According to the preset objective function, the target value of each initial sampling point is calculated; wherein, the objective function is constructed based on the number of seismic phases matched with the source parameters and the positioning reliability of the source parameters;

[0007] Based on the distribution of target values ​​at each initial sampling point, determine the local region where the initial sampling points with target values ​​greater than the preset target threshold are located;

[0008] Multiple local sampling points are generated within the local area;

[0009] Calculate the target value for each local sampling point. Based on the target values ​​of each local sampling point, find the optimal sampling point where the objective function reaches its maximum target value, and determine the source parameters and matching phase corresponding to the optimal sampling point.

[0010] Optionally, the method further includes:

[0011] Using the nearest station where P-waves were observed as the initial reference point, the predetermined spatial search range is expanded horizontally; the temporal search range is determined based on the average velocity of P-waves and the preset horizontal search range; and the depth search range of the seismic source is determined based on geological structural features and the distribution of seismic activity.

[0012] The search area for the source parameters is determined based on the spatial search range, temporal search range, and depth search range.

[0013] Optionally, multiple initial sampling points are generated within a preset source parameter search area, including:

[0014] The number of initial sampling points is determined based on the search area for the earthquake source parameters;

[0015] The Sobol sampling method is used to randomly generate initial sampling points with this number and source parameter dimensions.

[0016] Optionally, the target value for each initial sampling point is calculated according to a preset objective function, including:

[0017] Using each initial sampling point as the seismic source, determine the travel time residual between the seismic source and any station;

[0018] If the travel time residual is less than the preset residual threshold, the initial sampling point is determined to match the station.

[0019] Statistically analyze the matching results between all initial sampling points and each station;

[0020] Based on the matching results, the target value of the objective function is calculated.

[0021] Optionally, the objective function is:

[0022] f(x1,x2,x3,t)=W N ·S N (x1,x2,x3,t)+W L ·S L (x1,x2,x3,t)

[0023] Among them, S N (x1,x2,x3,t) represents the number of seismic phases matched under the source parameters (x1,x2,x3,t), S L (x1,x2,x3,t) represents the location reliability of the source parameters (x1,x2,x3,t), W N W L x1 is the weight, x2 is the latitude, x3 is the depth, and t is the time of the earthquake.

[0024] Optionally, the method for calculating the number of seismic phases is as follows:

[0025]

[0026] Where, N P N represents the number of P-waves matched at the sampling points. SW represents the number of S-waves matched to the sampling points. P W S For the weight, N station This refers to the total number of stations.

[0027] Optionally, the method for calculating the positioning reliability is as follows:

[0028]

[0029] Where n is the number of stations associated with the sampling point; Let be the travel time residual of the P-wave corresponding to the i-th station. W represents the travel time residual of the S-wave corresponding to the i-th station; i The weight of the i-th station is calculated based on its distance from the sampling point; MaxTolerance is the preset maximum tolerance window.

[0030] Optionally, based on the distribution of target values ​​at each initial sampling point, the local region where the target value is greater than a preset target threshold is determined, including:

[0031] The longitude, latitude, depth, and time of the initial sampling points where the target value is greater than the target threshold are within the quantile range of the longitude, latitude, depth, and time dimensions, respectively.

[0032] The local region is determined based on the quantile range.

[0033] Optionally, finding the optimal sampling point where the objective function reaches its maximum objective value based on the target values ​​of each local sampling point includes:

[0034] Based on the target value of each local sampling point, the sub-region where the local sampling point with the target value is greater than the target threshold is determined; within the sub-region, multiple sub-region sampling points are generated; the target value of each sub-region sampling point is calculated; the search area is gradually narrowed according to this process until the preset stopping condition is reached, and the maximum target value and the corresponding optimal sampling point are determined.

[0035] This application embodiment also provides a seismic correlation and earthquake location device, including:

[0036] The initial sampling point generation module is used to generate multiple initial sampling points within a preset source parameter search area;

[0037] The calculation module is used to calculate the target value of each initial sampling point according to a preset objective function; wherein, the objective function is constructed based on the number of seismic phases matched with the source parameters and the positioning reliability of the source parameters;

[0038] The local region determination module is used to determine the local region where the target value of the initial sampling point is greater than the preset target threshold, based on the distribution of the target values ​​of each initial sampling point.

[0039] A local sampling point generation module is used to generate multiple local sampling points within the local area;

[0040] The optimal solution determination module is used to calculate the target value of each local sampling point, and based on the target value of each local sampling point, find the optimal sampling point where the objective function reaches the maximum target value, and determine the source parameters and matching phase corresponding to the optimal sampling point.

[0041] As can be seen from the above description, the seismic correlation and earthquake location method and apparatus provided in this application generate multiple initial sampling points within the source parameter search area, calculate the target value of each initial sampling point according to the objective function, determine the local area where the target value of each initial sampling point is greater than a preset target threshold based on the distribution of the target values ​​of each initial sampling point, generate multiple local sampling points within the local area, calculate the target value of each local sampling point, and find the optimal sampling point where the objective function reaches its maximum target value based on the target values ​​of each local sampling point, and determine the source parameters and matching seismic phase corresponding to the optimal sampling point. The method of this application can quickly and efficiently obtain seismic correlation and location results, improve the estimation accuracy of source parameters, and enhance the accuracy and reliability of seismic correlation results. Attached Figure Description

[0042] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0043] Figure 1 This is a schematic diagram of the method flow of an embodiment of this application;

[0044] Figure 2 This is a block diagram of the device structure according to an embodiment of this application;

[0045] Figure 3 This is a block diagram of the electronic device structure according to an embodiment of this application. Detailed Implementation

[0046] To make the objectives, technical solutions, and advantages of this disclosure clearer, the following detailed description is provided in conjunction with specific embodiments and the accompanying drawings.

[0047] It should be noted that, unless otherwise defined, the technical or scientific terms used in the embodiments of this application should have the ordinary meaning understood by one of ordinary skill in the art to which this disclosure pertains. The terms "first," "second," and similar terms used in the embodiments of this application do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed after the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are only used to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0048] Among related technologies, grid-based seismic correlation methods divide the target search area into multiple regular grids and determine whether the correlation conditions are met by evaluating the number of associated seismic phases in each grid. This method is slow, time-consuming, and dependent on the grid size and the allowable error range related to the grid size. Neural network-based seismic correlation methods use millions of P- and S-wave arrival time series generated based on one-dimensional velocity models as training data to construct a convolutional neural network for seismic correlation tasks. Alternatively, graph neural networks can be used to model the station set and the source space separately, learning the spatial-temporal relationship between seismic phases through a joint graph structure to achieve automatic matching of sources and phases. However, this method performs poorly under high noise levels, and the correlation effect is related to the quantity and quality of training samples. Cluster-based seismic correlation methods use multivariate Gaussian distributions to model the seismic phase set corresponding to earthquake events and use the expectation-maximization algorithm to jointly estimate seismic correlation and source parameters such as epicenter location, time of origin, and magnitude under the maximum likelihood criterion. This method is also affected by noise levels and performs poorly under high noise levels.

[0049] In view of this, embodiments of this application provide a seismic correlation and earthquake location method. Based on the number of seismic correlations and the reliability of location, an objective function is constructed. Within the search area of ​​source parameters, the search range is gradually narrowed to find the optimal solution that can simultaneously maximize the number of seismic correlations and take into account the location accuracy. This method can quickly and efficiently obtain seismic correlation and location results, improve the estimation accuracy of source parameters, and enhance the accuracy and reliability of seismic correlation results.

[0050] The technical solution of this application will be further described in detail below through specific embodiments.

[0051] like Figure 1As shown, this application provides a seismic correlation and seismic location method, including:

[0052] S101: Generate multiple initial sampling points within the preset source parameter search area;

[0053] In this embodiment, within the defined source parameter search area, multiple approximately uniformly distributed initial sampling points are quickly generated using a preset sampling method. That is, several sampling points are initialized globally as an initial solution space, and the region where the source is located is explored based on the distribution of the initial sampling points in the initial solution space within the source parameter search area.

[0054] In some implementations, methods for dividing the source parameter search area include:

[0055] Using the nearest station where P-waves were observed as the initial reference point, the predetermined spatial search range is expanded horizontally. The temporal search range is determined based on the average velocity of P-waves and the preset horizontal search range. The depth search range of the seismic source is determined based on geological structural features and the distribution of seismic activity.

[0056] The search area for earthquake source parameters is determined based on the spatial search range, temporal search range, and depth search range.

[0057] In this embodiment, considering that P-waves are relatively stable in seismic phase identification and that the arrival phase of P-waves is closer to the time of earthquake occurrence than that of S-waves, the source parameter search area is divided based on P-waves. During the division, the location of the nearest station recording the arrival time of P-waves is used as the initial reference point. Spatially, a certain range is extended horizontally along the initial reference point to form a spatial search range. Temporally, the propagation time of P-waves is estimated by combining the average velocity of P-waves within the region and the certain search range along the horizontal direction. A temporal search range is formed based on the time before and after the arrival time of P-waves. In terms of depth, the depth search range of the source is determined based on the tectonic features of the study area and the depth distribution of existing seismic activity. Thus, the source parameter search area is determined based on the spatial search range, temporal search range, and depth search range.

[0058] Optionally, under one experimental condition, with the initial reference point as the center, a rectangular spatial search range is formed by extending 0.6° in each of the east, west, north, and south directions in the spatial dimension. In the time dimension, the time search range is from 10 seconds before to 1 second after the arrival of the P wave. In the depth dimension, the depth search range is 0–30 km, thus constructing the source parameter search area.

[0059] After delineating the source parameter search area, several approximately uniformly distributed initial sampling points are generated within this area. The method includes: determining the number of initial sampling points based on the source parameter search area, and then randomly generating initial sampling points with the specified number and source parameter dimensions using the Sobol sampling method. Specifically, the number of initial sampling points suitable for the area search is first determined based on the size of the source parameter search area, considering factors such as computational efficiency and search effectiveness. Then, the Sobol sampling method is used to randomly generate the specified number of approximately uniformly distributed initial sampling points within the source parameter search area. These initial sampling points are nearly uniformly distributed points in the four-dimensional space of the source parameter solutions for seismic correlation and earthquake location, possessing four dimensions: longitude, latitude, depth, and time. The random and approximately uniform sampling points generated using the Sobol sampling method not only avoid the shortcomings of area exploration but also allow for rapid distribution within the source parameter search area, effectively improving search speed and efficiency.

[0060] S102: Calculate the target value for each initial sampling point according to the preset objective function; wherein, the objective function is constructed based on the number of seismic phases matched with the source parameters and the positioning reliability of the source parameters;

[0061] In this embodiment, the seismic association and location of an earthquake are treated as a joint optimization problem. An objective function to be optimized is constructed based on the number of seismic phases matched with the source parameters and the location reliability of the source parameters. By finding the optimal sampling point within the source parameter search area, the objective function can reach the maximum target value under the source parameters corresponding to the optimal sampling point, thus obtaining the optimal source parameters and associated seismic phases.

[0062] In some approaches, the objective function is constructed as follows:

[0063] f(x1,x2,x3,t)=W N ·S N (x1,x2,x3,t)+W L ·S L (x1,x2,x3,t) (1)

[0064] Among them, S N (x1,x2,x3,t) represents the number of seismic phases matched under the source parameters (x1,x2,x3,t), S L (x1,x2,x3,t) represents the location reliability of the source parameters (x1,x2,x3,t), W N W L x1 is the weight, x2 is the latitude, x3 is the depth, and t is the time of the earthquake.

[0065] The method for calculating the number of seismic phases matching the source parameters (x1,x2,x3,t) is as follows:

[0066]

[0067] In the formula, N P N represents the number of P-waves that can be associated with this sampling point. S W represents the number of S-waves that can be associated with this sampling point. P and W S The weights are denoted as N, and the sum of the two weights is 1. station The total number of stations set.

[0068] Considering that the difference between the theoretical and observed travel times of seismic phases (i.e., travel time residuals) imposes constraints on the source parameters, and that the distance between stations affects the reliability of the travel time residuals, a location reliability term based on epicentral distance weighting that satisfies the travel time residual constraints is constructed, expressed as:

[0069]

[0070] In the formula, n is the number of stations associated with the sampling point; Let be the travel time residual of the P-wave corresponding to the i-th station. W represents the travel time residual of the S-wave corresponding to the i-th station; P and W S These are the weights of the travel time residuals for the P-wave and S-wave, respectively, and the sum of the two weights is 1; W i The weight of the i-th station is calculated based on its distance from the sampling point. Once the sampling point is determined, the distance between the station and the sampling point is determined according to the longitude and latitude of the sampling point and the station's location. The ratio of this distance to the defined maximum distance is calculated to obtain the weight of the station; MaxTolerance is the preset maximum tolerance window; S L ∈[0,1].

[0071] For the positioning reliability term shown in equation (3), distant stations are affected by factors such as low signal-to-noise ratio, large phase picking error, and large systematic error, resulting in lower reliability of their travel time residuals. Nearby stations are less affected by errors and other factors, leading to higher reliability of their travel time residuals. Therefore, the weight W corresponding to positioning reliability in the objective function... L The weight can be adjusted based on the distance between stations; for distant stations, the weight can be reduced, and for nearby stations, the weight can be increased.

[0072] In some implementations, the target value for each initial sampling point is calculated according to the constructed objective function, including:

[0073] Using each initial sampling point as the seismic source, determine the travel time residual between the seismic source and any station;

[0074] If the travel time residual is less than the preset residual threshold, the initial sampling point is determined to match the station.

[0075] Statistically analyze the matching results between all initial sampling points and each station;

[0076] Based on the matching results, calculate the target value of the objective function.

[0077] In this embodiment, for each initial sampling point, the travel time residual between the source and each station within the source parameter search area is determined, with the initial sampling point as the source. If the travel time residual is less than the residual threshold, the initial sampling point is considered to match the station, and the earthquake generated by the initial sampling point is associated with the station. The association results between all initial sampling points and all stations are determined in this way. Then, for the associated initial sampling points and stations, the matching results between all initial sampling points and each station are counted. The specific matching results include, but are not limited to: counting the number of P-waves and S-waves that match the source parameters of each initial sampling point with each station, counting the number of stations that match the source parameters of each initial sampling point, and recording the travel time residuals of the P-waves and S-waves matched by each station. Based on the matching results, the number of seismic phases that match the source parameters of each initial sampling point is calculated according to formula (2), the positioning reliability of the source parameters of each initial sampling point is calculated according to formula (3), and the target value of each initial sampling point is calculated according to the objective function shown in formula (1) based on the calculation results.

[0078] In some implementations, using each initial sampling point as the seismic source, the travel time residual between the seismic source and any station is determined, including:

[0079] Based on the source location of the initial sampling point and the epicentral distance of the station, the corresponding theoretical travel time is obtained from the preset phase travel time table;

[0080] Based on the theoretical travel time and the time of origin of the initial sampling point, the theoretical arrival time from the initial sampling point to the station is determined;

[0081] The travel time residual is calculated based on the theoretical arrival time and the observed arrival time of the seismic waves recorded by the station.

[0082] In this embodiment, for each initial sampling point, based on its longitude, latitude, depth, and epicentral distance from any station, a pre-constructed seismic phase travel time table is consulted to obtain the corresponding theoretical travel time, which is the time required for the seismic wave generated from the initial sampling point as the source to propagate to that station. Based on the seismic origin time of the initial sampling point and the obtained theoretical travel time, the theoretical arrival time of the seismic wave generated from the initial sampling point to that station is calculated. The observed arrival time of the seismic wave generated from the initial sampling point to that station is obtained from the station, and the difference between the theoretical arrival time and the observed time is calculated to obtain the travel time residual.

[0083] In some methods, the travel time table is constructed as follows: based on a preset one-dimensional velocity model, the travel time from any point at depth to any point on the surface is calculated, generating a travel time table of all points at all depths to all points on the surface. Since a one-dimensional velocity model is used, differences in horizontal position do not affect the arrival time calculation. A preset horizontal layered velocity model can be used to calculate the arrival time; the horizontal velocity value should be greater than the distance from the earthquake event to the farthest station. The above is only a general method for constructing the travel time table; this embodiment does not provide a detailed explanation of the construction principle and specific details of the travel time table.

[0084] S103: Based on the distribution of target values ​​at each initial sampling point, determine the local area where the initial sampling point with a target value greater than the preset target threshold is located;

[0085] In this embodiment, after calculating the target value of all initial sampling points, the local area where the target value is greater than the target threshold is determined according to the distribution of the target value. That is, after the initial round of search, a portion of the initial sampling points with higher target values ​​are selected from all the initial sampling points. Based on the location distribution of these initial sampling points, the local area where the earthquake occurred is preliminarily determined.

[0086] In some implementations, determining the local region where the initial sampling point where the target value is greater than the target threshold is located includes:

[0087] The longitude, latitude, depth, and time of the initial sampling points where the target value is greater than the target threshold are within the quantile range of the longitude, latitude, depth, and time dimensions, respectively.

[0088] The local region is determined based on the quantile range.

[0089] In this embodiment, based on the target values ​​of all initial sampling points, multiple initial sampling points with higher target values ​​are obtained. The local region for the next sampling is determined based on the quantiles of the source parameters of the high-value sampling points. Specifically, for initial sampling points with target values ​​greater than the target threshold, the upper and lower quantiles of these initial sampling points in the longitude, latitude, depth, and time dimensions are calculated based on their longitude, latitude, depth, and time values. Based on the upper and lower quantiles in each dimension, the quantile ranges in the longitude, latitude, depth, and time dimensions are determined. Then, based on the quantile ranges in the four dimensions, the local region containing the optimal solution of the objective function is determined.

[0090] S104: Generate multiple local sampling points within a local area;

[0091] S105: Calculate the target value of each local sampling point. Based on the target value of each local sampling point, find the optimal sampling point where the objective function reaches the maximum target value, and determine the source parameters and matching phase corresponding to the optimal sampling point.

[0092] In this embodiment, after determining the local region, several local sampling points are generated within the local region using the Sobol sampling method. These local sampling points serve as the solution space for a secondary search, and the optimal sampling point is then sought within this solution space. The method for finding the optimal sampling point is as follows: The target value of each local sampling point is calculated; based on the distribution of the target values ​​of each local sampling point, a sub-region is determined where the target value of the local sampling point is greater than the target threshold; multiple sub-region sampling points are generated within each sub-region; the target value of each sub-region sampling point is calculated; and the search area is gradually narrowed according to the above process until a preset stopping condition is reached, thus determining the maximum target value and the corresponding optimal sampling point.

[0093] This embodiment employs an iterative search method. First, approximately uniformly distributed initial sampling points are randomly generated within the global source parameter search area. The target value of each initial sampling point is calculated. Based on the distribution characteristics of the target value across the four dimensions of the source parameters, a local region where the optimal sampling point may exist is identified. Based on this local region, several new local sampling points are generated, and the target value of each local sampling point is calculated. The sub-region where the optimal sampling point is located is identified based on the distribution characteristics of the target value. Several sub-region sampling points are then generated within the sub-region, and the target value of each sub-region sampling point is calculated. The search range for the optimal sampling point is narrowed based on the distribution characteristics of the target value. Following this iterative search process, the search range is gradually narrowed. When the set iteration stopping condition is met, the search stops. Within the current search range, the maximum target value of the objective function, the optimal sampling point corresponding to the maximum target value, and the source parameters and seismic correlation results corresponding to the optimal sampling point are selected, thereby achieving seismic correlation and earthquake location.

[0094] In some methods, to improve search efficiency and accelerate the convergence process of iteration, the density of local sampling points generated in a local area is greater than the density of initial sampling points generated in the source parameter search area. In the subsequent search process of narrowing the search range, the number and density of sampling points suitable for the corresponding search area size can be determined according to the size of the search area.

[0095] In some approaches, to improve computational efficiency, the target value of the sampling points is calculated using a multi-dimensional matrix according to the objective function. By processing multiple sampling points in parallel using the PyTorch framework, computational efficiency can be significantly improved.

[0096] The seismic correlation and earthquake location method provided in this embodiment constructs an objective function based on the constraints of the number of associated seismic earthquakes and location reliability. This objective function not only focuses on the number of associated earthquakes but also adds constraints on earthquake location, enabling more accurate earthquake location estimates. During the seismic correlation and location process, the search range is gradually narrowed within a large source parameter search area to find the optimal solution of the objective function. This achieves the goal of exploring the global optimum while strengthening the search for local optima. Compared to existing methods, this application can automatically determine the search range based on data quality, flexibly addressing different correlation problems. For example, when data quality is poor, the target values ​​of sampling points are generally low, and some noisy points may be considered as sampling points with high target values. Sampling points with larger target values ​​are sparsely distributed, resulting in a larger search area. Conversely, when data quality is good, sampling points with larger target values ​​are more concentrated, resulting in a smaller search area. Experimental verification shows that the method of this application not only achieves a high recall rate but also provides more accurate and reliable earthquake location results.

[0097] In some embodiments, the objective function of this application may be subject to various constraints based on actual needs, conditions and objectives, thereby constructing diverse objective functions.

[0098] In one scenario, when a phase-picking model based on a neural network can be used to output the probability value of phase picking, and the obtained probability value is relatively reliable, this probability value can be used as a quantitative indicator of the confidence level of the phase result. A picking probability constraint is introduced into the objective function, and the phases of all associated stations are weighted and averaged using the picking probability to construct a phase quantity term that matches the source parameters, expressed as:

[0099]

[0100] in, Let W represent the phase pickup probability values ​​of the P-wave and S-wave at the i-th station, respectively; P and W S The weights are S, and the sum of the two weights is 1. N ∈[0,1].

[0101] The second scenario involves stations simultaneously associated with both P-waves and S-waves. In this case, the simultaneous association of a station with both P-waves and S-waves is used as a crucial criterion for reliability assessment. This involves introducing co-station phase constraints and constructing a phase quantity term that matches the source parameters, expressed as:

[0102]

[0103] Where, N B The number of stations that are simultaneously associated with both P-wave and S-wave.

[0104] The third scenario involves comprehensively considering both the pick probability constraint and the co-station facies constraint, constructing a facies quantity term that matches the source parameters, expressed as:

[0105]

[0106] The fourth scenario considers that amplitude helps improve correlation and positioning accuracy. An amplitude constraint is introduced, converting the amplitude into magnitude, and a positioning reliability term is constructed, expressed as:

[0107]

[0108] Among them, M i The magnitude of a single earthquake is calculated using conventional methods based on amplitude, epicentral distance, and gauge function. This is the average magnitude calculated based on the magnitudes of all individual earthquakes.

[0109] The effectiveness of the method of this application is illustrated below with reference to specific embodiments.

[0110] Using the method of this application, a model of M in Jishishan County, Gansu Province was constructed. L 6.2 Earthquake Aftershock Sequence. Seismic waveforms observed continuously for 14 days after the earthquake were collected. A phase picking model based on a neural network was used to automatically pick phases. Seismic correlations were performed using the method described in this application and compared with existing methods such as REAL, GaMMA, and PyOcto.

[0111] During this period, a total of 927 seismic events were manually identified, with REAL, PyOcto, GaMMA, and this application associating 1097, 1285, 1143, and 1200 earthquakes, respectively. Near the epicenter, a total of 909 events were manually recorded, with REAL, PyOcto, GaMMA, and this application associating 1064, 1184, 1080, and 1143 earthquakes, respectively. The number of seismic events associated by the four methods was significantly higher than that of the manual catalog. To further evaluate the consistency between the seismic catalogs constructed by each method, the recall rates of the seismic catalogs were calculated using both the manually constructed catalog and the catalogs constructed by the four methods as references, as shown in Table 1.

[0112] Table 1 Recall rates among different directories

[0113]

[0114] In Table 1, the horizontal axis represents the reference catalog, and the vertical axis represents the evaluation catalog. The values ​​represent the recall rate of the evaluation catalog for events in the reference catalog under the set matching conditions. The results show that this application achieves the highest recall rate when using manually constructed, REAL, and PyOcto-constructed seismic catalogs as references. In summary, the seismic catalog constructed in this application is not only numerous but also has a relatively high recall rate, indicating that the seismic catalog constructed in this application is relatively reliable.

[0115] Regarding the spatial distribution characteristics of earthquakes, since both REAL and PyOcto are based on grid search strategies, their distribution of earthquake events in the epicenter region exhibits a clear and regular arrangement. This application introduces focal location constraints during the earthquake correlation process, resulting in a more continuous and concentrated distribution of earthquake events, with the overall spatial morphology more closely resembling the actual characteristics of regional aftershock activity. To quantify the differences in focal parameter estimation among different earthquake catalogs, using an artificial earthquake catalog as a reference, the deviation distribution of earthquake catalogs constructed by the four methods in terms of origin time and focal location is calculated. In terms of time estimation, the standard deviations of the four methods are similar, but this application has the smallest average error, indicating that its origin time estimation is more accurate; in terms of focal location estimation, this application has the smallest error amplitude, outperforming other methods. It should be noted that this application achieved the highest recall rate in comparison with the artificial earthquake catalog, and most of the events in the multiple correlations were events with weak facies constraints. Therefore, compared to existing methods, this application can obtain more accurate focal parameters while maintaining high recall.

[0116] To evaluate the robustness of the proposed method under different noise environments, a synthetic test dataset was constructed based on the aftershock sequence on the first day after the 2023 Jishishan mainshock. This dataset contains 232 standard events and 10 seismic stations, including several events with occurrence times differing by only a few seconds, significantly increasing the challenge of seismic correlation tasks. Based on a regional velocity model, the theoretical arrival time of each seismic event to each station was calculated, constructing a standard phase set containing 2320 standard P-phases and S-phases. To simulate uncertainties in the actual phase picking process, different perturbations were introduced into the standard phases.

[0117] Based on the constructed test dataset, the performance of four methods was evaluated under different noise levels. Evaluation metrics included recall, accuracy, and F1 score compared to standard events; recall, accuracy, and F1 score compared to standard seismic phases; and deviation of source parameters from standard events. Experimental results show that this application exhibits excellent robustness under various noise conditions, with the recall and accuracy of its earthquake events and phases consistently approaching 100%, almost unaffected by disturbances. In contrast, the GaMMA and PyOcto methods are more sensitive to noise, with all metrics significantly decreasing as noise levels increase. Due to computational efficiency limitations, the REAL method used a large grid spacing in the tests, resulting in poor performance under high noise conditions.

[0118] Regarding source parameter estimation, this application consistently demonstrates lower overall deviation than other methods even with increasing noise levels, particularly in the estimation of epicenter location and time of origin, where errors remain at a low level, indicating strong anti-interference capabilities. Compared to existing methods, this application not only achieves more accurate seismic correlation but also possesses higher source parameter estimation accuracy. The improved positioning capability further enhances the accuracy and reliability of the correlation results.

[0119] Under noise levels of 0-5, the objective function was tested under different constraints. The results show that the recall of earthquake events and phases remains close to 100% under all constraints, stably associating all standard events with phases. Furthermore, the source parameter error remains relatively stable across all noise levels, without significant increase with noise enhancement, indicating good robustness under high noise conditions and different objective functions. When the noise level exceeds 2, introducing a picking probability constraint into the objective function effectively improves accuracy. Different constraints have little impact on recall; adding a picking probability constraint increases the reliability of the association, and adding seismic phases from the same station associates more earthquake events.

[0120] Regarding computational efficiency, in the practical application of the Jishishan earthquake sequence, the test platform was configured with a 20-core CPU. In the experiment, this application used a single-core CPU, while PyOcto and GaMMA were computed in parallel using 20 cores. Test results show that REAL took the longest time, while this application ran the fastest, reducing the time by approximately 17 times compared to its single-core CPU version. GaMMA, with its multi-core CPU, reduced the time by approximately 50%, and PyOcto by approximately 25%. Furthermore, the method in this application is insensitive to external parameters such as the number of sampling points, the optimal number of samples, the confidence interval, and the number of iterations, allowing for rapid migration and application to different regions, demonstrating good adaptability and versatility.

[0121] It should be noted that the method in this embodiment can be executed by a single device, such as a computer or server. The method can also be applied in a distributed scenario, where multiple devices cooperate to complete the task. In such a distributed scenario, one of these devices may execute only one or more steps of the method in this embodiment, and the multiple devices will interact with each other to complete the method described.

[0122] It should be noted that the above description describes specific embodiments of this specification. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims may be performed in a different order than that shown in the embodiments and still achieve the desired results. Furthermore, the processes depicted in the drawings do not necessarily require the specific or sequential order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0123] like Figure 2 As shown in the illustration, this application also provides a seismic correlation and earthquake location device, comprising:

[0124] The initial sampling point generation module is used to generate multiple initial sampling points within a preset source parameter search area;

[0125] The calculation module is used to calculate the target value of each initial sampling point according to a preset objective function; wherein, the objective function is constructed based on the number of seismic phases matched with the source parameters and the positioning reliability of the source parameters;

[0126] The local region determination module is used to determine the local region where the target value of the initial sampling point is greater than the preset target threshold, based on the distribution of the target values ​​of each initial sampling point.

[0127] The local sampling point generation module is used to generate multiple local sampling points within a local area.

[0128] The optimal solution determination module is used to calculate the target value of each local sampling point, and based on the target value of each local sampling point, find the optimal sampling point where the objective function reaches the maximum target value, and determine the source parameters and matching phase corresponding to the optimal sampling point.

[0129] For ease of description, the above devices are described in terms of function, divided into various modules. Of course, in implementing the embodiments of this application, the functions of each module can be implemented in one or more software and / or hardware.

[0130] The apparatus described above is used to implement the corresponding methods in the foregoing embodiments and has the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0131] Figure 3 This embodiment illustrates a more specific hardware structure of an electronic device, which may include a processor 1010, a memory 1020, an input / output interface 1030, a communication interface 1040, and a bus 1050. The processor 1010, memory 1020, input / output interface 1030, and communication interface 1040 are interconnected internally via the bus 1050.

[0132] The processor 1010 can be implemented using a general-purpose CPU (Central Processing Unit), microprocessor, application-specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of this specification.

[0133] The memory 1020 can be implemented in the form of ROM (Read Only Memory), RAM (Random Access Memory), static storage device, dynamic storage device, etc. The memory 1020 can store the operating system and other applications. When the technical solutions provided in the embodiments of this specification are implemented by software or firmware, the relevant program code is stored in the memory 1020 and is called and executed by the processor 1010.

[0134] The input / output interface 1030 is used to connect input / output modules to realize information input and output. Input / output modules can be configured as components within the device (not shown in the figure) or externally connected to the device to provide corresponding functions. Input devices may include keyboards, mice, touchscreens, microphones, various sensors, etc., while output devices may include displays, speakers, vibrators, indicator lights, etc.

[0135] The communication interface 1040 is used to connect a communication module (not shown in the figure) to enable communication between this device and other devices. The communication module can communicate via wired means (such as USB, Ethernet cable, etc.) or wireless means (such as mobile network, WIFI, Bluetooth, etc.).

[0136] Bus 1050 includes a pathway for transmitting information between various components of the device, such as processor 1010, memory 1020, input / output interface 1030, and communication interface 1040.

[0137] It should be noted that although the above-described device only shows the processor 1010, memory 1020, input / output interface 1030, communication interface 1040, and bus 1050, in specific implementations, the device may also include other components necessary for normal operation. Furthermore, those skilled in the art will understand that the above-described device may only include the components necessary for implementing the embodiments of this specification, and not necessarily all the components shown in the figures.

[0138] The electronic devices described above are used to implement the corresponding methods in the foregoing embodiments and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0139] The computer-readable medium of this embodiment includes permanent and non-permanent, removable and non-removable media, and information storage can be implemented by any method or technology. Information can be computer-readable instructions, data structures, program modules, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic magnetic disk storage or other magnetic storage devices, or any other non-transfer medium that can be used to store information accessible by a computing device.

[0140] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of this disclosure (including the claims) is limited to these examples; within the framework of this disclosure, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of the embodiments of this application as described above, which are not provided in the details for the sake of brevity.

[0141] Additionally, to simplify the description and discussion, and to avoid obscuring the embodiments of this application, the well-known power / ground connections to integrated circuit (IC) chips and other components may or may not be shown in the provided drawings. Furthermore, the apparatus may be shown in block diagram form to avoid obscuring the embodiments of this application, and this also takes into account the fact that the details of the implementation of these block diagram apparatuses are highly dependent on the platform on which the embodiments of this application will be implemented (i.e., these details should be fully understood by those skilled in the art). While specific details (e.g., circuits) have been set forth to describe exemplary embodiments of this disclosure, it will be apparent to those skilled in the art that the embodiments of this application can be implemented without these specific details or with variations thereof. Therefore, these descriptions should be considered illustrative rather than restrictive.

[0142] Although this disclosure has been described in conjunction with specific embodiments thereof, many substitutions, modifications, and variations of these embodiments will be apparent to those skilled in the art from the foregoing description. For example, other memory architectures (e.g., dynamic RAM (DRAM)) may be used with the embodiments discussed.

[0143] The embodiments of this application are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of the appended claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the embodiments of this application should be included within the protection scope of this disclosure.

Claims

1. A method of associating and locating earthquakes, characterized by, The method comprises the steps of: generating a plurality of initial sampling points within a preset search area of source parameters; calculating a target value of each initial sampling point according to a preset target function; wherein the target function is constructed according to the number of matched seismic phases and the positioning reliability of the source parameters; determining a local area where the initial sampling points with a target value greater than a preset target threshold value are located according to the distribution of the target values of the initial sampling points; generating a plurality of local sampling points within the local area; calculating the target value of each local sampling point, finding the optimal sampling point where the target function reaches the maximum target value based on the target values of the local sampling points, and determining the source parameters corresponding to the optimal sampling point and the matched seismic phases.

2. The method of claim 1, wherein, Further comprising: taking the nearest station observing P waves as an initial reference point, expanding a predetermined spatial search range in the horizontal direction, determining a time search range according to the average velocity of P waves and the preset search range in the horizontal direction, and determining a depth search range of the source according to the geological structure characteristics and the distribution of seismic activity; determining the source parameter search area according to the spatial search range, the time search range and the depth search range.

3. The method of claim 2, wherein, Within the preset search area of source parameters, a plurality of initial sampling points are generated, including: determining the number of initial sampling points according to the source parameter search area; generating initial sampling points with the number and having the dimension of source parameters by using the Sobol sampling method.

4. The method of claim 3, wherein, According to the preset target function, the target value of each initial sampling point is calculated, including: determining the travel time residual between the source at each initial sampling point and any station; if the travel time residual is less than a preset residual threshold, it is determined that the initial sampling point matches the station; counting the matching results of all initial sampling points and stations; calculating the target value of the target function according to the matching results.

5. The method of claim 4, wherein, The target function is: f(x1,x2,x3,t) = W N • S N (x1,x2,x3,t) + W L • S L (x1,x2,x3,t) where S N (x1,x2,x3,t) is the number of matched phases under the source parameter (x1,x2,x3,t), S L (x1,x2,x3,t) is the positioning reliability of the source parameter (x1,x2,x3,t), W N , W L is the weight, x1 is the longitude, x2 is the latitude, x3 is the depth, and t is the time of occurrence.

6. The method of claim 5, wherein, The calculation method of the number of seismic phases is: where N P is the number of P-waves matched to the sample point, N S is the number of S-waves matched to the sample point, W P , W S is the weight, N station is the total number of stations.

7. The method of claim 6, wherein, The calculation method of the positioning reliability is: wherein n is the number of stations associated with the sampling point; is the travel time residual of the P-wave corresponding to the i-th station, is the travel time residual of the S-wave corresponding to the i-th station; W i is the weight calculated by the i-th station according to its distance from the sampling point; and MaxTolerance is a preset maximum tolerance window.

8. The method of claim 1, wherein, According to the distribution of the target values of the initial sampling points, the local area where the initial sampling points with a target value greater than a preset target threshold value are located is determined, including: counting the longitude, latitude, depth and time of the initial sampling points with a target value greater than the target threshold value in the quantile range of the longitude dimension, the latitude dimension, the depth dimension and the time dimension; determining the local area according to the quantile range.

9. The method of claim 1, wherein, Based on the target values of the local sampling points, the optimal sampling point where the target function reaches the maximum target value is found, including: determining the sub-area where the local sampling points with a target value greater than the target threshold value are located based on the size of the target values of the local sampling points; generating a plurality of sub-area sampling points within the sub-area; calculating the target value of each sub-area sampling point; gradually reducing the search area according to the process until a preset stop condition is reached to determine the maximum target value and the corresponding optimal sampling point.

10. A device for associating and locating earthquakes, characterized in that, It comprises: an initial sampling point generation module for generating a plurality of initial sampling points within a preset search area of source parameters; a calculation module for calculating the target value of each initial sampling point according to a preset target function; wherein the target function is constructed according to the number of matched seismic phases and the positioning reliability of the source parameters; The local region determining module is configured to determine a local region in which the initial sampling points with the target values greater than the preset target threshold are located according to a distribution of the target values of the initial sampling points. The local sampling point generating module is configured to generate a plurality of local sampling points in the local region. The optimal solution determining module is configured to calculate the target value of each local sampling point, find an optimal sampling point in which the target function reaches a maximum target value based on the target values of the local sampling points, and determine a corresponding source parameter and a matched seismic phase of the optimal sampling point.

Images