Pre-earthquake geomagnetic anomaly detection method based on support vector machine

Abstract

The invention discloses a pre-earthquake geomagnetic anomaly detection method based on a support vector machine. Comprising the following steps: 1) obtaining a residual sequence of a geomagnetic vertical component; 2) calculating a geomagnetic daily change template; 3) calculating geomagnetic distortion characteristics monitored by the station; 4) constructing a data set; and 5) establishing a geomagnetic anomaly detection model based on a support vector machine. The method relates to the technical fields of earthquake forecasting, geomagnetic anomaly detection, machine learning and the like, has the advantages of being simple, small in calculation amount, labor-saving, easy to implement and the like, and is suitable for geomagnetic change anomaly detection.

Description

Technical Field

[0001] The present invention relates to the technical fields of earthquake prediction, geomagnetic anomaly detection and machine learning, and in particular to a method for detecting geomagnetic anomaly before an earthquake based on a support vector machine. Background Art

[0002] Earthquakes are vibrations caused by the rapid release of energy from the earth's crust. They are natural disasters that are extremely harmful to human society. A large number of studies have shown that many physical quantities will undergo abnormal changes before an earthquake occurs, which are called earthquake precursors. Exploring the correlation between abnormal precursor phenomena and earthquakes is of great significance for understanding and predicting earthquakes.

[0003] At present, the total number of operating seismic stations in China has exceeded 1,000, including strain gauges, magnetometers, geoelectric instruments, water level meters and other instruments. These instruments monitor changes in factors such as crustal deformation, geoelectricity, geomagnetism, gravity and underground fluids. The accumulated data show that there are abnormal phenomena before earthquakes. It is a meaningful exploration to extract possible precursor information based on data from a period of time before an earthquake occurs to assist in earthquake prediction.

[0004] The development and occurrence of an earthquake is an extremely complex process, which often causes abnormal changes in geomagnetic signals, known as seismomagnetic phenomena. Affected by the earthquake, the geomagnetic anomaly changes more dramatically in areas near the epicenter, while the geomagnetic anomaly changes relatively less in areas farther from the epicenter. At present, the main methods for determining geomagnetic anomalies before an earthquake include the loading and unloading response ratio method, the daily ratio method, the daily variation spatial correlation method, the low point displacement method, the daily value difference method, the harmonic amplitude ratio, etc. Summary of the invention

[0005] The purpose of the present invention is to use the geomagnetic monitoring data before an earthquake to train a model for geomagnetic anomaly detection. The method has the advantages of simple method, small amount of calculation, labor saving, easy implementation, etc., and is suitable for geomagnetic change anomaly detection.

[0006] The present invention comprises the following steps:

[0007] (1) Obtain the residual sequence of the geomagnetic vertical component:

[0008] According to y residual (t) = y(t) - y trend (t) Calculate the residual sequence y of the vertical component of the geomagnetic field residual (t), where y(t) is the geomagnetic vertical component data monitored by the seismic station, y trend (t) = a0 + a1t + a2t 2 +…+a s t s is the trend of the vertical component of the geomagnetic field. Based on y(t), the least square method is used to determine a0, a1, …, as The value of; s is the order of the polynomial, which is set manually; t=1,2,…,N, N is the number of data points of y(t);

[0009] (2) Calculation template for daily geomagnetic variation:

[0010] according to y residual (t)

[0011] Split by day; if the jth day satisfies Then the jth day is a magnetically quiet day; according to Calculate the daily variation template of geomagnetic field template , where m is the number of data points for daily geomagnetic monitoring, q = N / m, Y i is the geomagnetic monitoring data of magnetically quiet days, n is the total number of magnetically quiet days, j = 1, 2, ..., q, Dst j is the Dst index of the jth day, obtained from the World Geomagnetic Data Center (Kyoto); Dst thr is a parameter, set manually;

[0012] (3) Calculate the geomagnetic distortion characteristics of the station monitoring:

