Method and system for predicting earthquake in lithosphere magnetic field vector weakly-varying region

Abstract

The invention discloses a lithosphere magnetic field vector weakly-varying region earthquake prediction method and system, and belongs to the field of earthquake geomagnetic field precursor anomaly detection. Lithosphere magnetic field data is acquired, horizontal and vertical vectors are extracted, a vector mean value is calculated according to grids, a Molchan chart method is adopted to evaluate prediction efficiency and reversely extract a threshold value, a low-value region is screened to generate weakly-varying region data, and the weakly-varying region data is extracted. According to the method, the limitation of single component analysis is overcome, objective evaluation of prediction efficiency and automatic extraction of a threshold value are realized, the corresponding relation between the fault property and the vector feature is established, and the accuracy of earthquake occurrence site prediction is remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of earthquake geomagnetic precursor anomaly detection, specifically to an earthquake prediction method and system based on weakly variable zones of the lithosphere magnetic field vector. Background Technology

[0002] Earthquakes are natural phenomena caused by crustal activity. Disastrous earthquakes can cause secondary disasters such as building collapses, ground subsidence, and landslides, posing a significant threat to people's property and lives. With the advancement of digital seismic network construction and the standardization of seismic electromagnetic equipment, geomagnetic anomaly signals are considered reliable earthquake precursors, especially before major earthquakes, as they contain potential earthquake precursor information.

[0003] Chinese invention patent CN 116449412 B discloses a method for analyzing the spatiotemporal variation of the geomagnetic field for short-term earthquake monitoring. This method divides geomagnetic stations within a given latitude and longitude range into multiple rectangular grids. It preprocesses the multi-component data from each station's detection day and the previous 60 days into multiple standard time series. Then, it calculates the daily absolute average error of each component's predicted sequence and uses a sliding interquartile range (interquartile range) detection algorithm to detect anomalies in the difference sequence, obtaining the daily anomaly values ​​for each component. Next, it removes anomalies caused by geomagnetic storms based on the daily geomagnetic disturbance index. Finally, it generates a daily anomaly index matrix for each station based on the grid size. The final regional daily anomaly index matrix is ​​obtained by summing the anomaly index matrices from all stations' detection day and the previous 29 days, thus determining the latitude and longitude range of the anomalies.

[0004] However, the aforementioned existing technologies have the following shortcomings: First, this method only uses geomagnetic data for anomaly detection, failing to fully explore the joint characteristic information of horizontal and vertical vectors, and making it difficult to distinguish the differentiated magnetic field response characteristics corresponding to different fault types; Second, the anomaly detection threshold of this method is artificially set based on statistical methods, lacking an objective quantitative evaluation mechanism for prediction effectiveness, and unable to dynamically optimize the threshold parameters based on historical earthquake case retrospective results; Third, the anomaly area output by this method is relatively large, only providing a geographical reference range for possible earthquakes, making it difficult to accurately pinpoint the earthquake origin; Fourth, this method does not establish a correspondence between fault properties and magnetic field vector characteristics, and cannot make targeted epicenter location predictions based on the different characteristics of reverse faults or strike-slip faults. Therefore, existing earthquake geomagnetic anomaly detection and analysis methods still need improvement in terms of prediction accuracy and objective quantification. Summary of the Invention

[0005] To address the aforementioned problems in existing technologies, this invention aims to provide a method and system for earthquake prediction in areas with weak variations in lithospheric magnetic field vectors. By jointly analyzing horizontal and vertical vectors, combining the Molchan chart method to establish an objective evaluation mechanism for prediction effectiveness, and introducing fault property discrimination rules, accurate prediction of earthquake occurrence locations can be achieved.

[0006] The first aspect of this invention provides a method for earthquake prediction in areas with weakly variable lithospheric magnetic field vectors, comprising the following steps: acquiring lithospheric magnetic field data from multiple mobile geomagnetic measuring points within a preset time period in a target study area, and extracting horizontal and vertical vector data from each measuring point; dividing the target study area into multiple grid points according to a preset grid spacing, and calculating the mean horizontal vector and the mean comprehensive vector of each grid point; constructing a spatiotemporal alarm model using the Molchan chart method, calculating the missed alarm rate and spatiotemporal occupancy rate, generating a prediction performance evaluation curve, and extracting the mean vector threshold in reverse; filtering grid points in low-value areas based on the mean vector threshold to generate spatial data of the weakly variable area; determining the fault nature based on the directional characteristics of the vertical vector on both sides of the fault, and generating a fault nature discrimination result; and combining the spatial data of the weakly variable area and the fault nature discrimination result to determine the spatial range of the potential seismic source area and output the prediction result.

[0007] The second aspect of the present invention provides an earthquake prediction system for weakly variable lithospheric magnetic field vector zones, including a data acquisition module, a vector mean calculation module, an effectiveness evaluation module, a weakly variable zone extraction module, a fault discrimination module, and a source zone output module. Each module is connected in sequence to form a data processing link, realizing a complete processing flow from raw lithospheric magnetic field data to prediction results of potential source zones.

[0008] Compared with existing technologies, this invention has the following advantages: First, it integrates the joint analysis of horizontal and vertical vectors, making full use of all vector information of the lithosphere magnetic field, overcoming the limitations of single-component analysis, and improving the reliability of anomaly identification; Second, it introduces the Molchan chart method for objective evaluation of prediction effectiveness, verifies prediction capabilities through historical earthquake examples, and extracts the optimal threshold based on spatiotemporal occupancy, realizing the transformation from subjective experience judgment to objective numerical prediction; Third, it establishes the correspondence between fault properties and vector characteristics. When an earthquake occurs on a reverse fault, the vertical vector shows a boundary reversal characteristic, while when an earthquake occurs on a strike-slip fault, the horizontal vector weakens significantly and the vertical vector spreads in the same direction, realizing a differentiated epicenter location strategy; Fourth, it adopts the optimal spatiotemporal occupancy threshold of predicting 60% of earthquakes with 20% of the area, significantly narrowing the prediction area and improving the accuracy of earthquake location prediction. Attached Figure Description

