Earthquake event reliability assessment method

Abstract

The invention relates to the technical field of seismic data processing, in particular to a seismic event reliability assessment method. The method comprises the following steps: firstly, constructing a holographic seismic event sample, dividing a seismic phase into a near field and a far field by using a near-stage influence range model, executing differential travel time constraint, and eliminating physical noise; and then, performing multi-dimensional coupling cleaning on the seismic phases by adopting iterative mean deviation analysis, linear regression fitting and spatial topology detection technologies, and generating a high-confidence data set. And on the basis, calculating a spatial distribution score of the event, performing reliability judgment in combination with a dynamic threshold model based on station sparse degree perception, and introducing a dispersion analysis mechanism based on a root-mean-square error to perform secondary verification on the edge event. According to the method, the problems that existing single threshold evaluation is poor in adaptability to a non-uniform station network and sensitive to abnormal values are solved, false earthquake events detected by artificial intelligence can be effectively filtered, and the accuracy and robustness of an earthquake directory are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of earthquake data processing technology, specifically to a method for assessing the reliability of earthquake events. Background Technology

[0002] With the rapid development of earthquake monitoring technology, the amount of earthquake data is growing exponentially. Currently, automatic earthquake data processing systems based on artificial intelligence technology have been widely used in earthquake event detection and analysis, greatly improving data processing efficiency.

[0003] However, in practical applications, existing automatic processing technologies still have the following significant defects and shortcomings: First, insufficient ability to identify false earthquake events. While AI-based automatic earthquake data processing systems can detect weak signals with high sensitivity, their automatically generated preliminary earthquake catalogs may contain some non-real events caused by noise mispicking or unreasonable correlation, resulting in false information in the final generated earthquake event catalog. This seriously affects the accuracy and reliability of subsequent earthquake research (such as earthquake location, mechanism inversion, etc.). Second, fixed evaluation criteria and lack of environmental adaptability. Existing technologies typically use traditional evaluation methods based on a single threshold when assessing the reliability of detected events. However, the actual earthquake monitoring environment is extremely complex, with uneven distribution of station density (such as dense and sparse areas), varying earthquake event sizes, and inconsistent data quality. Such fixed evaluation criteria are difficult to adapt to these changes, leading to the possibility of missing real events in sparse areas due to excessively high thresholds, or failing to effectively filter false events in dense areas due to excessively low thresholds. Third, sensitivity to outliers and weak anti-interference ability. Existing evaluation methods are quite sensitive to outliers in the data. When monitoring data contains noise interference, incorrect correlation of seismic phases, or anomalies in data from individual stations, a single evaluation index often fluctuates significantly, leading to misjudgments of the overall reliability of the seismic event. Fourth, there is insufficient tolerance for marginal events. For seismic events occurring at the edge of the monitoring network, or events with slightly missing data due to equipment failure but with a reliable overall trend, the existing one-size-fits-all evaluation mechanism lacks necessary tolerance. This results in many valuable real-world events being incorrectly excluded simply because they are at the edge of the monitoring network or have slightly higher data dispersion.

[0004] In summary, overcoming the limitations of existing single fixed threshold evaluation standards that are sensitive to uneven station distribution and data anomalies, and thus accurately identifying and effectively filtering false earthquake events detected by artificial intelligence in complex monitoring environments, is a technical problem that urgently needs to be solved in the field of earthquake data processing.

[0005] Therefore, a method for assessing the reliability of earthquake events is proposed. Summary of the Invention

[0006] The purpose of this invention is to provide a method for assessing the reliability of seismic events. First, a holographic seismic event sample is constructed. Using a near-station influence range model, seismic phases are divided into near-field and far-field regions, and differentiated travel time constraints are applied to eliminate physical noise. Subsequently, iterative mean deviation analysis, linear regression fitting, and spatial topology detection techniques are employed to perform multi-dimensional coupling cleaning of the seismic phases, generating a high-confidence dataset. Based on this, a spatial distribution score for the event is calculated, and a reliability determination is made using a dynamic threshold model based on station sparsity awareness. Furthermore, a dispersion analysis mechanism based on root mean square error is introduced to perform secondary verification of marginal events. This invention solves the problems of poor adaptability to non-uniform station networks and sensitivity to outliers in existing single-threshold assessments. It can effectively filter out false seismic events detected by artificial intelligence-based automatic seismic data processing systems, significantly improving the accuracy and robustness of the seismic catalog.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method for assessing the reliability of seismic events, comprising: Acquire seismic event data, including spatiotemporal parameters, seismic phases, and station information; calculate epicentral distance and observation travel time based on seismic event data; A near-site influence range model is constructed based on epicentral distance and near-site determination threshold. Seismic phases are divided into near-field and far-field phases according to epicentral distance and observation travel time. The near-field and far-field phases are filtered for consistency by matching the corresponding travel time window constraints. After removing noise, an effective set of phases and travel time residuals are generated. A cascaded screening process is performed on the effective seismic phase set, and outlier seismic phases are removed through iterative mean deviation analysis. A regression model of travel time residuals and epicentral distance is established to eliminate abnormal seismic phases that deviate from the regression trend line. A station spatial coverage topology is constructed to detect the gradient change of the epicentral distance reception ratio, and far-end seismic phases with a sudden drop in ratio are removed to generate a high-confidence seismic phase subset. Calculate the spatial distribution score of the high-confidence seismic phase subset; construct a dynamic threshold model using the number of seismic phases as an inverse variable, combined with the upper limit of the maximum threshold and the lower limit of the minimum threshold, to generate the pass threshold; determine the reliability of the earthquake event based on the spatial distribution score and the pass threshold.

[0008] Preferably, event information, station information, and seismic phase information are obtained and integrated by querying the business database or parsing earthquake reports and monitoring waveform files to construct a standardized dataset containing seismic phase objects; for each seismic phase object in the standardized dataset, the following parameters are calculated: Extract the longitude and latitude of the station to which the seismic phase belongs, as well as the longitude and latitude of the earthquake event, and calculate the surface arc length between the station location and the epicenter of the event as the epicentral distance; extract the actual arrival time of the seismic phase and the time of occurrence of the earthquake event, and calculate the time difference between the actual arrival time and the time of occurrence as the observation travel time.

[0009] Preferably, the specific steps for consistency filtering using the near-site influence range model include: setting a near-site determination threshold; for each seismic phase in the dataset, if the epicentral distance is greater than the near-site determination threshold, it is defined as a far-field seismic phase, and the observed travel time meets the following numerical range constraints: the observed travel time value is greater than the ratio of the epicentral distance to the upper limit of the P-wave propagation velocity, and less than the ratio of the epicentral distance to the lower limit of the P-wave propagation velocity; if the epicentral distance is not greater than the near-site determination threshold and any condition within a predetermined time is met, it is defined as a near-field seismic phase, and a relaxed near-field time window is adopted to allow the existence of direct waves and first waves; seismic phases falling outside the travel time window are removed; if the number of remaining seismic phases after removal is less than the number of associated stations, the event is marked as a false event and the evaluation is terminated.

[0010] Preferably, the calculation of the travel time residuals of the remaining seismic phases includes: Theoretical travel time acquisition: For each seismic phase, the theoretical arrival time is queried from a pre-set travel time table based on the epicentral distance and focal depth; if it exceeds the coverage of the travel time table, the TauPyModel calculation module is called, the Earth velocity model is selected, and the corresponding theoretical arrival time is dynamically calculated using the ray tracing method. Phase type correction: If multiple theoretical phases have the same theoretical travel time at the same epicentral distance, the theoretical phase with the highest physical propagation probability at that distance is selected as the benchmark phase according to the phase dominance interval of the epicentral distance. Residual calculation: The travel time residual of a single seismic phase is obtained by calculating the difference between the observed travel time and the theoretical arrival time; Overall correction: The arithmetic mean of the absolute values ​​of the residuals of all individual seismic phases is used as the overall residual benchmark for the event.