[0013] (a) Calculate the time offset d of the geomagnetic low point r,1 =[d1(r,1),d1(r,2),...,d1(r,q)], where

[0014] d1(r,j)=|T(min(y j ))-T(min(y template ))|,T(min(y j )) is the time of the geomagnetic low point at station r on the jth day, T(min(y template ) is the low point time of the geomagnetic daily variation template of station r; r = 1, 2, ..., num, num is the number of stations;

[0015] (b) Calculate the Chebyshev distance d r,2 =[d2(r,1),d2(r,2),…,d2(r,q)], where d2(r,j)=max|y j -y template | is the Chebyshev distance of station r on the jth day;

[0016] (c) Calculate the degree of daily geomagnetic variation distortion d r,3 =[d3(r,1),d3(r,2),…,d3(r,q)], where, g(r,j)=|(max(y j )-min(y j ))-(max(ytemplate )-min(y template ))|, max() and min() represent the maximum and minimum values ​​respectively; a, b and c are parameters set manually; (4) Construct the data set:

[0017] (a) The input features of the dataset are Among them, μ k and σ k f k The mean and standard deviation of the geomagnetic diurnal variation are shown in Table 1. k = 1, 2, 3; f1 = [f1(1), f1(2), …, f1(num)] is the time offset of the geomagnetic low point, f2 = [f2(1), f2(2), …, f2(num)] is the Chebyshev distance, and f3 = [f3(1), f3(2), …, f3(num)] is the degree of geomagnetic diurnal variation distortion; f1(r) and f2(r) are d r,1 and d r,2 The maximum deviation value of f3(r) is d r,3 The maximum value of

[0018] (b) The output of the dataset is the station label; according to R = 10 g×M Calculate the earthquake impact range R. If the distance from the station to the epicenter is less than R, the station label is 1, otherwise, the station label is 0, where M is the magnitude and g is a parameter set manually;

[0019] (5) Establish a geomagnetic anomaly detection model based on support vector machine. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 The original geomagnetic data and fitted trend components monitored by Taian Station;

[0021] Figure 2 is the geomagnetic residual sequence of Taian station;

[0022] Figure 3 This is the change diagram of the Dst index, an indicator of the magnitude of the magnetic storm before and after the earthquake;

[0023] Figure 4 This is the daily variation template of the geomagnetic vertical component at Taian Station;

[0024] Figure 5 This is the time variation diagram of the geomagnetic low point offset at Taian Station;

[0025] Figure 6 This is the change diagram of geomagnetic Chebyshev distance at Taian Station;

[0026] Figure 7 This is the change diagram of the degree of distortion of the geomagnetic diurnal variation at Taian Station;

[0027] Figure 8 and Fig. 9 The geomagnetic anomaly detection model was established. DETAILED DESCRIPTION

[0028] This embodiment uses the geomagnetic data monitored by 20 seismic stations near the epicenter of the 8.6 plain earthquake in 2023, and establishes a geomagnetic anomaly detection model based on a support vector machine to detect whether the geomagnetic changes are abnormal.

[0029] The specific steps are as follows:

[0030] (1) Obtain the residual sequence of the geomagnetic vertical component:

[0031] According to y residual (t) = y(t) - y trend (t) Calculate the residual sequence y of the vertical component of the geomagnetic field residual (t), where y(t) is the geomagnetic vertical component data monitored by the seismic station, y trend (t) = a0 + a1t + a2t 2 +…+a s t s is the trend of the vertical component of the geomagnetic field. Based on y(t), the least square method is used to determine a0, a1, …, a s The value of ; s is the order of the polynomial, which is set manually; t=1,2,…,N, N is the number of data points of y(t).

[0032] In this embodiment, the magnitude 8.6 plain earthquake in 2023 is taken as an example, and the earthquake intensity is 5.5. Taking Tai'an station as an example, N = 100800, s = 6, a0 = 42671.0224, a1 = -3.3679, a2 = -0.2708, a3 = 5.6471, a4 = -0.5699, a5 = -1.3338, a6 = 0.3149, the fitted trend component and the monitored original geomagnetic data are as follows Figure 1 As shown, the geomagnetic residual sequence is as follows Figure 2 shown.

[0033] (2) Calculation template for daily geomagnetic variation:

[0034] according to y residual (t)

[0035] Split by day; if the jth day satisfies Then the jth day is a magnetically quiet day; according to Calculate the daily variation template of geomagnetic field template , where m is the number of data points for daily geomagnetic monitoring, q = N / m, Y i is the geomagnetic monitoring data of magnetically quiet days, n is the total number of magnetically quiet days, j = 1, 2, ..., q, Dst jis the Dst index of the jth day, obtained from the World Geomagnetic Data Center (Kyoto); Dst thr It is a parameter set manually.

[0036] In this embodiment, taking Tai'an Station as an example, the geomagnetic change has a periodicity of 24 hours, the sampling interval is 1 minute, m=1440, q=70, Dst thr = -30. If the jth day satisfies Then this day is a magnetically quiet day. In this embodiment, the total number of magnetically quiet days is n=42. Figure 3 is the change of Dst index, Figure 4 This is the template for the daily geomagnetic variation of Tai'an Station.

[0037] (3) Calculate the geomagnetic distortion characteristics of the station monitoring:

[0038] (a) Calculate the time offset d of the geomagnetic low point r,1 =[d1(r,1),d1(r,2),…,d1(r,q)], where d1(r,j)=|T(min(y j ))-T(min(y template ))|,T(min(y j )) is the time of the geomagnetic low point at station r on the jth day, T(min(y template ) is the low point time of the geomagnetic diurnal variation template of station r; r = 1, 2, …, num, num is the number of stations;

[0039] (b) Calculate the Chebyshev distance d r,2 =[d2(r,1),d2(r,2),…,d2(r,q)], where d2(r,j)=max|y j -y template | is the Chebyshev distance of station r on the jth day;

[0040] (c) Calculate the degree of daily geomagnetic variation distortion d r,3 =[d3(r,1),d3(r,2),…,d3(r,q)], where, g(r,j)=|(max(y j )-min(y j ))-(max(y template )-min(y template ))|, max() and min() represent the maximum and minimum values ​​respectively; a, b and c are parameters, which are set manually;

[0041] In this embodiment, a=120, b=15, c=1. Figure 5-Figure 7They are the time variation graph of geomagnetic low point offset, the variation graph of geomagnetic Chebyshev distance and the variation graph of geomagnetic diurnal variation distortion degree at Tai'an station.

[0042] (4) Constructing the dataset:

[0043] (a) The input features of the dataset are Among them, μ k and σ k f k The mean and standard deviation of the geomagnetic diurnal variation are shown in Table 1. k = 1, 2, 3; f1 = [f1(1), f1(2), …, f1(num)] is the time offset of the geomagnetic low point, f2 = [f2(1), f2(2), …, f2(num)] is the Chebyshev distance, and f3 = [f3(1), f3(2), …, f3(num)] is the degree of geomagnetic diurnal variation distortion; f1(r) and f2(r) are d r,1 and d r,2 The maximum deviation value of f3(r) is d r,3 The maximum value of

[0044] (b) The output of the dataset is the station label; according to R = 10 g×M Calculate the earthquake impact range R. If the distance from the station to the epicenter is less than R, the station label is 1, otherwise, the station label is 0, where M is the magnitude and g is a parameter set manually;

[0045] In this embodiment, M=5.5, g=0.43, and R is calculated to be 232 kilometers.

[0046] (5) Establish a geomagnetic anomaly detection model based on support vector machine.

[0047] The data set is divided into a training set and a test set. The support vector machine method is used to learn the training set data and establish a geomagnetic anomaly detection model. The established geomagnetic anomaly detection model is used to detect the geomagnetic changes monitored by the station to determine whether there are anomalies.

[0048] In this embodiment, the data monitored by 20 stations are used for modeling, of which the data monitored by 14 stations are used for training, and the data monitored by the remaining 6 stations are used for testing. The kernel function of the support vector machine selects the radial basis function, and the penalty parameter value is 10. Based on the established model, the geomagnetic changes monitored by the stations are detected for anomalies. Figure 8 is the detection result of geomagnetic anomaly based on the two features f1 and f2. Fig. 9 This is the result of geomagnetic change anomaly detection based on the two features f2 and f3.

Claims

1. A method for detecting geomagnetic anomalies before earthquakes based on support vector machine. Its characteristics include the following specific steps: (1) Obtain the residual sequence of the geomagnetic vertical component: According to y residual (t) = y(t) - y trend (t) Calculate the residual sequence y of the vertical component of the geomagnetic field residual (t), where y(t) is the vertical component data of geomagnetic field monitored by the seismic station, y trend (t) = a0 + a1t + a2t 2 +…+a s t s is the trend of the vertical component of the geomagnetic field. Based on y(t), the least square method is used to determine a0, a1, …, a s The value of; s is the order of the polynomial, which is set manually; t=1,2,...,N, N is the number of data points of y(t); (2) Calculation template for daily geomagnetic variation: according to y residual (t) Split by day; if the jth day satisfies Then the jth day is a magnetically quiet day; according to Calculate the daily variation template of geomagnetic field template , where m is the number of data points for daily geomagnetic monitoring, q = N / m, Y i is the geomagnetic monitoring data of magnetically quiet days, n is the total number of magnetically quiet days, j = 1, 2, ..., q, Dst j is the Dst index of the jth day, obtained from the World Geomagnetic Data Center (Kyoto); Dst thr is a parameter, set manually; (3) Calculate the geomagnetic distortion characteristics of the station monitoring: (a) Calculate the time offset d of the geomagnetic low point r,1 =[d1(r,1),d1(r,2),...,d1(r,q)], where d1(r,j)=|T(min(y j ))-T(min(y template ))|,T(min(y j )) is the time of the geomagnetic low point at station r on the jth day, T(min(y template ) is the low point time of the geomagnetic diurnal variation template of station r; r = 1, 2, …, num, num is the number of stations; (b) Calculate the Chebyshev distance d r,2 =[d2(r,1),d2(r,2),...,d2(r,q)], where d2(r,j)=max|y j -y template | is the Chebyshev distance of station r on the jth day; (c) Calculate the degree of distortion of the geomagnetic diurnal variation d r,3 =[d3(r,1),d3(r,2),...,d3(r,q)], where, g(r,j)=|(max(y j )-min(y j ))-(max(y template )-min(y template ))|, max() and min() represent the maximum and minimum values ​​respectively; a, b and c are parameters, which are set manually; (4) Constructing the dataset: (a) The input features of the dataset are Among them, μ k and σ k f k The mean and standard deviation of the geomagnetic diurnal variation are shown in Table 1. k = 1, 2, 3; f1 = [f1(1), f1(2), …, f1(num)] is the time offset of the geomagnetic low point, f2 = [f2(1), f2(2), …, f2(num)] is the Chebyshev distance, and f3 = [f3(1), f3(2), …, f3(num)] is the degree of geomagnetic diurnal variation distortion; f1(r) and f2(r) are d r,1 and d r,2 The maximum deviation value of f3(r) is d r,3 The maximum value of (b) The output of the dataset is the station label; according to R = 10 g×M Calculate the earthquake impact range R. If the distance from the station to the epicenter is less than R, the station label is 1, otherwise, the station label is 0, where M is the magnitude and g is a parameter set manually; (5) Establish a geomagnetic anomaly detection model based on support vector machine.

Images