[0009] Figure 1 This is a flowchart of the earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to the present invention.

[0010] Figure 2 This is an architectural diagram of an earthquake prediction system for weakly variable lithosphere magnetic field vector zones according to the present invention.

[0011] Figure 3 This is a spatial distribution map showing the annual changes in the horizontal vector of the lithosphere magnetic field in the North-South Seismic Belt from 2010 to 2024.

[0012] Figure 4 This is a spatial distribution map showing the annual changes in the vertical vector of the lithosphere magnetic field in the North-South Seismic Belt from 2010 to 2024.

[0013] Figure 5 This is a graph showing the effectiveness assessment results of the Molchan chart method for the H value of the lithosphere magnetic field in the North-South Seismic Belt.

[0014] Figure 6 This is a graph showing the effectiveness assessment results of the Molchan chart method for the F-values ​​of the lithosphere magnetic field in the North-South Seismic Belt.

[0015] Figure 7 It is a spatial distribution map of the predicted region based on threshold extraction using the Molchan chart method.

[0016] Figure 8 This is a three-stage evolution distribution map of the horizontal vector of changes in the lithosphere magnetic field before the Lushan 7.0 magnitude earthquake.

[0017] Figure 9 This is a spatial distribution map of potential seismic source areas for long-term location prediction using multiple observation methods in Yunnan Province. Detailed Implementation

[0018] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. Those skilled in the art should understand that the following embodiments are only used to illustrate the technical solution of the present invention and are not intended to limit the scope of protection of the present invention.

[0019] See Figure 1 This invention provides a method for earthquake prediction in areas with weakly variable lithospheric magnetic field vectors. The core idea of ​​this method is to jointly analyze the horizontal and vertical vectors of the lithospheric magnetic field, use the vector mean method to numerically characterize the weakly variable areas, establish an objective evaluation mechanism for prediction effectiveness using the Molchan chart method, and extract the optimal threshold in reverse. Simultaneously, fault property discrimination rules are introduced to differentiate between reverse faults and strike-slip faults, ultimately determining the spatial extent of potential seismic source areas. In one embodiment of this invention, the method specifically includes the following steps: Step S1: Acquisition and vector extraction of lithosphere magnetic field data.

[0020] In this step, it is first necessary to acquire lithospheric magnetic field data from multiple mobile geomagnetic measuring points within the target study area over a preset time period. In one embodiment of the present invention, the target study area is taken as the North-South Seismic Belt, which covers a region with latitudes from 21° to 37° and longitudes from 97° to 107°. This region is one of the areas in China with the most frequent strong earthquake activity and has a good sample of earthquake cases. The preset time period is preferably a variation cycle between two consecutive years, such as June 2010 to June 2011 as a complete observation cycle. This setting can eliminate the influence of seasonal variations and highlight the interannual scale tectonic magnetic field variation characteristics.

[0021] Acquisition of lithospheric magnetic field data relies on a mobile geomagnetic observation network, which consists of multiple mobile geomagnetic measuring points distributed throughout the target study area. Each measuring point is equipped with a high-precision magnetometer capable of measuring the three components of the geomagnetic field. Extracting the horizontal and vertical vector data from each measuring point from the acquired lithospheric magnetic field data is a crucial step in this process. In one embodiment of the invention, the horizontal vector data includes the northward component X and the eastward component Y of the magnetic field, which together constitute the vector components in the horizontal plane; the vertical vector data is the vertical component Z of the magnetic field, characterizing the variation of the magnetic field in the vertical direction. During the data preprocessing stage, the raw observation data needs to undergo processes such as geomagnetic storm removal, diurnal variation correction, and long-term variation separation to obtain a clean lithospheric magnetic field anomaly signal. Preferably, when the geomagnetic disturbance index Dst is less than -40 nT, the observation data for that day is marked as affected by a geomagnetic storm and removed in subsequent analysis.

[0022] In one embodiment of this invention, the temporal resolution of the data is the annual variation, that is, the difference in magnetic field components at the same measuring point in two adjacent years is calculated as the annual variation vector of that measuring point. Taking the horizontal vector as an example, the change in the horizontal vector of a measuring point in year t and year t+1 can be expressed as the vector difference between the observed value of that year and the observed value of the previous year. This differential processing method can effectively extract magnetic field variation information related to tectonic activity, while eliminating the influence of the long-term variation trend of the geomagnetic field. Through the above steps, the horizontal and vertical vector data of all measuring points in the target study area in each observation period can be obtained, providing a data foundation for subsequent gridded calculations.

[0023] Step S2: Mesh generation and vector mean calculation.

[0024] After data acquisition, the target study area needs to be divided into multiple grid points according to a preset grid spacing, and the vector mean of each grid point needs to be calculated. The purpose of gridding is to convert discretely distributed measurement point data into regularly distributed spatial field data, which facilitates subsequent spatial analysis and threshold extraction. In one embodiment of the present invention, the preset grid spacing is preferably 0.1 degrees, that is, one grid point is set every 0.1 degrees along the longitude and latitude directions. This grid density can ensure spatial resolution and find a sufficient number of measurement points around most grid points for statistical calculation.

[0025] For each grid point, the number of measuring points within a preset radius around that grid point is counted. Subsequent vector mean calculations are only performed when the number of measuring points meets a preset minimum number of measuring points. In one embodiment of this invention, the preset radius is preferably 100 km, and the preset minimum number of measuring points is preferably 3. This setting is based on the following: a 100 km radius can cover the regional-scale magnetic field anomaly characteristics while avoiding signal smoothing effects caused by an excessively large area; the requirement of a minimum of 3 measuring points ensures the reliability of the statistical results and avoids random errors caused by too few measuring points. For areas with sparsely distributed measuring points, when the number of measuring points within a 100 km radius around a grid point is less than 3, that grid point is not included in subsequent calculations and is displayed as a blank area in the result graph.