[0011] Preferably, the specific method for removing outlier phases through iterative mean deviation analysis is as follows: Initialize the phase set and calculate the mean of the residuals of all phases in the set; traverse each phase in the set and calculate the mean of the residuals of the remaining set after assuming the removal of outlier phases; calculate the proportion of mean change caused by the removal of phases, which is the ratio of the absolute value of the difference between the original mean and the new mean to the absolute value of the original mean; determine whether the following two removal conditions are met simultaneously: Condition 1: The proportion of mean change is greater than a preset mean deviation threshold; Condition 2: The absolute value of the travel time residual of the phase is greater than a preset absolute residual threshold; If met, the phase is removed, and the next round of iterative calculation is performed until there are no more phases in the set that meet the removal conditions.

[0012] Preferably, establishing a regression model for travel time residuals and epicentral distance specifically includes: for the remaining seismic phases after mean deviation analysis, constructing a linear regression equation with epicentral distance as the independent variable and residual as the dependent variable, and determining the slope and intercept of the regression line; calculating the corrected residual for each seismic phase, wherein the corrected residual is the difference between the actual observed residual of the seismic phase and the predicted residual calculated according to the linear regression equation; calculating the mean of all corrected residuals and the average of the absolute deviations of all corrected residuals relative to the mean; defining the elimination criteria: first, detecting the number of seismic phase samples participating in the regression; if the number of samples is less than a preset statistical significance threshold (e.g., 5), then only the absolute threshold discrimination is performed: if the absolute value of the corrected residual of the seismic phase is greater than the preset absolute tolerance upper limit, then the seismic phase is eliminated; if the number of samples is not less than the statistical significance threshold, then statistical discrimination is performed: if the deviation of the corrected residual of the seismic phase from the mean is greater than a specified multiple of the preset standard deviation (e.g., 3 times), and the absolute value of the original residual of the seismic phase is greater than the preset absolute residual threshold, then the seismic phase is eliminated.

[0013] Preferably, constructing the station spatial coverage topology and detecting gradient changes includes: constructing an epicentral distance sequence: extracting a list of epicentral distances for all online stations and a list of epicentral distances for all actually detected seismic phases, and sorting them in ascending order; calculating a reception ratio sequence: for each epicentral distance point in the list of actually detected seismic phases, defining a forward-covering sliding window (e.g., a window length of 50km or containing the nearest N stations), and calculating the local reception ratio within the window, wherein the local reception ratio is the ratio of the number of actually detected seismic phases within the window distance range to the total number of online stations within the range; Sudden drop detection: If a distance point exists in the received ratio sequence where the corresponding local received ratio is less than a preset local density threshold, and the distance interval between this point and the preceding phase is significantly greater than the theoretical attenuation interval of seismic waves in that area, then all phases with an epicentral distance not less than the stated distance point are removed. Joint abrupt change detection: Calculate the difference in received ratios and the difference in epicentral distance intervals between two adjacent epicentral distances; if the decrease in the received ratio of adjacent phases exceeds a preset change threshold, and the interval between adjacent epicentral distances is greater than the average interval of all stations in the network, then a spatial fault is identified, and all distant phases following the fault are removed.

[0014] Preferably, the dynamic threshold model and its judgment logic are specifically expressed as follows: the spatial distribution score is the arithmetic mean of the local reception ratios of all earthquakes in the seismic phase set; the calculation logic of the dynamic threshold model is as follows: calculate the difference between the maximum upper threshold and the minimum lower threshold, multiply the difference by the convergence coefficient, divide the difference by the total number of seismic phases associated with the current event, and add the quotient to the minimum lower threshold; wherein, the convergence coefficient is used to adjust the rate at which the threshold changes with the number of seismic phases; the value of the maximum upper threshold is positively correlated with the sparsity of stations in the monitoring area; the value of the minimum lower threshold is negatively correlated with the density of stations in the monitoring area; the preliminary reliable judgment condition is: the spatial distribution score is not less than the dynamic threshold; calculate the arithmetic square root of the ratio of the coverage area of ​​the seismic event network to the number of stations to obtain the average station spacing; extract the value corresponding to the average station spacing based on the preset parameter mapping relationship, wherein the average station spacing value is positively correlated with the extracted maximum upper threshold and minimum lower threshold values.

[0015] Preferably, the specific implementation method for determining the reliability of an earthquake event based on the spatial distribution score and the pass threshold is as follows: when the spatial distribution score is less than the dynamic threshold, a secondary verification process is initiated; the root mean square error of the spatial distribution score sequence is calculated, where the root mean square error is the square root of the arithmetic mean of the sum of the squares of the differences between the reception ratio of each seismic phase and the overall spatial distribution score of the event; verification must pass if the following two conditions are met simultaneously: Condition 1: the root mean square error is less than half the difference between the upper limit of the maximum threshold and the lower limit of the minimum threshold, indicating that the station distribution is relatively uniform and it is not an extreme abnormal event; Condition 2: the spatial distribution score is not less than the relaxed threshold, where the relaxed threshold is defined as the ratio of the dynamic threshold minus the lower limit of the minimum threshold to the convergence coefficient; if both conditions are met simultaneously, the reliability status of the event is corrected to "reliable"; otherwise, the "false event" judgment is maintained.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. This invention effectively addresses the high false alarm rate of existing AI-based automatic earthquake data processing systems. While existing AI detection models boast high sensitivity, they often incorporate noise interference or falsely correlated events. This invention does not rely on a single indicator but constructs a rigorous "multi-dimensional filter" encompassing physical constraints and statistical cleaning. It introduces a near-station influence range model for initial physical consistency screening, eliminating obvious noise that violates travel time patterns. Furthermore, it utilizes machine learning algorithms such as mean deviation analysis, linear regression fitting, and spatial topology detection to deeply clean the seismic phase data, accurately identifying and removing "disguised" seismic phases with statistically abnormal characteristics. This progressive filtering mechanism significantly reduces the residual rate of false earthquake events, ensuring the final generated earthquake event catalog has extremely high physical reliability, providing a high-quality data foundation for subsequent scientific research such as earthquake location, focal mechanism inversion, and imaging of the Earth's internal structure.

[0017] 2. This invention overcomes the limitations of traditional methods that use a fixed threshold for a one-size-fits-all approach, greatly enhancing the system's environmental adaptability. In actual monitoring, the distribution of stations is often extremely uneven (e.g., dense and sparse areas coexist), and fixed standards easily lead to "missed reports in sparse areas" or "false reports in dense areas." This invention innovatively proposes a dynamic threshold model based on seismic phase density sensing. The model can automatically adjust the threshold based on the number of seismic phases associated with an event and the station density in the area. In sparsely populated areas, due to the scarcity of data sources, the algorithm automatically raises the threshold, requiring extremely high consistency in single-station detection to ensure the authenticity of the event; while in densely populated areas, the algorithm automatically lowers the threshold, allowing for the absence of local seismic phases, but relying on a large number of stations for statistical confirmation. This mathematical dynamic adjustment mechanism enables the method to be compatible with monitoring networks of different densities and scales simultaneously, achieving unbiased and accurate assessment of the reliability of earthquake events across the entire network.

[0018] 3. This invention significantly improves the system's fault tolerance for atypical real events by introducing a secondary verification mechanism based on root mean square error. In traditional assessments, earthquakes occurring at the edge of the seismic network or weak earthquakes with slightly missing data are often mistakenly "falsely identified" due to insufficient spatial coverage scores. This invention designs a "revival" logic: even if the overall score of an event is slightly low, if the dispersion of its score distribution is extremely low, indicating that although the data is scarce, its internal consistency is extremely high, the system will still determine it as a reliable event. This mechanism effectively solves the problem of identifying edge events caused by limitations in monitoring geometry, and maximizes the preservation of valuable real earthquake records while ensuring that the overall data quality is not degraded, demonstrating the dialectical unity between rigor and inclusiveness in the algorithm. Attached Figure Description

[0019] Figure 1A flowchart of a seismic event reliability assessment method provided in an embodiment of the present invention; Figure 2 A flow chart of a multi-dimensional cleaning process based on seismic phases provided in an embodiment of the present invention; Figure 3 The flowchart of the secondary verification logic provided in the embodiment of the present invention. Detailed Implementation