[0026] The calculation of the mean horizontal vector is one of the core aspects of this step. In one embodiment of the present invention, the method for calculating the mean horizontal vector H is as follows: First, collect the horizontal vector change data of all measuring points within a 100km radius around a certain grid point; then, sum the horizontal vectors of these measuring points to obtain the resultant vector; finally, divide the magnitude of the resultant vector by the number of measuring points to obtain the mean horizontal vector of the grid point. Its mathematical expression is: , Where H is the mean horizontal vector, in nT, which represents the average intensity of the change in the horizontal vector of the measuring points around the grid point, and the value range is usually 0 to 20 nT; The horizontal vector change of the m-th measuring point is a two-dimensional vector composed of the changes in the X and Y components of the measuring point, with the unit being nT; m is the number of measuring points within 100km around the grid point that meet the conditions, with the unit being the number of points, and the value is not less than 3; the vertical line symbol in the numerator indicates the operation of taking the vector magnitude value, that is, calculating the magnitude of the vector summation result.

[0027] The physical meaning of the mean horizontal vector H is as follows: when the horizontal vectors at all measuring points within a certain area change in the same direction and with similar magnitudes, the sum of the vectors is larger, and the resulting H value after dividing by the number of measuring points is also larger, indicating that the area is in a state of uniform change. Conversely, when the horizontal vectors at all measuring points within a certain area change in a scattered and mutually canceling direction, the sum of the vectors is smaller, and the corresponding H value is also smaller, indicating that the area is in a state of weak change, i.e., the so-called "weakly changing zone". Retrospective analysis of earthquake cases shows that earthquake epicenters are often located near or inside the boundary of weakly changing zones, and this spatial correspondence constitutes the physical basis of the prediction method of this invention.

[0028] The calculation method for the comprehensive vector mean is similar to that for the horizontal vector mean, the difference being that the input data includes not only horizontal vectors but also vertical vectors. In one embodiment of the present invention, the calculation method for the comprehensive vector mean F is as follows: , Where F is the mean value of the composite vector, in nT, which represents the average intensity of the combined changes of the horizontal and vertical vectors of the measuring points around the grid point, and its value range is usually from 0 to 25 nT; The comprehensive vector change of the m-th measuring point is a three-dimensional vector composed of the changes in the X, Y, and Z components of that measuring point, with units of nT; the meanings of the other symbols are the same as those in the formula for the mean value of the horizontal vector.

[0029] Compared to the horizontal vector mean H, the composite vector mean F incorporates additional information from the vertical vector, thus providing a more comprehensive reflection of the three-dimensional variations in the lithospheric magnetic field. Retrospective analysis of earthquake cases shows that, within the testing period prior to 2019, the predictive efficacy of the composite vector mean F was superior to using the horizontal vector mean H alone, indicating that the vertical vector contains important precursory information. In one embodiment of this invention, both the H and F vector means are calculated simultaneously, and subsequent efficacy evaluations and threshold extractions are performed separately to obtain more reliable prediction results.

[0030] Through the above gridded calculation process, the discretely distributed measurement point data can be converted into H-value and F-value matrices covering the entire target study area. For example... Figure 3 As shown, the spatial distribution of the mean horizontal vector values ​​of the North-South Seismic Belt from 2010 to 2024 exhibits a clear spatiotemporal evolution characteristic, with the location and extent of low-value areas changing from year to year. Figure 4 As shown, the spatial distribution of vertical vectors during the same period also exhibits characteristic changes related to seismic activity. By marking the epicenters of earthquakes of magnitude 5 or above that will occur within the next year on the spatial distribution map of the vector mean for that year, the spatial correspondence between the epicenter and the boundary of the low-value area can be clearly seen.

[0031] Step S3: Predictive performance evaluation and threshold extraction based on Molchan chart method After calculating the vector mean, the effectiveness of the prediction method needs to be objectively evaluated, and the optimal vector mean threshold should be extracted based on the evaluation results. This invention uses the Molchan chart method to construct a spatiotemporal alarm model. This method is a widely recognized standard for evaluating prediction effectiveness in the international seismological community, and it can quantitatively describe the trade-off between the prediction method's spatiotemporal occupancy rate and the missed detection rate.

[0032] The core idea of ​​the Molchan chart method is to plot prediction performance curves under different threshold conditions, with the spatiotemporal occupancy rate τ as the horizontal axis and the false negative rate v as the vertical axis. In one embodiment of this invention, the formula for calculating the false negative rate v is: , Where, v is the false negative rate, dimensionless, ranging from 0 to 1, representing the proportion of earthquakes that were predicted to occur but actually did; N is the total number of actual earthquakes occurring in the target study area during the test period, measured in times, and includes earthquakes of magnitude 5.0 and above; h is the number of earthquakes that were predicted to occur, measured in times, i.e., the number of earthquakes whose epicenters fell within the prediction area; and Nh is the number of earthquakes that were missed, i.e., the number of earthquakes whose epicenters did not fall within the prediction area. The smaller the false negative rate v, the higher the hit rate of the prediction method and the better the prediction performance.

[0033] The formula for calculating the spacetime occupancy rate τ is: , Where τ is the spatiotemporal occupancy rate, a dimensionless value ranging from 0 to 1, representing the proportion of the predicted region to the total area of ​​the study area; A is the area of ​​the anomaly zone, in km², which is the spatial range covered by grid points whose vector mean is lower than the current threshold; R is the total area of ​​the study area, in km², which is the entire spatial range of the target study area. The smaller the spatiotemporal occupancy rate τ, the more concentrated the predicted region and the higher the spatial resolution.

[0034] In the Molchan diagram, the diagonal line τ+v=1 represents the prediction result of random guessing. Any prediction curve above this diagonal line indicates that the prediction effect is better than random guessing, that is, the prediction method has actual predictive ability. In one embodiment of the present invention, the formula for calculating the area skill score ASS is: , Wherein, ASS is the area skill score, dimensionless, ranging from 0 to 1; τ is the spacetime occupancy rate, integrated from 0 to 1 as an integral variable; v(τ) is the false negative rate function value corresponding to a spacetime occupancy rate of τ; (1-τ-v(τ)) represents the vertical distance between the prediction curve and the random line. A larger area skill score (ASS) indicates that the prediction curve is further away from the random line, and the prediction performance is better. When ASS is greater than 0.5, the prediction performance is considered to pass the test, and the prediction results at this threshold have practical application value; when ASS is less than or equal to 0.5, the prediction performance is considered to fail the test, and the prediction results at this threshold are no better than random guessing.

[0035] In one embodiment of the present invention, a 5-year period is used as a complete prediction and testing cycle, and a sliding test is performed with a 1-year step size, which can obtain multiple sets of performance evaluation results for different time periods. For example... Figure 5-6 As shown, the Molchan chart method performance evaluation results for the H and F values ​​of the lithospheric magnetic field in the North-South seismic belt during the period of 2011-2023 indicate that: for the H value, the prediction curves for the three test periods of 2013-2017, 2014-2018, and 2015-2019 are mostly above the random line, indicating that the prediction performance of using the low-value area of ​​the horizontal vector for earthquakes of magnitude 5 and above is better than random guessing during the period of 2013-2019; for the F value, the prediction curves for the four test periods of 2011-2015, 2012-2016, 2013-2017, and 2014-2018 are all above the random line, and about half of the prediction curves for 2015-2019 are above the random line. Before 2019, the overall prediction performance of the comprehensive vector mean F is better than that of using the horizontal vector mean H alone, because the F value usually has a higher detection rate under conditions of lower spatiotemporal occupancy.

[0036] Based on the effectiveness evaluation results of the Molchan chart method, this invention employs a reverse extraction method to determine the optimal vector mean threshold. In one embodiment of this invention, the preset spatiotemporal occupancy threshold is preferably 0.2, i.e., τ=0.2 is selected as the benchmark point for threshold extraction. When the spatiotemporal occupancy τ equals 0.2, the corresponding false negative rate v is read from the prediction effectiveness curve, and the vector mean threshold that makes the low-value area proportion exactly 20% is calculated in reverse. The advantage of this threshold extraction method is that: using 20% ​​of the area can predict about 60% of earthquakes of magnitude 5 and above, achieving the optimal balance between prediction accuracy and coverage; the determination of the threshold is entirely based on the statistical analysis results of historical earthquake cases, avoiding the arbitrariness brought about by subjective human setting.

[0037] like Figure 7As shown, the spatial distribution of the prediction area, plotted according to the threshold for each year, indicates that the epicenters of the two largest earthquakes in Sichuan Province during the period of 2010-2023—the magnitude 7.0 earthquake in Lushan and the magnitude 7.0 earthquake in Jiuzhaigou—both fell within or near the boundaries of the prediction area. Specifically, the epicenter of the magnitude 7.0 earthquake in Lushan fell within the prediction area for 2011-2012; there was a small local prediction area near the epicenter of the magnitude 7.0 earthquake in Jiuzhaigou, and the epicenter location of Jiuzhaigou was also covered in the prediction area of ​​the previous year. In addition, although the epicenters of strong earthquakes such as the magnitude 6.6 earthquakes in Minxian and Zhangxian, the magnitude 6.4 earthquakes in Kangding, the magnitude 6.5 earthquakes in Ludian, the magnitude 6.4 earthquakes in Yangbi, and the magnitude 6.1 earthquake in Lushan had low coverage in the current prediction area, their coverage was significantly higher in the prediction area of ​​the previous year. This suggests that extending the prediction time window to 18 months has a higher success rate.

[0038] Step S4: Generation of spatial data for weakly variable regions and identification of spatiotemporal evolution characteristics.

[0039] After threshold extraction, it is necessary to filter low-value grid points within the target study area based on the vector mean threshold to generate spatial data of weakly variable regions. In one embodiment of this invention, a weakly variable region is defined as follows: in the annually changing regional lithosphere magnetic field, the area of ​​grid points where the horizontal vector mean H or the comprehensive vector mean F is lower than the corresponding threshold is a weakly variable region. Numerically, the vector mean value within a weakly variable region is in a relatively small range; in terms of image features, the vector arrows within a weakly variable region are scattered, inconsistent, and have a small magnitude close to zero. Compared with the surrounding stable and uniformly changing regions, weakly variable regions have clearly identifiable differences.

[0040] The generation process of spatial data for weakly variable regions is as follows: First, read the H-value matrix or F-value matrix for the current year; then, traverse each grid point and determine whether its vector mean is lower than the corresponding threshold; finally, connect all grid points that meet the conditions to form a closed region, thus forming the spatial boundary of the weakly variable region. In one embodiment of the present invention, the spatial data of the weakly variable region is stored in vector format, containing the latitude and longitude coordinate sequence of the boundary vertices, which facilitates subsequent spatial overlay analysis and visualization.

[0041] In addition to extracting weakly variable regions in a single year, this invention also focuses on the spatiotemporal evolution characteristics of these regions. For example... Figure 8 As shown, taking the Lushan 7.0 magnitude earthquake as an example, the pre-earthquake weak change zone went through three typical evolution stages: The first stage was 2 to 3 years before the earthquake, when a large area of ​​weak change zone appeared in the horizontal vector of the lithosphere magnetic field change, covering a wide range with unclear boundaries; the second stage was 2 to 7 months before the earthquake, when the weak change zone shrank to a local small area, while the vector direction of the surrounding area showed an overall convergence characteristic, forming a clear directional field pattern; the third stage was 1 to 1 month before the earthquake, when the weak change zone evolved into a large area again, and the vector direction of the area near the epicenter showed reverse changes or weak changes.

[0042] In one embodiment of the present invention, the shrinkage ratio of the weakly variable zone area is used as a criterion for entering the short-term prediction stage. Specifically, when the area of ​​the weakly variable zone in the second stage is no more than one-third of the area in the first stage, it is determined that the earthquake-inducing process in the source area has entered a critical stage, requiring enhanced monitoring and tracking. This determination method based on area evolution characteristics has a clear physical basis: a large-area weakly variable zone corresponds to the early stage of stress accumulation, a localized small-scale weakly variable zone corresponds to the earthquake-inducing stage of stress concentration, and a further expanding weakly variable zone corresponds to the stress release preparation stage before an earthquake.