[0020] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] Please see Figures 1 to 3 This invention provides a method for assessing the reliability of earthquake events, the technical solution of which is as follows: A method for assessing the reliability of seismic events, comprising: Acquire seismic event data, including spatiotemporal parameters, seismic phases, and station information; calculate epicentral distance and observation travel time based on seismic event data; A near-site influence range model is constructed based on epicentral distance and near-site determination threshold. Seismic phases are divided into near-field and far-field phases according to epicentral distance and observation travel time. The near-field and far-field phases are filtered for consistency by matching the corresponding travel time window constraints. After removing noise, an effective set of phases and travel time residuals are generated. A cascaded screening process is performed on the effective seismic phase set, and outlier seismic phases are removed through iterative mean deviation analysis. A regression model of travel time residuals and epicentral distance is established to eliminate abnormal seismic phases that deviate from the regression trend line. A station spatial coverage topology is constructed to detect the gradient change of the epicentral distance reception ratio, and far-end seismic phases with a sudden drop in ratio are removed to generate a high-confidence seismic phase subset. Calculate the spatial distribution score of the high-confidence seismic phase subset; construct a dynamic threshold model using the number of seismic phases as an inverse variable, combined with the upper limit of the maximum threshold and the lower limit of the minimum threshold, to generate the pass threshold; determine the reliability of the earthquake event based on the spatial distribution score and the pass threshold.

[0022] Example 1 This invention provides a method for assessing the reliability of earthquake events, which runs on the automated processing server of a provincial or national earthquake network center. Deep learning-based automatic earthquake monitoring systems generate thousands of earthquake event trigger records daily in real time. Due to the deep learning model's misidentification of noise and unreasonable correlations, the automatically generated catalog contains approximately 30% false events. This embodiment, as a post-processing module of an artificial intelligence-based automatic earthquake data processing system, aims to automatically clean the catalog, retaining high-confidence, genuine earthquake events.

[0023] As one embodiment of the present invention, refer to Figure 1 A flowchart of a method for assessing the reliability of earthquake events, referring to... Figure 2 A multi-dimensional cleaning process flow chart based on seismic phases, refer to Figure 3 Flowchart of secondary verification logic.

[0024] Furthermore, by querying the business database or parsing earthquake reports and monitoring waveform files, event information, station information, and seismic phase information are obtained and integrated to construct a standardized dataset containing seismic phase objects; for each seismic phase object in the standardized dataset, the following parameters are calculated: Extract the longitude and latitude of the station to which the seismic phase belongs, as well as the longitude and latitude of the earthquake event, and calculate the surface arc length between the station location and the epicenter of the event as the epicentral distance; extract the actual arrival time of the seismic phase and the time of occurrence of the earthquake event, and calculate the time difference between the actual arrival time and the time of occurrence as the observation travel time.

[0025] Specifically, based on the preliminary automatic rapid reporting results, the initial epicenter location was 104.60°E, 34.88°N. After the earthquake, the seismic correlation program successfully identified the event, extracting 134 valid seismic phases, including the actual arrival times of P-waves and S-waves. For these 134 seismic phases, using the built-in Earth ellipsoid model, the surface arc length between the coordinates of each station and the initial epicenter coordinates was calculated one by one, accurately obtaining the epicentral distance of each station. Simultaneously, the machine-picked arrival time of each seismic phase was read, and the time of occurrence (21 hours, 21 minutes, and 54 seconds) was subtracted to obtain a standardized observation travel time sequence. This process transforms the original waveform metadata into a unified "distance-time" two-dimensional feature vector, constructing a standardized data foundation for subsequent physical and statistical analysis, ensuring the comparability of data recorded by different stations and instruments in the same dimension.

[0026] This invention achieves the transformation from multi-source heterogeneous monitoring data to a unified assessment object through standardized data acquisition and parameter calculation steps. In actual monitoring, station coordinate systems and map projection methods often differ, and direct processing can easily introduce geometric errors. This method uniformly uses the surface arc length to calculate the epicentral distance, eliminating the influence of projection deformation; at the same time, it converts absolute arrival time into relative observation travel time, so that data characteristics no longer depend on the absolute time axis, but focus on the physical propagation law. This preprocessing mechanism not only improves the computational efficiency of subsequent algorithms, but more importantly, eliminates systematic errors caused by clock synchronization differences or coordinate system confusion, providing a solid data foundation for subsequent high-precision physical consistency filtering and statistical cleaning, ensuring the objectivity and accuracy of the assessment results.

[0027] Coarse screening: Statistical analysis is performed on each seismic phase and the corresponding station in the seismic phase object set. If the number of seismic phases is less than the number of stations, the event is directly defined as a false event.

[0028] Further, the specific steps for consistency filtering using the near-site influence range model include: setting a near-site judgment threshold; for each seismic phase in the dataset, if the epicentral distance is greater than the near-site judgment threshold, it is defined as a far-field seismic phase, and the observed travel time meets the following numerical range constraints: the observed travel time value is greater than the ratio of the epicentral distance to the upper limit of the P-wave propagation velocity, and less than the ratio of the epicentral distance to the lower limit of the P-wave propagation velocity; if the epicentral distance is not greater than the near-site judgment threshold and any condition within a predetermined time is met, it is defined as a near-field seismic phase, and a relaxed near-field time window is adopted to allow the existence of direct waves and first waves; seismic phases falling outside the travel time window are removed; if the number of remaining seismic phases after removal is less than the number of associated stations, the event is marked as a false event and the evaluation is terminated.

[0029] To ensure the robustness of this method in actual network monitoring, especially addressing the edge false alarm problem that is prone to occur in artificial intelligence systems, the parameter settings of the near-field influence range model should follow the engineering principle of "near-field inclusion and far-field strict control." The specific configuration rules are as follows: First, the setting of the near-site judgment threshold (model boundary): This threshold serves as the switching point for the physical filtering strategy and is mainly determined based on the statistical characteristics of historical earthquake travel time residuals in the monitored area. In regional seismic networks, automatic seismic data processing systems based on artificial intelligence technology often incorrectly correlate seismic signals from different events, similar azimuths, or the same time period, generating false events. To effectively suppress such long-distance miscorrelation noise, strict travel time velocity constraints need to be implemented in the far-field region. However, in the near-field range, due to the lateral non-uniformity of the crustal structure, the influence of low-velocity overburden, and the superposition of direct waves, surface reflected waves (such as virtual reflections), and shallow refracted waves, the actual travel time often deviates significantly from the theoretical prediction of the regional one-dimensional reference model, resulting in larger travel time residuals. If strict far-field travel time tolerance is applied too early in this range, it is very easy to mistakenly reject true near-field phases, impairing positioning accuracy and detection capabilities.

[0030] Therefore, for stations in the near field, a lenient near field travel time tolerance window is adopted to accommodate local velocity structure deviations and waveform complexity; for stations in the far field, a strict far field physical filtering strategy is switched to strongly clean up long-distance noise such as cross-event false correlations.

[0031] Second, far-field velocity constraint parameters: For far-field phases exceeding the judgment threshold, due to the long propagation path, the actual P-wave travel time should conform to the stable propagation law of the crustal medium. The lower limit of the P-wave apparent velocity is set at 2.9 km / s to exclude abnormally slow travel times caused by insufficient modeling or picking errors of low-velocity sedimentary layers; the upper limit is set at 6.5 km / s to eliminate ultra-short travel times caused by phase mislabeling, station clock synchronization errors, or non-seismic transient interference.

[0032] Third, adaptation of the near-field time window: In the near-field region within the judgment threshold, considering the high heterogeneity of local velocity structures and the aliasing of first-arrival waveforms, it is not advisable to rely on a single velocity model for rigid segmentation. A dynamic travel time tolerance strategy is adopted: using the theoretical P-wave travel time calculated by the regional velocity model as the center, typical picking errors (e.g., ±1.5 seconds) are superimposed, or the upper and lower limits of the travel time are inferred based on a reasonable velocity range (e.g., 2.9–6.5 km / s), and the wider interval of the two is taken as the final window. This strategy ensures that, even if the source location is not yet accurate, the core near-station data is not mistakenly deleted during the initial event localization stage, thereby guaranteeing the reliability of subsequent precise localization and magnitude inversion.

[0033] The formula for the near-field influence range model is as follows: Far-field phase conditions: epicentral distance > near-site determination threshold, requiring observation travel time to meet the following requirements: Near-field phase conditions: epicentral distance The near-station determination threshold requires the observation travel time to meet the following conditions: or Where dis represents the epicentral distance, and 2.9 km / s and 6.5 km / s represent the lower and upper limits of the reasonable propagation speed of the P wave, respectively.