[0043] Retrospective analysis of earthquakes of magnitude 6 or higher in the North-South Seismic Belt since 2011 verified the universality of the aforementioned three-stage evolutionary characteristics. Multiple earthquakes, including the Lushan 7.0 earthquake, the Minxian-Zhangxian 6.6 earthquake, and the Yangbi 6.4 earthquake, exhibited similar evolutionary patterns, indicating that this characteristic evolution is prevalent before strong earthquakes and possesses operability as a predictive indicator.

[0044] Step S5: Fault nature identification and differential epicenter location.

[0045] After extracting the weakly variable zone, it is necessary to further analyze the directional characteristics of the vertical vector on both sides of the fault, determine the fault nature, and generate differentiated epicenter location strategies. A key innovation of this invention lies in discovering the close relationship between the regional characteristics of the vertical vector and the fault nature; this correspondence provides crucial evidence for accurately pinpointing the epicenter location.

[0046] In one embodiment of the present invention, the fault nature discrimination rule is as follows: when the vertical vectors exhibit opposite directions on both sides of the fault, it is determined to be a reverse fault-related area. Taking the Lushan 7.0 magnitude earthquake as an example, this earthquake occurred in the southern segment of the Longmenshan fault and is a typical reverse earthquake. The spatial distribution of vertical vectors one year before the earthquake shows that, with the Longmenshan fault as the boundary, the vertical vector arrows point upwards on the east side and downwards on the west side, exhibiting obvious reverse characteristics on both sides of the fault. The physical mechanism of this reverse characteristic is closely related to the kinematic characteristics of the reverse fault: reverse fault activity leads to relative uplift of the hanging wall and relative subsidence of the footwall, causing differential changes in the magnetization intensity of the lithosphere, which in turn results in the observed reverse distribution of vertical vectors on the surface.

[0047] When vertical vectors exhibit a consistent direction around a fault, while horizontal vectors show significant weakening, the area is identified as a strike-slip fault region. Taking the Jiuzhaigou 7.0 magnitude earthquake as an example, this earthquake occurred near the Minjiang Fault and is a left-lateral strike-slip earthquake. The pre-earthquake spatial distribution of vertical vectors shows that the vertical vector directions around the fault are basically consistent, with no obvious reverse characteristics; simultaneously, the horizontal vectors in this region show a significant weakening, with H values ​​at low levels. The physical mechanism of this combination of characteristics is related to the kinematic characteristics of strike-slip faults: strike-slip fault activity mainly manifests as horizontal shearing motion, with a smaller vertical component, thus the vertical vectors maintain a consistent direction; while the accumulation of stress in the horizontal direction leads to changes in the directional arrangement of magnetic minerals, manifested as a weakening of the horizontal vectors.

[0048] In one embodiment of the present invention, differentiated epicenter location strategies are employed for different fault types. For regions associated with reverse faults, the predicted epicenter is located at the intersection of the vertical vector reversal boundary and the boundary of the weakened zone. This location strategy can narrow the epicenter area from the entire weakened zone to a strip-shaped region near the boundary. For regions associated with strike-slip faults, the predicted epicenter is located within the weakened zone where the horizontal vector weakening is most significant. This location strategy emphasizes the degree of weakening rather than the boundary location.

[0049] The effectiveness of the above-mentioned discrimination rule was verified through joint analysis of fault properties and vector characteristics of earthquakes of magnitude 6 and above in the North-South Seismic Belt. As shown in Table 1, the vertical vectors of thrust earthquakes such as Lushan (magnitude 7.0), Menyuan (magnitude 6.4), Lushan (magnitude 6.1), and Jishishan (magnitude 6.2) all showed significant reverse characteristics along the seismogenic fault before the earthquake; while the vertical vectors of strike-slip earthquakes such as Luding (magnitude 6.8), Jiuzhaigou (magnitude 7.0), and Ludian (magnitude 6.5) showed a unidirectional distribution. This correspondence was repeated in multiple earthquake cases, demonstrating high reliability and operability.

[0050] Step S6: Comprehensive determination and prediction results of potential seismic source areas.

[0051] After extracting the weakly variable zone and determining the fault properties, it is necessary to comprehensively analyze the results from both aspects to determine the spatial range of the potential seismic source area and output the final prediction result. In one embodiment of the present invention, the determination of the potential seismic source area follows the guiding principle of "finding the source by the field", that is, starting from the spatial distribution characteristics of the regional magnetic field, gradually narrowing the prediction range, and finally locking the possible location of the earthquake.

[0052] The specific process for comprehensively determining potential seismic source areas is as follows: First, the spatial data of weakly variable areas extracted using the vector mean threshold is used as the background area, which accounts for about 20% of the total area of ​​the study area; then, the fault property discrimination results are superimposed on the weakly variable areas, and the range is further narrowed according to the different characteristics of reverse faults or strike-slip faults; finally, combined with the spatiotemporal evolution characteristics of the weakly variable areas, when the weakly variable areas shrink from a large area in the first stage to a small area in the second stage, the small area of ​​weakly variable areas is determined as potential seismic source areas.

[0053] In one embodiment of the present invention, the determination of the potential seismic source region can also be achieved by integrating the results of multiple observation methods for spatial overlay analysis. For example... Figure 9 As shown, taking Yunnan Province as an example, spatial b-value distribution data, GNSS strain rate field data, and low-value area data of the lithospheric magnetic field horizontal vector are overlaid to extract the overlapping anomaly areas jointly indicated by the three types of data. Specifically, spatial b-values ​​are used for medium- to long-term predictions from 2024 to 2028, identifying areas with low seismic activity and high stress accumulation; GNSS strain rate fields are used for continuous observations from 2022 to 2024, identifying areas with abnormal crustal deformation rates; and the 20% low-value area of ​​the lithospheric magnetic field horizontal vector is used to identify areas with weak changes. The three types of data share common delineation areas in regions such as the southern part of the Yuanmou-Luzhijiang Fault, the central part of the Qujiang-Jianshui Fault, and the central part of the Honghe Fault. These areas have been identified as potential medium- to long-term seismic source areas.

[0054] After identifying potential seismic source areas, it is necessary to further track the migration of minor earthquakes within the region as a triggering condition for short-term prediction. In one embodiment of the present invention, a minor earthquake monitoring window is set up in the potential seismic source area and its boundary region to statistically analyze the frequency and spatial distribution of minor earthquakes in each time period. When the frequency of minor earthquakes shows a sudden increase relative to the previous average level, the corresponding area is marked as a key area for short-term prediction, and an intensive monitoring and early warning program is initiated. This triggering mechanism based on minor earthquake activity can further narrow the time window based on the medium- and long-term prediction area, improving the timeliness of the prediction.

[0055] The prediction results include: spatial boundary coordinates of the potential seismic source area, predicted fault type, prediction time window, and prediction confidence level. In one embodiment of the present invention, the prediction time window is preferably within 12 to 18 months after the data acquisition is completed. Based on the aforementioned analysis, an 18-month prediction window has a higher success rate. The prediction confidence level is graded according to the effectiveness evaluation results of the Molchan chart method: when the ASS is greater than 0.7, it is considered high confidence; when the ASS is between 0.5 and 0.7, it is considered medium confidence; and when the ASS is less than 0.5, no prediction result is output.

[0056] See Figure 2This invention also provides an earthquake prediction system for areas with weak lithospheric magnetic field vector variation. This system, corresponding to the aforementioned method embodiments, can realize a complete processing flow from raw lithospheric magnetic field data to prediction results for potential seismic source areas. In one embodiment of this invention, the system includes the following modules: Data acquisition module 1 is used to acquire lithospheric magnetic field data from multiple mobile geomagnetic measuring points within the target study area, and to extract horizontal and vertical vector data. In one embodiment of the present invention, data acquisition module 1 includes a data interface unit and a preprocessing unit. The data interface unit is responsible for establishing a data connection with the mobile geomagnetic observation network and automatically acquiring observation data from each measuring point according to a preset time period. The preprocessing unit is responsible for quality control of the raw data, including magnetic storm removal, diurnal variation correction, outlier detection, etc., and outputting horizontal and vertical vector data that meet the requirements of subsequent analysis. The output of data acquisition module 1 serves as the input of vector mean calculation module 2, and the two modules are connected through a standardized data interface.

[0057] The vector mean calculation module 2 is connected to the data acquisition module 1 and is used to divide the target study area into grid points and calculate the horizontal vector mean H and the comprehensive vector mean F of each grid point. In one embodiment of the present invention, the vector mean calculation module 2 includes a grid division unit and a mean calculation unit. The grid division unit divides the study area into a regularly distributed array of grid points according to a preset grid spacing. The mean calculation unit traverses each grid point, counts the number of surrounding measurement points, and when the minimum number of measurement points is met, calls the aforementioned vector mean calculation formula and outputs the H value and F value of that grid point. The output of the vector mean calculation module 2 is an H value matrix and an F value matrix covering the entire study area, which serve as the input of the performance evaluation module 3.

[0058] The performance evaluation module 3 is connected to the vector mean calculation module 2. It is used to construct a spatiotemporal alarm model using the Molchan chart method, generate a prediction performance evaluation curve, and extract the vector mean threshold based on a preset spatiotemporal occupancy threshold. In one embodiment of the invention, the performance evaluation module 3 includes a performance calculation unit and a threshold extraction unit. The performance calculation unit reads the vector mean data and corresponding seismic catalog data for multiple test periods, calculates the missed detection rate v and spatiotemporal occupancy τ according to different threshold conditions, plots the τ-v curve, and calculates the area skill score (ASS). The threshold extraction unit locates the point corresponding to the preset spatiotemporal occupancy threshold on the τ-v curve, calculates the vector mean threshold in reverse, and determines whether the prediction performance for that test period passes the test. The output of the performance evaluation module 3 is the performance-verified vector mean threshold, which serves as the input to the weak variable area extraction module 4.

[0059] The weakly variable region extraction module 4 is connected to the performance evaluation module 3 and is used to filter low-value grid points based on a vector mean threshold to generate spatial data of the weakly variable region. In one embodiment of the present invention, the weakly variable region extraction module 4 includes a threshold filtering unit and a boundary generation unit. The threshold filtering unit traverses each grid point of the H-value matrix or F-value matrix and marks grid points with a vector mean lower than the threshold as weakly variable region grid points. The boundary generation unit connects adjacent weakly variable region grid points into closed regions and uses a boundary tracking algorithm to generate a spatial boundary coordinate sequence of the weakly variable region. In addition, the weakly variable region extraction module 4 also includes an evolution analysis unit, used to compare the area changes of the weakly variable region in different years and identify the three-stage spatiotemporal evolution characteristics. The output of the weakly variable region extraction module 4 is the spatial data of the weakly variable region, which also serves as the input to the fault discrimination module 5 and the source region output module 6.

[0060] The fault discrimination module 5 is connected to the weakly variable zone extraction module 4 and is used to determine the fault nature based on the directional characteristics of vertical vector data on both sides of the fault, generating a fault nature discrimination result. In one embodiment of the present invention, the fault discrimination module 5 includes a fault information database and a feature matching unit. The fault information database stores the spatial distribution information, geometric parameters, and historical activity records of known active faults within the target study area. The feature matching unit reads the vertical vector direction information from the spatial data of the weakly variable zone, spatially superimposes it with the fault location in the fault information database, and determines whether the vertical vector exhibits reverse or unidirectional characteristics on both sides of the fault, thereby outputting a discrimination result for reverse faults or strike-slip faults. The output of the fault discrimination module 5 is the fault nature discrimination result, which serves as the input to the source area output module 6.