[0034] Seismic phases are first classified into near-field and far-field categories based on the epicentral distance and a preset near-station determination threshold. Each type of seismic phase only needs to meet the travel time tolerance window of its corresponding region to be retained, without needing to meet both near-field and far-field conditions simultaneously. The near-field time window is based on the theoretical P-wave travel time calculated from the regional one-dimensional velocity model and is dynamically expanded in conjunction with typical picking errors (such as ±1.5s) to accommodate travel time deviations caused by surface wave aliasing, low-velocity overburden, and station positioning uncertainties.

[0035] Specifically, based on the crustal velocity model of the northwestern region where this earthquake occurred, a near-station determination threshold radius of 350 kilometers was set. Centered on the initial epicenter (104.60°E, 34.88°N), all available stations within a 350-kilometer radius were scanned and identified. For near-field phases within this epicenter range, considering the complexity of direct waves and crustal reflections, a strict single velocity cutoff was not adopted; instead, a wider travel time window was used. For far-field phases beyond this range, their apparent velocities were strictly limited to between 2.9 km / s and 6.5 km / s. During this process, the initial 134 phases were bidirectionally verified using the IASP91 theoretical travel time model, and uncorrelated stations were scanned to select high-quality signals with a measured arrival time deviation within ±2 seconds of the theoretical arrival time. This step not only verified the physical rationality of the original phases but also successfully added 18 new phases conforming to physical laws, bringing the total number of phases involved in the assessment to 152.

[0036] This invention cleverly solves the problem of inconsistent propagation characteristics of seismic waves in the near and far fields by constructing a near-site influence range model. In the near field, an inclusive strategy is used to protect the true signal whose travel time is offset due to waveform complexity; in the far field, a high-confidence initial screening defense line is constructed based on the physical boundaries of the crustal velocity structure. This method significantly improves the accuracy and robustness of seismic phase identification, while greatly reducing the computational load and false alarm risk of subsequent algorithms.

[0037] Furthermore, the calculation of the travel time residuals of the remaining seismic phases includes: Theoretical travel time acquisition: For each seismic phase, the theoretical arrival time is queried from a pre-set travel time table based on the epicentral distance and focal depth; if it exceeds the coverage of the travel time table, the TauPyModel calculation module is called, the Earth velocity model is selected, and the corresponding theoretical arrival time is dynamically calculated using the ray tracing method. Phase type correction: If multiple theoretical phases have the same theoretical travel time at the same epicentral distance, given the epicentral distance and depth, the typical phase type that conforms to the theoretical travel time at that distance shall be selected as the reference phase. Residual calculation: The travel time residual of a single seismic phase is obtained by calculating the difference between the observed travel time and the theoretical arrival time; Overall event travel time residual assessment: Calculate the arithmetic mean of the absolute values ​​of the residuals of all individual seismic phases as the overall residual benchmark for the event.

[0038] Specifically, for the set of seismic phases compiled after initial physical screening, the theoretical arrival time of each phase was calculated based on the IASP91 model and focal depth. Considering the unique crustal thickness in Northwest China, a phase type correction was automatically applied during the theoretical value calculation to ensure that the benchmark phase conforms to the regional propagation patterns of Pg or Pn waves. Subsequently, the difference between the actual observed travel time and the theoretical arrival time of each phase was calculated. Statistical analysis revealed that the residuals of some phases in the initial screening set exhibited anomalies. At this point, the arithmetic mean of the absolute values ​​of the residuals of all phases was calculated to construct the overall residual benchmark for this event. This benchmark not only reflects the overall positioning deviation level of the event but also provides a dynamic reference standard for subsequently identifying individual outliers exceeding the reasonable error range, ensuring that the residual analysis is based on a unified and accurate physical model.

[0039] This invention standardizes the calculation path for theoretical travel time and residuals, effectively improving the scientific rigor of the evaluation benchmark. By introducing TauPyModel ray tracing and phase dominance interval correction, it solves the accuracy problem of the lookup table method under complex geological structures, avoiding misjudgments caused by errors in theoretical value calculations. Using the mean of all individual residuals as the overall benchmark allows the system to adaptively perceive the positioning quality level of the current event. This residual analysis method based on a dynamic benchmark, compared to traditional methods using fixed absolute values, more objectively reflects the true fit of seismic events under different geological environments, providing a reliable quantitative basis for subsequent accurate anomaly data removal.

[0040] Furthermore, the specific method for removing outlier phases through iterative mean deviation analysis is as follows: initialize the phase set, and calculate the mean of the residuals of all phases in the set. Iterate through each seismic phase in the set and calculate the hypothesis of removing outlier phases. Mean of the residuals of the remaining set ; Calculate the percentage change in mean caused by removing the seismic phase, which is the ratio of the absolute value of the difference between the original mean and the new mean to the absolute value of the original mean, i.e. Determine whether the following two elimination conditions are met simultaneously: Condition 1: The change ratio of the mean is greater than the preset mean deviation threshold; Condition 2: The absolute value of the travel time residual of the seismic phase is greater than the preset absolute residual threshold; If met, the seismic phase is eliminated, and the next round of iteration calculation is performed until there are no more seismic phases in the set that meet the elimination conditions.

[0041] The setting of the two rejection thresholds in this step follows the following robust statistical principles: 1. Mean change percentage threshold (10%): Based on the "outlier" theory in robust statistics. With a sample size greater than 30, if the removal of a single data point causes a mean shift exceeding 10%, that point is considered a disruptive outlier with a 95% confidence level. For small samples (<10), this threshold is relaxed to 20%-25%.

[0042] 2. Residual Absolute Threshold (3.0 seconds): This value is set to k times the baseline of the overall residuals of the event (e.g., k=3). Based on the assumption of normal distribution of errors, 3 times the standard deviation covers 99.7% of the normal fluctuation range, and residuals exceeding this range belong to low-probability abnormal events.

[0043] 3. Logical Relationship: The "AND" logic is used to reduce the false positive rate. A large mean change alone may be due to actual travel time jumps caused by geological structures, and a high residual may be due to special seismic phases; only when both conditions are met can it be confirmed as noise that is "statistically abnormal and physically unexplainable".

[0044] Specifically, a second round of in-depth quality control was initiated. First, iterative mean deviation analysis was performed, traversing all seismic phases and calculating the impact of removing a single point on the overall residual mean. The analysis revealed that 11 seismic phases had absolute residuals exceeding 3.0 seconds, and removing these phases significantly improved the stability of the set's mean. These 11 phases were determined to exceed a reasonable error threshold and did not meet the requirements for high-precision positioning, and were therefore removed. Subsequently, a linear regression model of travel time residuals and epicentral distance was established to check the consistency of the remaining seismic phases. Evaluation revealed that 11 other seismic phases, while their absolute residuals did not exceed the limit, significantly deviated from the overall linear trend line in the regression plot. These phases were determined to potentially have human picking errors or be affected by special path effects, disrupting spatial consistency, and were also removed. After these two steps of cleaning, a total of 22 abnormal seismic phases were removed, ultimately retaining a subset of 130 high-quality, highly consistent seismic phases.

[0045] This invention employs a cascaded screening strategy to achieve in-depth cleaning from statistical distribution to spatial trends. Iterative mean deviation analysis utilizes robust statistical principles to accurately identify and eliminate large residual outliers that have a "destructive" impact on overall accuracy. The linear regression model further captures systematic biases hidden within the normal residual range, eliminating noise that, while numerically small, violates the laws of physical propagation. This "two-pronged" cleaning mechanism completely solves the problem of a single indicator failing to account for both gross errors and systematic biases, ensuring that the final retained seismic phase set is not only numerically accurate but also spatially conforms to the linear attenuation law of seismic waves, greatly improving data purity.

[0046] Furthermore, establishing a regression model for travel time residuals and epicentral distance specifically includes: for the remaining seismic phases after mean deviation analysis, constructing a linear regression equation with epicentral distance as the independent variable and residual as the dependent variable, and determining the slope and intercept of the regression line; calculating the corrected residual for each seismic phase, wherein the corrected residual is the difference between the actual observed residual of the seismic phase and the predicted residual calculated according to the linear regression equation; calculating the mean of all corrected residuals and the average of the absolute deviations of all corrected residuals from the mean; defining the exclusion criterion: first, detecting the number of seismic phase samples participating in the regression. If the number of samples is less than the preset statistical significance threshold (set to 5), then only the absolute threshold discrimination is performed: if the absolute value of the corrected residual of the seismic phase is greater than the preset absolute tolerance upper limit, then the seismic phase is removed; if the number of samples is not less than the statistical significance threshold, then statistical discrimination is performed: if the deviation of the corrected residual of the seismic phase from the mean is greater than a specified multiple of the preset standard deviation (set to 3 times), and the absolute value of the original residual of the seismic phase is greater than the preset absolute residual threshold, then the seismic phase is removed.

[0047] It is worth noting that although the Earth's medium exhibits heterogeneity, in the far-field range exceeding the near-site determination threshold, travel time residuals and distance typically show a first-order linear correlation. To overcome the nonlinear effects of extreme geological environments, this step preferably performs an applicability test before model establishment: calculating the Pearson correlation coefficient between the residuals and the epicentral distance; if the absolute value of the correlation coefficient is less than 0.3, the linear elimination step is skipped. Furthermore, the random sampling consensus algorithm used in this embodiment does not rely on a strict global linear assumption, but rather iteratively finds the largest subset of consistent in-situ points. Therefore, even in the presence of local nonlinear structures, it can still accurately identify outliers that significantly deviate from the physical framework, avoiding the false deletion caused by outliers skewing the regression line in ordinary least squares methods.

[0048] Specifically, after removing gross errors, a refined linear consistency test was performed on the remaining phases. A linear regression model was constructed with epicentral distance as the X-axis and travel time residuals as the Y-axis. Although most phases met the requirements in terms of absolute error, the regression analysis detected 11 phases whose distributions significantly deviated from the overall linear regression trend line. Calculations showed that the corrected residuals of these phases (i.e., the deviation relative to the regression line) were significantly greater than three times the standard deviation of the mean of the corrected residuals. This means that although their absolute residual values ​​may not exceed the hard threshold of 3.0 seconds, they exhibit abnormal nonlinear characteristics in their spatial propagation patterns (possibly due to path effects or phase misjudgment). Based on this, these 11 "latent anomalies" were determined to have disrupted spatial consistency and were removed. Thus, a total of 22 anomalous phases (11 gross errors + 11 linear anomalies) were removed, retaining a subset of 130 high-confidence phases.

[0049] This invention utilizes a linear regression model to capture the systematic trend of residuals changing with distance, identifying and eliminating hidden outliers. In some cases, erroneous seismic phases may fall precisely near the mean, but spatially violate the propagation laws of seismic waves. By calculating corrected residuals and combining them with regression trend lines for screening, this step effectively removes statistically normal but physically illogical seismic phases, further improving data purity, preventing systematic biases from misleading the final scoring and location, and ensuring the spatial linear consistency of the retained data.

[0050] Furthermore, constructing the station spatial coverage topology and detecting gradient changes includes: constructing an epicentral distance sequence: extracting a list of epicentral distances for all online stations and a list of epicentral distances for all actually detected seismic phases, and sorting them in ascending order; calculating a reception ratio sequence: for each epicentral distance point in the list of actually detected seismic phases, defining a forward-covering sliding window (e.g., a window length of 50km or containing the nearest N stations), and calculating the local reception ratio within this window, whereby the local reception ratio is the ratio of the number of actually detected seismic phases within the window distance range to the number of online stations within the range; the local reception ratio formula is as follows: in, This represents the local reception ratio. This indicates the total number of all detected seismic phases. This indicates the number of all online stations. This represents the j-th epicentral distance of the actually detected seismic phase. This represents the distance from the k-th epicenter among the online stations. Indicates the range of window distances. It is an indicator function that returns 1 when the condition is true and 0 otherwise.

[0051] Perform sudden drop detection: If there is a distance point in the received ratio sequence whose corresponding local received ratio is less than the preset local density threshold, and the distance interval between this point and the preceding phase is significantly greater than the theoretical attenuation interval of the seismic waves in this area, then all phases with an epicentral distance not less than the distance point are removed.

[0052] Perform joint mutation detection: calculate the difference in the reception ratio and the difference in the epicentral distance between two adjacent epicentral distances; if the decrease in the reception ratio of adjacent epicentral distances exceeds the preset change threshold and the interval between adjacent epicentral distances is greater than the average interval of all stations in the network, then a spatial fault is determined to have occurred, and all far-end seismic phases after the fault are removed.

[0053] Specifically, a spatial topological scan was performed on the remaining seismic phase data after linear regression cleaning to eliminate false far-end correlations caused by network coverage faults. First, an ascending sequence of epicentral distances was constructed, and the station locations corresponding to the current 130 earthquakes were extracted. Combined with information from all online stations within a 350km radius of the epicenter, the local reception ratio within each distance window was calculated. For the Gansu Dingxi event, a forward-covering sliding window was defined, and the gradient change in density of each node was calculated sequentially. During the detection process, the difference in reception ratios and distance intervals between adjacent epicentral distances were analyzed. Calculations showed that these 130 seismic phases exhibited reasonable physical attenuation characteristics in their spatial distribution. The reception ratio sequence decreased smoothly with increasing distance, and no "spatial fault" phenomenon was detected, where the decrease in adjacent reception ratios exceeded a preset threshold (set at 20%) and was accompanied by large-span distance intervals (exceeding the average spacing of the entire network). This means that the seismic phase set is continuous in spatial topology, and there are no false signal chains accidentally pieced together by random noise from distant locations. Based on this, gradient detection was passed, and these 130 seismic phases were confirmed as the final high-confidence seismic phase subset, providing rigorous spatial geometric assurance for subsequent scoring calculations.

[0054] This invention introduces topological consistency detection from the perspective of station spatial coverage, effectively defending against "spatial fault type" spurious events. Existing AI-based automatic seismic data processing systems often mistakenly associate random noise from two distant stations with the distant waveform of a near-seismic event because the noise happens to be temporally aligned. These spurious phases typically cause a sharp drop in the reception ratio spatially. This algorithm, by detecting a sudden drop in the gradient of the epicentral distance reception ratio, can keenly identify this spatial discontinuity and automatically dismantle unreliable distant phase chains. This not only eliminates highly concealed far-field noise but also ensures that the retained phases have a reasonable physical envelope from a geometrical perspective, greatly enhancing the reliability of event localization and ensuring that the final generated phase subset has high physical self-consistency in both spatiotemporal dimensions.

[0055] Furthermore, the specific expression of the dynamic threshold model and its judgment logic is as follows: the spatial distribution score is the arithmetic mean of the local reception ratios of all earthquakes in the seismic phase set; the calculation logic of the dynamic threshold model is: calculate the difference between the maximum upper threshold and the minimum lower threshold, multiply the difference by the convergence coefficient, divide the difference by the total number of seismic phases associated with the current event, and then add the quotient to the minimum lower threshold; the threshold calculation formula is: in, To generate a pass threshold, The threshold is the total number of seismic phases associated with the current event; A is the convergence coefficient used to adjust the rate at which the threshold changes with the number of seismic phases. The value of A is typically set to the rounded-to-the-element of 40% of the average number of historical reliable seismic phases in the target monitoring area. The 40% ratio is chosen based on the principle of statistical significance: when the real-time number of seismic phases is less than half the average level, it is still in a small sample fluctuation zone, and the value of A ensures that the threshold remains at a relatively high level; only when the number of seismic phases exceeds this critical value does the threshold begin to decrease significantly. If historical data is missing, data from adjacent similar seismic networks can be referenced, or expert empirical values ​​can be used. For example, if the average number of seismic phases in the area's historical events is 20, the convergence coefficient can be set to 8 to ensure that the threshold curve smoothly transitions from high to low stringency within the statistically significant sample size range. This is the maximum threshold, and its value is positively correlated with the sparsity of stations within the monitoring area; The minimum threshold is set as the lower limit, and its value is negatively correlated with the density of monitoring stations within the monitoring area; the preliminary reliable judgment condition is that the spatial distribution score is not less than the passing threshold. and The value is determined by: calculating the ratio of the number of stations in the seismic event's network to the actual coverage area of ​​the network to obtain the station density; extracting the value corresponding to the station density based on a preset parameter mapping relationship, wherein the station density value is negatively correlated with the extracted maximum threshold and minimum threshold values, that is, taking the larger value in sparsely populated areas. and In areas with dense stations, take smaller values and This allows for adaptive evaluation of monitoring networks with different densities. To prevent threshold overflow under small sample sizes, a truncation operation is performed on the basic increment to ensure its value does not exceed (maximum threshold upper limit - minimum threshold lower limit). The truncated increment is then added to the minimum threshold lower limit to ensure that the final generated pass threshold never exceeds the maximum threshold upper limit, preserving necessary tolerance for edge events. Specifically, when the total number of seismic phases is less than 3, the pass threshold is directly locked to the maximum threshold upper limit, and the above division operation is no longer performed to avoid numerical instability.