[0061] The source area output module 6 is connected to the weakly variable area extraction module 4 and the fault discrimination module 5. It is used to determine the spatial extent of the potential source area by comprehensively analyzing the spatial data of the weakly variable area and the fault property discrimination results, and outputs the potential source area data. In one embodiment of the invention, the source area output module 6 includes a spatial overlay unit and a result output unit. The spatial overlay unit performs spatial overlay analysis on the weakly variable area boundary, the fault property discrimination results, and other optional observation data (such as spatial b-values, GNSS strain rate fields, etc.) to extract the anomalous overlapping areas jointly indicated by multi-source data. The result output unit formats and outputs information such as the spatial boundary coordinates of the potential source area, the predicted fault type, the prediction time window, and the prediction confidence level, supporting visualization and early warning information dissemination.

[0062] The six modules described above are sequentially connected to form a complete data processing chain. The output of data acquisition module 1 flows to vector mean calculation module 2; the output of vector mean calculation module 2 flows to performance evaluation module 3; the output of performance evaluation module 3 flows to weak variation zone extraction module 4; the output of weak variation zone extraction module 4 simultaneously flows to fault discrimination module 5 and source area output module 6; and the output of fault discrimination module 5 also flows to source area output module 6. This modular system architecture has good scalability and maintainability, facilitating functional upgrades and parameter optimization according to actual needs.

[0063] It should be noted that the above embodiments are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Equivalent substitutions or alternatives made based on the above technical solutions, such as using different grid spacings, different radius ranges, or different spatiotemporal occupancy thresholds, should all be included within the scope of the claims of the present invention.

[0064] Table 1. Correspondence between the properties and vector characteristics of faults with earthquake magnitude 6 or above in the North-South Seismic Belt

[0065]

[0066] As can be seen from Table 1, there are obvious differences between reverse fault type earthquakes and strike-slip fault type earthquakes in terms of vertical and horizontal vector characteristics. These differences provide a reliable criterion for the fault nature identification in this invention.

[0067] An evaluation of the predictive effectiveness using the Molchan chart method for the H and F values ​​of the lithospheric magnetic field in the North-South Seismic Belt from 2011 to 2023 reveals a significant annual variation in predictive effectiveness. During 2013-2018, both H and F values ​​demonstrated good predictive effectiveness, with an area skill score (ASS) greater than 0.5, indicating that using the low-value areas of the vector mean was more effective than random guessing for predicting earthquakes of magnitude 5 and above. However, after 2019, the predictive effectiveness of both vector means declined significantly, with the prediction curves even falling below the random line for some test periods.

[0068] The annual variation in predictive effectiveness may be related to the following factors: First, the methods and data processing procedures for mobile geomagnetic observations may have been adjusted around 2019, affecting data consistency; second, seismic activity itself exhibits periodic variations, and the seismic activity characteristics of the North-South seismic belt may have changed after 2019; third, the background field of the lithospheric magnetic field evolves over time, and changes in long-term trends may affect the accuracy of extracting annual variations. Therefore, in practical applications, it is recommended to use the threshold corresponding to the verification period that has passed the effectiveness test in the past five years for prediction, and to continuously update the threshold parameters as new data accumulates.

[0069] Compared with CN 116449412 B, this invention has significant technical advantages in the following aspects: First, regarding the anomaly detection mechanism, existing technologies use LSTM-AE neural networks for time series prediction, and detect anomalies through daily absolute average error and moving interquartile range method, with the threshold set manually based on statistical methods; while this invention uses the Molchan chart method for objective evaluation of prediction effectiveness, and determines the optimal threshold through historical earthquake case backtracking, realizing the transformation from subjective experience judgment to objective numerical prediction.

[0070] Second, in terms of vector information utilization, existing technologies calculate the outliers of each component separately and then simply sum them up, failing to fully utilize the directional information of the vector; while this invention calculates the vector mean through a vector summation method, preserving the directional characteristics of the vector, and can identify the degree of consistency of changes in each measuring point within the region, making the physical meaning of the weak variation zone clearer.

[0071] Third, in terms of spatial prediction accuracy, existing technologies output a large range of anomalous areas, which can only provide a geographical reference range where earthquakes may occur; while this invention significantly narrows the prediction area by using the optimal spatiotemporal occupancy threshold of predicting 60% of earthquakes with 20% of the area, and further achieves accurate location of the epicenter by combining fault property discrimination.

[0072] Fourth, in terms of fault type identification, existing technologies have not established a correspondence between fault properties and magnetic field anomalies; however, this invention proposes for the first time the difference between reverse faults and strike-slip faults in the vertical vector direction characteristics, providing differentiated epicenter location strategies for different types of earthquakes.

[0073] The earthquake prediction method and system for weakly variable lithosphere magnetic field vector zones provided by this invention can be applied to the following scenarios: First, medium- to long-term earthquake hazard analysis. By continuously monitoring the mean value of the regional lithosphere magnetic field vector, the spatial distribution and spatiotemporal evolution characteristics of weakly variable zones are identified, providing spatial references for areas where strong earthquakes may occur in the next 1-3 years.

[0074] Second, the delineation of key earthquake monitoring and defense zones. By combining multi-source observation data such as spatial b-values ​​and GNSS strain rate fields, potential seismic source areas jointly indicated by multiple methods are determined through spatial overlay analysis, providing a scientific basis for government departments to delineate key earthquake monitoring and defense zones.

[0075] Third, determining the triggering conditions for short-term forecasting. By tracking the three-stage spatiotemporal evolution characteristics of weakly variable zones and the migration of minor earthquakes, signals indicating the entry into the short-term forecasting phase are identified, providing decision support for the initiation of intensive monitoring and early warning procedures.

[0076] Fourth, the effectiveness evaluation of earthquake prediction research. The Molchan chart method used in this invention is a widely recognized standard for evaluating prediction effectiveness in the international seismological community, and can provide a reference framework for evaluating the effectiveness of other earthquake prediction methods.

[0077] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of the claims of the present invention.

Claims