[0056] The parameter mapping relationships are stored in a lookup table in the configuration file. For example, using the average station spacing (unit: kilometers) as a metric: when the average station spacing is greater than 150 kilometers (sparse region), the maximum threshold is set to 0.95 and the minimum threshold to 0.40; when the average station spacing is between 60 and 150 kilometers (medium region), the maximum threshold is set to 0.85 and the minimum threshold to 0.30; when the average station spacing is less than 60 kilometers (dense region), the maximum threshold is set to 0.75 and the minimum threshold to 0.20; the threshold parameters corresponding to the non-node station spacing are calculated using linear interpolation.

[0057] Specifically, for the 130 high-confidence seismic phases retained after the previous cleaning steps, the spatial distribution score of this set is first calculated. These 130 seismic phases are traversed, their corresponding local reception ratio values ​​are extracted, and the average is calculated to quantify the overall coverage quality of the event. Then, the dynamic threshold generation process begins. The network parameters of the monitoring area in Northwest China are read. Although the overall station density in this area is generally low, for this magnitude 3.x earthquake, due to the successful association of a large number of seismic phases (130), threshold parameters suitable for medium-to-high density scenarios are matched based on the parameter mapping relationship: the minimum threshold is set to 0.20, the maximum threshold is set to 0.75, and the convergence coefficient is set to 10. During the calculation, the total number of seismic phases N (130) acts as an inverse variable, greatly diluting the incremental term in the numerator ((0.75-0.20)×10÷130≈0.04). The final generated pass threshold (0.20+0.04=0.24) is very close to the minimum lower limit. The calculated spatial distribution score (e.g., 0.65) was compared with the dynamic threshold (0.24), and the result significantly met the standard, thus satisfying the preliminary reliability condition. This process demonstrates that for significant events with abundant and highly consistent data, the entry threshold can be automatically lowered to ensure their smooth passage.

[0058] This invention constructs an adaptive threshold model based on dual sensing of "number of seismic phases" and "station density," overcoming the bottleneck of traditional fixed thresholds in non-uniform seismic networks. Existing technologies typically employ a "one-size-fits-all" scoring standard, which can easily lead to false alarms due to noise in sparsely populated areas or missed alarms due to local gaps in densely populated areas. This model utilizes inverse proportional function logic to achieve intelligent threshold floating: when seismic phases are scarce, the threshold automatically approaches the upper limit to enforce a strict admission standard; when seismic phases are abundant, the threshold automatically lowers to the lower limit to establish a more lenient admission standard. This mechanism enables the system to be compatible with different density environments simultaneously, achieving accurate assessment of the reliability of seismic events across the entire network.

[0059] Further, based on the spatial distribution score and the pass threshold, the specific implementation method for determining the reliability of an earthquake event is as follows: when the spatial distribution score is less than the pass threshold, a secondary verification process is initiated; the root mean square error (RMSE) of the spatial distribution score sequence is calculated, where the RMSE is the square root of the arithmetic mean of the sum of the squares of the differences between the reception proportions of each seismic phase and the overall spatial distribution score of the event; verification success requires the simultaneous fulfillment of the following two conditions: Condition 1 (low dispersion): the RMSE is less than half the difference between the upper limit of the maximum threshold and the lower limit of the minimum threshold, indicating a relatively uniform station distribution and non-extreme anomaly; Condition 2 (relaxed threshold): the spatial distribution score is not less than the relaxed threshold, which is defined as the pass threshold minus the lower limit of the minimum threshold and the convergence coefficient, i.e., the formula: ;in, This is the slack quantity. If both of the above conditions are met simultaneously, the reliability status of the event is corrected to "reliable"; otherwise, the "false event" determination is maintained.

[0060] Specifically, although the 3.x magnitude earthquake in Northwest China directly passed the dynamic threshold assessment due to its extremely high data quality (130 seismic phases), a secondary verification logic was still executed synchronously in the background to ensure system robustness. The root mean square error (RMSE) of the received proportion sequence of these 130 seismic phases was calculated. Because of rigorous physical and statistical cleaning, the remaining seismic phases are spatially uniformly distributed, with minimal differences in detection rates across different distance segments. The calculated RMS error (e.g., 0.08) is far less than the system's set dispersion tolerance (i.e., (0.75-0.20)÷2=0.275). This extremely low dispersion strongly demonstrates the high consistency of the event's internal structure, ruling out the possibility of random noise coupling. Even if the event score is slightly lower, this excellent low dispersion characteristic can still satisfy the relaxed threshold conditions. Finally, based on these 130 high-quality data points that have undergone multiple verifications, the hypo2000 algorithm was used for relocation, correcting the initial epicenter (104.60°E, 34.88°N) to 104.60°E, 34.87°N.

[0061] The secondary verification mechanism based on root mean square error introduced in this invention adds a "safety net" to the evaluation system. In traditional evaluations, many real earthquakes located at the edge of the seismic network or in areas with weak monitoring capabilities are often mistakenly deleted due to their overall low scores. This method, by analyzing the dispersion of scores, can accurately identify high-quality events that, although their scores are low due to objective conditions, have extremely strong internal data consistency, and give them a chance to be re-reported with threshold relaxation. This not only effectively reduces the false negative rate of weak earthquakes at the edge, but also reflects the algorithm's balance between rigorous screening and scientific error tolerance, preserving earthquake records with research value to the greatest extent possible.

[0062] This invention completely eliminates the drawbacks of traditional single fixed threshold evaluation models by constructing a closed-loop evaluation system of "physical constraints - statistical cleaning - dynamic decision-making," organically combining geophysical propagation laws with data-driven statistical methods. This multi-dimensional evaluation strategy not only effectively identifies and eliminates false earthquake events generated by artificial intelligence models, but also significantly reduces sensitivity to noisy data while maintaining high accuracy. Ultimately, this method can output a high-purity earthquake catalog with strong physical consistency, significantly enhancing the practical value of automated processing systems.

[0063] Example 2 This embodiment demonstrates a method for seismic event reliability assessment and data cleaning deployed on the central server of a seismic monitoring network in Northwest my country. The monitoring area has complex topography and active geological structures. The distribution of monitoring stations is geographically constrained, exhibiting typical non-uniformity, and is frequently affected by wind noise and human activities. In actual operation, while AI-based automatic detection systems have greatly improved detection sensitivity, they have also brought problems such as high false trigger rates and large positioning errors. Especially when processing small to medium magnitude events, errors in phase picking or false associations often lead to a decline in catalog quality. This embodiment details how this invention utilizes a strategy combining physical constraints and statistical cleaning to accurately extract effective information from massive amounts of monitoring data and achieve an objective determination of event reliability by conducting a full-process assessment of a magnitude 3.x earthquake event occurring in Northwest China.

[0064] As one embodiment of the present invention, refer to Figure 1 A flowchart of a method for assessing the reliability of earthquake events, referring to... Figure 2 A multi-dimensional cleaning process flow chart based on seismic phases, refer to Figure 3 Flowchart of secondary verification logic.

[0065] First, triggering events to be evaluated are captured from the real-time stream processing queue via a dedicated data interface. Monitoring data shows that a magnitude 3.x earthquake occurred in Northwest China on the evening of January 3rd of a certain year. The initial epicenter location given by the automatic rapid reporting system is locked at 104.60°E, 34.88°N. In the initial state, the seismic correlation program identifies and extracts 134 valid seismic phases, including the actual arrival time information of P-waves and S-waves.

[0066] To eliminate discrepancies between instrument responses and time references at different stations, the 134 seismic phases were first standardized. Using the built-in WGS84 Earth ellipsoid model, the surface arc length between the coordinates of each associated station and the initial epicentral coordinates was calculated one by one, accurately obtaining the epicentral distance for each station. Simultaneously, the machine-picked time for each seismic phase was read, and the time of occurrence was subtracted to obtain a standardized observation travel time sequence. This process transforms the original waveform metadata into a unified "distance-time" two-dimensional feature vector, eliminating systematic errors caused by clock synchronization discrepancies or coordinate system confusion, and providing a solid data foundation for subsequent high-precision physical consistency filtering and statistical cleaning.

[0067] Given the significant crustal thickness and lateral variations in velocity structure in Northwest China, a single travel-time residual standard is insufficient to accommodate the complex wavefield characteristics. Therefore, the first round of physical consistency optimization was initiated. Based on the crustal velocity model for this region, a threshold radius of 350 kilometers was set for determining the influence range of near-stations.

[0068] A full network scan was conducted on all available stations within a 350-kilometer radius of the initial epicenter (104.60°E, 34.88°N). Within this physical space, a differentiated constraint strategy was implemented: for the near-field region within this range, considering the complex superposition of direct waves and crustal reflections, a strict single velocity cutoff was not adopted. Instead, a wider travel time window (e.g., a tolerance range of ±2 seconds) was matched using the IASP91 theoretical travel time model to accommodate reasonable travel time discrepancies caused by crustal heterogeneity; for the far-field region beyond this range, the apparent velocity was strictly limited to between 2.9 km / s and 6.5 km / s to strongly block interference from non-seismic signals.

[0069] This process not only verified the physical validity of the original 134 seismic phases, but also used the aforementioned physical model to conduct a "recall" scan of unrelated stations. After screening, 18 high-quality new seismic phases with measured arrival times deviating from theoretical arrival times within ±2 seconds were successfully added. This step increased the total number of seismic phases participating in the evaluation to 152. The newly added seismic phases significantly enhanced the azimuth coverage and distance gradient of the station distribution, filled some azimuth monitoring blind spots, and laid a richer physical foundation for subsequent precise positioning.

[0070] After the initial physical screening, a second round of in-depth quality control was initiated for the 152 collected seismic phases. This stage aimed to use statistical methods to eliminate "noise" that, while conforming to broad physical laws, exhibited abnormal statistical distribution.

[0071] First, an iterative mean deviation analysis was performed to eliminate gross errors. Based on the epicentral distance and focal depth, the TauPyModel module was used to calculate the theoretical arrival time for each seismic phase, resulting in an initial residual set. Subsequently, the algorithm entered an iterative loop, calculating the mean stability change of the remaining set after removing each seismic phase. Statistical analysis revealed that 11 seismic phases in the set had absolute residuals exceeding the reasonable error threshold of 3.0 seconds. Calculations showed that the presence of these seismic phases caused a significant shift in the overall residual mean, constituting statistically destructive outliers. These 11 seismic phases were determined to severely lower the overall positioning accuracy and not meet the requirements for high-precision positioning; therefore, they were removed as the first batch of outlier data.

[0072] Secondly, a linear consistency test was performed on the residual phases to eliminate latent anomalies. A linear regression model was constructed with epicentral distance as the independent variable and travel time residuals as the dependent variable. Although most phases met the requirements in terms of absolute error, 11 additional phases were detected in the regression analysis spectrum whose distributions significantly deviated from the overall linear regression trend line. Calculations showed that the corrected residuals of these phases (i.e., the deviations relative to the regression line) were significantly greater than three times the standard deviation of the mean of the corrected residuals. This means that although their absolute residual values ​​may not exceed the hard threshold of 3.0 seconds, they exhibit abnormal nonlinear characteristics in spatial propagation (most likely due to systematic deviations caused by special path effects or phase misjudgments). Based on this, these 11 "latent anomalies" were determined to have disrupted the spatial consistency of seismic wave propagation and were therefore eliminated.

[0073] After the two rigorous statistical cleaning steps described above, a total of 22 anomalous phases (11 gross errors + 11 linear anomalies) were removed, ultimately retaining a high-quality, highly consistent subset of phases. Subsequently, a spatial topological scan was performed on these 130 phases, confirming that their reception ratio changed smoothly with distance, without any abrupt drops in gradient, thus verifying the spatial continuity and integrity of the phase set.

[0074] Moving into the decision-making stage, a final assessment of the event's reliability is required. For the 130 high-confidence seismic phases ultimately retained, the average reception ratio across each distance segment is first calculated to obtain the event's spatial distribution score.

[0075] To address the issue that traditional fixed thresholds cannot adapt to non-uniform seismic networks, a dynamic threshold model is employed. First, based on the scanning radius (350 km) and the number of effective seismic phases (130), an inversion estimation is performed, calculating that the equivalent average inter-station spacing of the area affected by this earthquake is approximately 50-60 km, which falls within the "dense region" in the parameter mapping relationship. Based on this, the corresponding threshold parameters are automatically matched and extracted: a minimum threshold lower limit of 0.20, a maximum threshold upper limit of 0.75, and a convergence coefficient of 10.

[0076] In the threshold generation calculation, the total number of seismic phases N (130) associated with this event played a crucial role as an inverse proportional variable. According to the formula's logic, the large denominator significantly diluted the incremental term in the numerator, resulting in a very small dynamic increment. The final generated pass threshold was very close to the minimum threshold lower limit (approximately 0.24). This means that, given the abundant data and its high consistency after cleaning, the law of large numbers in statistics automatically lowered the pass threshold, acknowledging its reliability. Comparing the calculated spatial distribution score with this dynamically generated pass threshold showed that the score was significantly higher than the threshold. Therefore, the event was directly determined to meet the preliminary reliability criteria.