1. A method for earthquake prediction in areas with weak variations in lithospheric magnetic field vectors, characterized in that: include: Acquire lithosphere magnetic field data from multiple mobile geomagnetic measuring points within a preset time period in the target study area, and extract horizontal and vertical vector data of each measuring point from the lithosphere magnetic field data. The horizontal vector data includes the northward and eastward magnetic field components, and the vertical vector data is the vertical magnetic field component. The target study area is divided into multiple grid points according to a preset grid spacing. For each grid point, the number of measurement points within a preset radius around the grid point is counted. When the number of measurement points meets the preset minimum number of measurement points, the horizontal vector mean and the comprehensive vector mean of the grid point are calculated respectively. Based on the mean horizontal vector and the mean comprehensive vector within multiple testing periods, a spatiotemporal alarm model is constructed using the Molchan chart method. The false alarm rate and spatiotemporal occupancy rate under different threshold conditions are calculated, and a prediction performance evaluation curve is generated. When the spatiotemporal occupancy rate is equal to the preset spatiotemporal occupancy threshold, the corresponding vector mean threshold is extracted in reverse. Based on the vector mean threshold, grid point regions within the target study area whose horizontal vector mean or comprehensive vector mean is lower than the vector mean threshold are selected to generate weak variable area spatial data. The spatial distribution information and geometric parameters of known active faults within the target study area are obtained. Based on the directional characteristics of the vertical vector data in the weakly variable spatial data on both sides of the fault, the nature of the fault is determined and the fault nature discrimination result is generated. When the vertical vectors on both sides of the fault show opposite directional characteristics, it is determined to be a reverse fault related area. When the vertical vectors around the fault show consistent directional characteristics and the horizontal vectors are significantly weakened, it is determined to be a strike-slip fault related area. By combining the spatial data of the weakly variable zone and the fault property discrimination results, the spatial range of the potential seismic source area is determined and the potential seismic source area data is output, thus completing the prediction of the earthquake occurrence location.

2. The earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to claim 1, characterized in that: The mean horizontal vector is obtained by summing the horizontal vector data of all surrounding measuring points, taking the modulus, and then dividing by the number of measuring points. The mean comprehensive vector is obtained by summing the horizontal and vertical vector data of all surrounding measuring points, taking the modulus, and then dividing by the number of measuring points. The preset spatiotemporal occupancy threshold is 0.2, corresponding to a 20% proportion of the weak variable area spatial data to the total area of ​​the target study area. The prediction performance evaluation curve has a prediction hit rate of no less than 60% when the spatiotemporal occupancy value is 0.

2.

3. The earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to claim 1, characterized in that: The preset grid spacing is 0.1 degrees, the preset radius range is 100km, and the preset minimum number of measuring points is 3.

4. The earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to claim 1, characterized in that: The preset time period is the change period between two consecutive years, and the inspection period is 5 years, with a 1-year step size for sliding inspection.

5. The earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to claim 1, characterized in that, The calculation process of the missed detection rate and the spatiotemporal occupancy rate includes: obtaining the missed detection rate by dividing the difference between the actual total number of earthquakes and the predicted number of earthquakes by the actual total number of earthquakes; obtaining the spatiotemporal occupancy rate by dividing the area of ​​the anomalous area by the total area of ​​the study area; and calculating the area skill score based on the prediction performance evaluation curve. When the area skill score is greater than 0.5, the prediction performance is determined to have passed the test.

6. The earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to claim 1, characterized in that, The spatial data of the weak change zone also includes the identification of spatiotemporal evolution characteristics, specifically including: identifying the characteristics of large-area weak change zones that appear 2 to 3 years before the earthquake; identifying the characteristics of local small-scale weak change zones that appear 2 to 7 months before the earthquake with the overall convergence of the outer vector directions; identifying the characteristics of weak change zones that evolve into large-scale zones again 1 to 1 month before the earthquake with the vector directions near the epicenter reversed or changed weakly; and using the shrinkage ratio of the weak change zone area not less than one-third of the area of ​​the previous stage as a criterion for entering the short-term prediction stage.

7. The earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to claim 1, characterized in that, Determining the spatial range of the potential seismic source area also includes: spatially overlaying the spatial data of the weakly variable area with the spatial b-value distribution data and GNSS strain rate field data, extracting the abnormal overlapping area jointly indicated by the three types of data, and tracking the migration path and frequency of small earthquakes within the abnormally overlapping area. When the frequency of small earthquakes shows a sudden increase relative to the previous average level, the corresponding area is marked as a key area for short-term prediction.

8. The earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to claim 1, characterized in that: For the inverse fault-related region, the predicted earthquake origin is located at the intersection of the vertical vector reversal boundary and the boundary of the weak variable zone spatial data; for the strike-slip fault-related region, the predicted earthquake origin is located within the weak variable zone spatial data where the horizontal vector weakening is most significant.

9. The earthquake prediction method for weakly variable lithosphere magnetic field vector zones according to claim 1, characterized in that: The acquisition window for the lithosphere magnetic field data is from June of each year to June of the following year, and the prediction window for the earthquake location is within 12 to 18 months after the data acquisition ends.

10. A lithospheric magnetic field vector weak variation zone earthquake prediction system, used to implement the lithospheric magnetic field vector weak variation zone earthquake prediction method according to any one of claims 1-9, characterized in that, include: The data acquisition module is used to acquire lithosphere magnetic field data from multiple mobile geomagnetic measuring points within the target study area, and to extract horizontal and vertical vector data. The vector mean calculation module, connected to the data acquisition module, is used to divide the target study area into grid points and calculate the horizontal vector mean and the comprehensive vector mean of each grid point. The performance evaluation module is connected to the vector mean calculation module. It is used to construct a spatiotemporal alarm model using the Molchan chart method, generate a predicted performance evaluation curve, and extract the vector mean threshold in reverse according to the preset spatiotemporal occupancy threshold. The weak variable region extraction module is connected to the performance evaluation module and is used to filter low value region grid points according to the vector mean threshold to generate weak variable region spatial data. The fault discrimination module, connected to the weak variation zone extraction module, is used to determine the fault nature based on the directional features of the vertical vector data on both sides of the fault and generate fault nature discrimination results. The source area output module is connected to the weak variation zone extraction module and the fault discrimination module. It is used to determine the spatial range of the potential source area by combining the spatial data of the weak variation zone and the fault property discrimination results, and output the potential source area data.

Images