[0077] Although the event passed the initial assessment due to its extremely high data quality, a secondary verification based on the root mean square error (RMSE) was still performed concurrently in the background to ensure system robustness. The RMSE values ​​for the received proportion sequences of these 130 seismic phases were calculated. The results showed that the RMSE values ​​were extremely low (far below the system's set dispersion tolerance), strongly demonstrating the high consistency of the event's internal structure and completely ruling out the possibility of random noise coupling. This mechanism indicates that even for marginal events with smaller magnitudes and slightly lower scores, as long as they possess low-dispersion physical characteristics, they can still be accurately identified using a relaxed threshold, thereby effectively reducing the false negative rate.

[0078] Finally, based on these 130 high-quality seismic phase data points that had undergone physical constraints, statistical cleaning, and spatial verification, the hypo2000 localization algorithm was used to relocate the event. The localization results showed that the epicenter location was corrected to 104.60°E, 34.87°N. Compared with the initial automatic rapid report, the corrected epicenter location showed a slight adjustment in latitude, and the range of the error ellipse was significantly reduced. The event was marked as "reliable" and officially added to the earthquake catalog, providing high-precision basic data for subsequent focal mechanism calculations and seismogenic tectonic studies.

[0079] In summary, this embodiment successfully extracted 22 abnormal signals from the original monitoring data containing noise by constructing a closed-loop system of "physical screening - statistical cleaning - dynamic judgment". It also verified the physical and statistical dual reset reliability of the remaining data and achieved accurate assessment and location of the 3.x magnitude earthquake in Northwest China.

[0080] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for assessing the reliability of seismic events, characterized in that, include: Acquire earthquake event data, including spatiotemporal parameters, seismic phases, and station information; Calculate epicentral distance and observation travel time based on earthquake event data; A near-site influence range model is constructed based on epicentral distance and near-site determination threshold. Seismic phases are divided into near-field and far-field phases according to epicentral distance and observation travel time. The near-field and far-field phases are filtered for consistency by matching the corresponding travel time window constraints. After removing noise, an effective set of phases and travel time residuals are generated. A cascaded screening process is performed on the effective seismic phase set, and outlier seismic phases are removed through iterative mean deviation analysis. A regression model of travel time residuals and epicentral distance is established to eliminate abnormal seismic phases that deviate from the regression trend line. A station spatial coverage topology is constructed to detect the gradient change of the epicentral distance reception ratio, and far-end seismic phases with a sudden drop in ratio are removed to generate a high-confidence seismic phase subset. Calculate the spatial distribution score of the high-confidence seismic phase subset; construct a dynamic threshold model using the number of seismic phases as an inverse variable, combined with the upper limit of the maximum threshold and the lower limit of the minimum threshold, to generate the pass threshold; determine the reliability of the earthquake event based on the spatial distribution score and the pass threshold.

2. The earthquake event reliability assessment method according to claim 1, characterized in that, By querying the business database and parsing earthquake reports and monitoring waveform files, event information, station information, and seismic phase information are obtained and integrated to construct a standardized dataset containing seismic phase objects. For each seismic phase object in the standardized dataset, the following parameters are calculated: Extract the longitude and latitude of the station to which the seismic phase belongs, as well as the longitude and latitude of the earthquake event, and calculate the surface arc length between the station location and the epicenter of the event as the epicentral distance; extract the actual arrival time of the seismic phase and the time of occurrence of the earthquake event, and calculate the time difference between the actual arrival time and the time of occurrence as the observation travel time.

3. The earthquake event reliability assessment method according to claim 1, characterized in that, The specific steps for consistency filtering using the near-site influence range model include: setting a near-site determination threshold; for each seismic phase in the dataset, if the epicentral distance is greater than the near-site determination threshold, it is defined as a far-field seismic phase, and the observed travel time meets the following numerical range constraints: the observed travel time value is greater than the ratio of the epicentral distance to the upper limit of the P-wave propagation velocity, and less than the ratio of the epicentral distance to the lower limit of the P-wave propagation velocity; if the epicentral distance is not greater than the near-site determination threshold and any condition within a predetermined time is met, it is defined as a near-field seismic phase, and a relaxed near-field time window is adopted to allow the existence of direct waves and first waves; seismic phases falling outside the travel time window are removed; if the number of remaining seismic phases after removal is less than the number of associated stations, the event is marked as a false event and the evaluation is terminated.

4. The earthquake event reliability assessment method according to claim 1, characterized in that, The calculation of the travel time residuals of the remaining seismic phases includes: Theoretical travel time acquisition: For each seismic phase, the theoretical arrival time is queried from a pre-set travel time table based on the epicentral distance and focal depth; if it exceeds the coverage of the travel time table, the TauPyModel calculation module is called, the Earth velocity model is selected, and the corresponding theoretical arrival time is dynamically calculated using the ray tracing method. Phase type correction: If multiple theoretical phases have the same theoretical travel time at the same epicentral distance, the theoretical phase with the highest physical propagation probability at that distance is selected as the benchmark phase according to the phase dominance interval of the epicentral distance. Residual calculation: The difference between the observed travel time and the theoretical arrival time is the travel time residual of a single seismic phase; Overall correction: The arithmetic mean of the absolute values ​​of the residuals of all individual seismic phases is used as the overall residual benchmark for the event.

5. The earthquake event reliability assessment method according to claim 1, characterized in that, The specific method for removing outlier phases through iterative mean deviation analysis is as follows: Initialize the phase set and calculate the mean of the residuals of all phases in the set; traverse each phase in the set and calculate the mean of the residuals of the remaining set after assuming the removal of outlier phases; calculate the percentage change in mean caused by the removal of phases, which is the ratio of the absolute value of the difference between the original mean and the new mean to the absolute value of the original mean; determine whether the following two removal conditions are met simultaneously: Condition 1: The percentage change in mean is greater than the preset mean deviation threshold; Condition 2: The absolute value of the travel time residual of the phase is greater than the preset absolute residual threshold; If met, the phase is removed, and the next round of iteration calculation is performed until there are no more phases in the set that meet the removal conditions.

6. The method for assessing the reliability of earthquake events according to claim 1, characterized in that, The specific steps for establishing a regression model between travel time residuals and epicentral distance include: for the remaining seismic phases after mean deviation analysis, constructing a linear regression equation with epicentral distance as the independent variable and residuals as the dependent variable, and determining the slope and intercept of the regression line; calculating the corrected residual for each seismic phase, where the corrected residual is the difference between the actual observed residual and the predicted residual calculated based on the linear regression equation; calculating the mean of all corrected residuals and the average of the absolute deviations of all corrected residuals relative to the mean; and defining the elimination criteria: first, checking the number of seismic phase samples participating in the regression; if the number of samples is less than a preset statistical significance threshold, then only absolute threshold discrimination is performed: if the absolute value of the corrected residual of a seismic phase is greater than a preset absolute tolerance upper limit, then the seismic phase is eliminated; if the number of samples is not less than the statistical significance threshold, then statistical discrimination is performed: if the deviation of the corrected residual of a seismic phase from the mean is greater than a specified multiple of the preset standard deviation, and the absolute value of the original residual of the seismic phase is greater than a preset absolute residual threshold, then the seismic phase is eliminated.

7. The method for assessing the reliability of earthquake events according to claim 1, characterized in that, The process of constructing the spatial coverage topology of the stations and detecting gradient changes includes: constructing an epicentral distance sequence, extracting the epicentral distance lists of all online stations and all actually detected seismic phases, and sorting them in ascending order; calculating the reception ratio sequence, for each epicentral distance point in the list of actually detected seismic phases, defining a forward-covering sliding window, and calculating the local reception ratio within the window, which is the ratio of the number of actually detected seismic phases within the window distance range to the total number of online stations within the range; Perform abrupt drop detection: If there is a distance point in the received ratio sequence where the local received ratio corresponding to the distance point is less than a preset local density threshold, and the distance interval between the distance point and the preceding phase is greater than the theoretical attenuation interval of the regional seismic wave, then all phases with an epicentral distance not less than the distance point are removed; perform joint abrupt change detection, calculate the difference in received ratios and the difference in epicentral distance intervals between two adjacent distance points; if the decrease in the received ratio of adjacent phases exceeds a preset change threshold, and the distance interval between adjacent epicentral distances is greater than the average interval of all stations in the network, then a spatial fault is determined to have occurred, and all distant phases after the fault are removed.

8. The method for assessing the reliability of earthquake events according to claim 1, characterized in that, The specific expression of the dynamic threshold model and its judgment logic is as follows: The spatial distribution score is the arithmetic mean of the local reception ratios of all earthquake phases in the seismic phase set; the calculation logic of the dynamic threshold model is as follows: calculate the difference between the maximum upper threshold and the minimum lower threshold, multiply the difference by the convergence coefficient, divide by the total number of seismic phases associated with the current event, and add the quotient to the minimum lower threshold to obtain the passing threshold; wherein, if the obtained quotient is greater than the difference between the maximum upper threshold and the minimum lower threshold, the difference between the maximum upper threshold and the minimum lower threshold is taken as the quotient used for calculation; the convergence coefficient is used to adjust... The rate at which the threshold changes with the number of seismic phases; the upper limit of the maximum threshold is positively correlated with the sparseness of stations in the monitoring area; the lower limit of the minimum threshold is negatively correlated with the density of stations in the monitoring area; a preliminary reliable judgment condition is that the spatial distribution score is not less than the passing threshold; the arithmetic square root of the ratio of the coverage area of ​​the seismic event network to the number of stations is calculated to obtain the average station spacing; the value corresponding to the average station spacing is extracted based on the preset parameter mapping relationship, wherein the average station spacing value is positively correlated with the extracted values ​​of the upper limit of the maximum threshold and the lower limit of the minimum threshold.

9. The method for assessing the reliability of earthquake events according to claim 1, characterized in that, The specific implementation method for determining the reliability of an earthquake event based on the spatial distribution score and the pass threshold is as follows: When the spatial distribution score is less than the pass threshold, a secondary verification process is initiated; the root mean square error of the spatial distribution score sequence is calculated, where the root mean square error is the square root of the arithmetic mean of the sum of the squares of the differences between the reception proportions of each seismic phase and the overall spatial distribution score of the event; the verification must pass simultaneously to meet the following two conditions: Condition 1: The root mean square error is less than half of the difference between the upper limit of the maximum threshold and the lower limit of the minimum threshold; Condition 2: The spatial distribution score is not less than the relaxed threshold, where the relaxed threshold is defined as the pass threshold minus the lower limit of the minimum threshold and the convergence coefficient; If both conditions are met, the reliability status of the event will be corrected to "reliable"; otherwise, the "false event" determination will be maintained.

Images