A method for analyzing the spatio-temporal variation of earthquake geomagnetic field for short-term and impending earthquake monitoring

Abstract

The present invention discloses a method for analyzing the spatio-temporal variation of the earthquake geomagnetic field for short-term and impending earthquake monitoring, which is as follows. First, the geomagnetic stations within a given range of longitude and latitude are divided into multiple rectangular grids according to longitude and latitude. The multi-component data of each station within the range of longitude and latitude on the detection day and the previous 60 days are preprocessed into multiple segments of standard time series. Then, the daily absolute mean error of each component prediction sequence is taken, and the difference sequence is subjected to anomaly detection through a sliding interquartile range detection algorithm to obtain the daily anomaly value of each component. Then, according to the daily geomagnetic disturbance index, the anomalies caused by magnetic storms are removed. Finally, the daily anomaly index matrix of each station is generated according to the grid size, and the sum of the anomaly index matrices of all stations on the detection day and the previous 29 days is obtained to get the final regional daily anomaly index matrix, and the range of anomalous longitude and latitude is obtained. This method monitors the anomalous changes of the geomagnetic field within the region and provides a reference range for the geographical locations where earthquakes may occur.

Description

Technical Field

[0001] The present invention relates to a method for analyzing the spatio-temporal variation of the earthquake geomagnetic field for short-term and impending earthquake monitoring, belonging to the field of detecting precursor anomalies of the earthquake geomagnetic field. Background Art

[0002] An earthquake is a natural phenomenon caused by vibrations on or beneath the Earth's surface, such as crustal activity and volcanic activity. Disastrous earthquakes are generally caused by sudden movements of the Earth's crust. When a disastrous earthquake occurs without warning, it causes great harm. The building collapses, ground subsidence, and landslides caused by earthquakes pose a huge threat to the property and lives of the people. With the advancement of the digital construction of the seismic network and the standardization of seismic electromagnetic equipment, and the in-depth study of existing earthquake cases, geomagnetic anomaly signals are considered to be a reliable earthquake precursor signal. Especially before a major earthquake, the geomagnetic anomaly signal contains potential earthquake precursor information. Currently, the existing methods for extracting geomagnetic field anomalies are all based on the good synchronization of the daily variation of the vertical component of the geomagnetic field within a certain spatial range. By removing or suppressing the normal daily variation of the geomagnetic field, earthquake anomalies are extracted, and the spatio-temporal variation characteristics of the earthquake anomalies and their relationship with earthquakes are analyzed and studied. These methods cannot make full use of the existing data, and the calculation time span of the anomalies is long, and the timeliness is poor. Moreover, most of the current geomagnetic analysis methods are based on single stations for analysis, without considering the spatial relationship of the anomalies. Therefore, the existing methods for detecting and analyzing earthquake geomagnetic field anomalies need to be improved. Summary of the Invention

[0003] Aiming at the above problems existing in the prior art, the present invention aims to provide a method for analyzing the spatio-temporal variation of the earthquake geomagnetic field for short-term and impending earthquake monitoring. By detecting anomalies in the multi-component data of geomagnetic stations, spatio-temporal anomaly monitoring of multiple geomagnetic stations within a given range of longitude and latitude is realized, and on this basis, a geographical longitude and latitude range where an earthquake is most likely to occur is provided.

[0004] To achieve the object of the present invention, the technical solution adopted by the present invention is as follows:

[0005] A method for analyzing the spatio-temporal variation of the earthquake geomagnetic field for short-term and impending earthquake monitoring, the method comprising the following steps:

[0006] Step 1: Divide the longitude and latitude at equal intervals and generate an array of longitude coordinates and an array of latitude coordinates;

[0007] Step 2: Preprocess the multi-component data of a single station within the grid to generate multiple segments of standard time series;

[0008] Step 3: Input the time series into a prediction algorithm based on a long short-term neural network and an autoencoder to obtain a predicted time series;

[0009] Step 4: Calculate the daily absolute mean error of the predicted time series and perform sliding interquartile range anomaly detection to obtain the outlier time series;

[0010] Step 5: Remove the anomalies on the days of magnetic storms from the anomaly results and calculate the station anomaly index matrix;

[0011] Step 6: Repeat the process of Step 2 to Step 5 until all stations within the longitude and latitude range are calculated. Then sum the anomaly index matrices of all stations to obtain the cumulative anomaly index matrix, and output the geomagnetic anomaly heat map of this region.

[0012] As an improvement of the present invention, in the said Step 1, the method for equally spaced segmentation of longitude and latitude and generating the longitude coordinate array and latitude coordinate array is:

[0013] Divide the meridians and parallels by 0.1 degrees. Through the upper and lower limits of the input meridians and parallels, a rectangular grid composed of multiple small squares with a side length of 0.1 degrees can be obtained. Then generate a longitude coordinate array and a latitude coordinate array corresponding to the number of longitude and latitude points in the rectangular grid according to the number of longitude and latitude points in the grid.

[0014] As an improvement of the present invention, in the said Step 2, the single-component preprocessing method for a single station is:

[0015] Take the time series data of the vertical component, horizontal component, total amount, and magnetic declination of the station on the current day and the previous 60 days, and respectively downsample the time series data of each component. Divide a day into 96 time windows with a 15-minute time window, and then calculate the average value of the observation points within the window. Through downsampling, the geomagnetic data of a day is changed to 96 points. Then normalize the data to the interval [0, 1]. Set a sliding window with a length of 96, and slide on the normalized data with a step of 1 point to obtain multiple segments of input time series, so as to obtain the input data of each component.

[0016] As an improvement of the present invention, in the said Step 3, the prediction algorithm model based on the long short-term neural network and the autoencoder is:

[0017] The model consists of three parts: an encoder, a decoder, and a fully connected layer network, namely the LSTM-AE model. The encoder consists of an LSTM layer and a vector replication module. The input is encoded into a vector output through the LSTM network, and then the vector is replicated through the replication module to achieve the same dimension as the input requirement. The decoder consists of an LSTM layer identical to the encoder and a fully connected layer. The LSTM layer and the fully connected layer are used to convert the encoded vector into a reconstructed vector with the same dimension as the input sequence. Finally, a linear fully connected layer is used to convert the vector output by the decoder into a predicted value. The model adopts a semi-supervised learning method. By allowing the LSTM-AE model to learn the time series without earthquakes, it predicts the values of geomagnetic field data under normal conditions. Each component trains its LSTM-AE model separately.

[0018] The calculation process of the predicted time series for a single component is as follows:

[0019] First, the preprocessed data of each component is sequentially input into the trained LSTM-AE model. Then, the predicted time series data for the current day and the previous 59 days is obtained. Finally, the predicted data is de-normalized using the normalization parameters during model training to obtain the final predicted time series.

[0020] After each component data executes the above process, the time series of the vertical component can be obtained. The time series of the horizontal component The time series of the total amount The time series of the magnetic declination

[0021] As an improvement of the present invention, in step 4, the method for calculating the daily absolute mean error of the predicted time series is as follows:

[0022] Calculate the absolute mean error within a single day for each component in chronological order. The calculation formula is as follows:

[0023]

[0024] Where is the average absolute error on the t-th day, represents the predicted value at the i-th moment on the t-th day, and y i represents the observed value at the i-th moment on the t-th day. Then, the time series of the daily differences of the vertical component can be obtained. The time series of the daily differences of the horizontal component The time series of the daily differences of the total amount The time series of the daily differences of the magnetic declination Then, the sliding interquartile range method is used to detect each difference sequence to obtain the abnormal sequence of each component. The sliding interquartile range algorithm is as follows:

[0025] The sliding window size is 30. The data within the window is sorted into an ascending array from smallest to largest, and then the first quartile Q 1 and the third quartile Q 3 and the median Q 2 are obtained. Then, Q 3 is subtracted from Q 1 to obtain the IQR. The magnetic declination uses 2 times the IQR as the threshold for anomaly detection, with a 99.3% confidence probability of detecting anomalies. For the remaining components, 2.5 times the IQR is used as the threshold for anomaly detection, with a 99.9% confidence probability of detecting anomalies. The upper and lower limits for outliers are determined as shown in the following formula:

[0026]

[0027] where up is the upper limit, low is the lower limit, and a is the IQR detection multiplier. If the data t at the next step of the sliding window is within the upper and lower limits, it is considered normal. If it exceeds the upper and lower limits, it is considered abnormal data, and the exceeded part is the outlier Δt. The formula is as follows:

[0028]

[0029] Then, with 1 day as the step size, the calculation of the outliers Δt for each component for the current day and the previous 29 days is completed by sliding, and the value of Δt is used as the detection result for that day.

[0030] As an improvement of the present invention, in step 5, the method for removing the anomaly calculation on the day of magnetic storm in the anomaly result is as follows:

[0031] Calculate the minimum value of the geomagnetic disturbance index every day. If the minimum value of the geomagnetic index on that day is less than -40 nT, it is considered that the anomaly is caused by a magnetic storm, and the anomaly value of each component on that day is set to 0. For the remaining days with an anomaly value greater than 0, if the geomagnetic disturbance index on that day is greater than -40 nT, the anomaly value is set to 1 to convert it into an anomaly index. The purpose is to eliminate the anomalies caused by magnetic storms and eliminate the differences between the anomaly values of each component.

[0032] The calculation method for the station anomaly index matrix is as follows:

[0033] First, calculate the sum of the anomaly indices of the vertical component Z, horizontal component H, total amount F, and magnetic declination D of a single station on a single day as the anomaly index of the station on that day, that is, the maximum single-day anomaly index is 4 and the minimum is 0. Then, use the haversine formula to calculate the distances between all longitude and latitude positions within the rectangular grid and the station. The calculation formula is as follows:

[0034]

[0035] where lat 1 and lat2 represent the latitudes of two points respectively, lon 1 and lon 2 represent the longitudes of two points respectively, and R is the radius of the earth. If the distance is less than or equal to 200 km, then add the abnormal index value of the station on that day to the value at the corresponding position of the longitude and latitude matrix of this location, so that the single-day abnormal index matrix X of the single station within the grid can be obtained t , and then add all the single-day abnormal index matrices to obtain the cumulative abnormal index matrix of the single station The overall single-station calculation algorithm in step 5 is shown as follows:

[0036]

[0037]

[0038] As an improvement of the present invention, in step 6, repeat the process of step 2 to step 5 until all stations within the longitude and latitude range are calculated, then the cumulative abnormal index matrix of all stations in the region can be obtained, and then sum the abnormal index matrices of all stations to obtain the cumulative abnormal index matrix, as shown in the following formula:

[0039]

[0040] where n is the total number of stations, is the final cumulative abnormal index matrix, and the longitude and latitude corresponding to the position with the largest abnormal index in the matrix is the position most likely to have an earthquake, and according to draw a geomagnetic anomaly heat map of this longitude and latitude region as the output.

[0041]

[0042] Compared with other methods, the method proposed by the present invention has the following beneficial effects:

[0043] First, it fuses the abnormal detection results of multiple components of geomagnetic field stations, rather than only analyzing with a single component of geomagnetism, making full use of all the information of geomagnetism, and the method is more universal. On the one hand, the abnormal detection method uses more historical observation information for abnormal detection, with stronger anti-interference ability and greater robustness. On the other hand, the abnormal detection method gives anomalies based on statistical probability rather than the thresholds set based on experience in the past, making it more interpretable.

[0044] Second, it fuses the spatio-temporal relationship of geomagnetism, uses the spatial abnormal information of all stations within a given longitude and latitude range and the abnormal information of each station in time for earthquake precursor analysis, and can reflect the process of geomagnetic anomaly change in the region before an earthquake to a certain extent, making it more interpretable.

[0045] Thirdly, it is possible to monitor the spatio-temporal variations of the geomagnetic field anomalies within a given range of longitude and latitude, and provide a reference range for the geographical locations where earthquakes may occur. Description of the Drawings

[0046] Figure 1 It is a flow chart of a method for analyzing spatio-temporal variations of earthquake geomagnetic fields for short-term earthquake monitoring according to the present invention. Detailed Embodiments

[0047] The technical solution of the present invention will be further described below with reference to the drawings and embodiments.

[0048] Embodiment 1: Refer to Figure 1 A method for analyzing spatio-temporal variations of earthquake geomagnetic fields for short-term earthquake monitoring, the method comprising the following steps:

[0049] Step 1: Divide the longitude and latitude at equal intervals and generate a longitude coordinate array and a latitude coordinate array;

[0050] Step 2: Preprocess the multi-component data of a single station within the grid to generate multiple segments of standard time series;

[0051] Step 3: Input the time series into a prediction algorithm based on a long short-term neural network and an autoencoder to obtain a predicted time series;

[0052] Step 4: Calculate the daily absolute mean error of the predicted time series and perform a moving interquartile range anomaly detection to obtain an anomaly value time series;

[0053] Step 5: Remove the anomalies on the days of magnetic storms in the anomaly results and calculate the station anomaly index matrix;

[0054] Step 6: Repeat the process of Step 2 to Step 5 until all stations within the longitude and latitude range are calculated. Then sum the anomaly index matrices of all stations to obtain an accumulated anomaly index matrix, and output the geomagnetic anomaly heat map of this region.

[0055] In the said Step 1, the method for dividing the longitude and latitude at equal intervals and generating a longitude coordinate array and a latitude coordinate array is:

[0056] Divide the longitude lines and latitude lines by 0.1 degree. Through the upper and lower limits of the input longitude and latitude, a rectangular grid composed of multiple small squares with a side length of 0.1 degree can be obtained. Then generate a longitude coordinate array and a latitude coordinate array corresponding to the number of longitude and latitude points within the rectangular grid according to the number of longitude and latitude points within the grid.

[0057] In the said Step 2, the preprocessing method for a single component of a single station is:

[0058] Retrieve the time series data of the vertical component, horizontal component, total amount, and magnetic declination for the current day and the previous 60 days of the station. Then, downsample the time series data of each component. Divide a day into 96 time windows with a 15-minute time window, and then calculate the average value of the observation points within the window. Through downsampling, the geomagnetic data for a day is changed to 96 points. Then, normalize the data to the range [0, 1]. Set a sliding window of length 96 and slide it with a step size of 1 point on the normalized data to obtain multiple segments of input time series, thereby obtaining the input data for each component.

[0059] In step 3, the prediction algorithm model based on the long short-term neural network and autoencoder is as follows:

[0060] This model consists of three parts: an encoder, a decoder, and a fully connected layer network, that is, the LSTM-AE model. The encoder consists of an LSTM layer and a vector replication module. The input is encoded into a vector output through the LSTM network, and then through the replication module, this vector is replicated to achieve the same dimension as the input requirement. The decoder consists of an LSTM layer identical to the encoder and a fully connected layer. The LSTM layer and this fully connected layer are used to convert the encoded vector into a reconstructed vector with the same dimension as the input sequence. Finally, a linear fully connected layer is used to convert the vector output by the decoder into a predicted value. The model adopts a semi-supervised learning method. By allowing the LSTM-AE model to learn the time series during earthquake-free periods to predict the values of the geomagnetic field data under normal conditions, the LSTM-AE model for each component is trained separately.

[0061] The calculation process of the prediction time series for a single component is as follows:

[0062] First, input the preprocessed data of each component into the trained LSTM-AE model in sequence, then obtain the prediction time series data for the current day and the previous 59 days, and finally, denormalize the prediction data using the normalization parameters during model training to obtain the final prediction time series.

[0063] After each component data has completed the above process, the time series of the vertical component can be obtained. The time series of the horizontal component The time series of the total amount The time series of the magnetic declination

[0064] In step 4, the method for calculating the daily absolute mean error of the prediction time series is as follows:

[0065] Calculate the absolute mean error within a single day for each component in chronological order. The calculation formula is as follows:

[0066]

[0067] where is the average absolute error on the t-th day, represents the predicted value at the i-th moment on the t-th day, y i represents the observed value at the i-th moment on the t-th day. Then, the time series of the daily differences of the vertical component can be obtained The time series of the daily differences of the horizontal component The time series of the daily differences of the total amount The time series of the daily differences of the magnetic declination Then, the sliding interquartile range method is used to detect each difference sequence to obtain the abnormal sequences of each component. The sliding interquartile range algorithm is as follows:

[0068] The sliding window size is 30. The data within the window is arranged in ascending order to form an ascending array, and then the first quartile Q 1 and the third quartile Q 3 and the median Q 2 are obtained. Then, Q 3 is subtracted from Q 1 to obtain the IQR. For the magnetic declination, 2 times the IQR is used as the threshold for anomaly detection, with a 99.3% confidence probability of detecting anomalies. For the other components, 2.5 times the IQR is used as the threshold for anomaly detection, with a 99.9% confidence probability of detecting anomalies. The upper and lower limits for determining the outliers are shown in the following formula:

[0069]

[0070] where up is the upper limit, low is the lower limit, and a is the IQR detection magnification factor. If the data t at the next step of the sliding window is within the upper and lower limits, it is considered normal; if it exceeds the upper and lower limits, it is considered abnormal data, and the exceeded part is the outlier Δt. The formula is as follows:

[0071]

[0072] Then, with 1 day as the step size, the calculation of the outliers Δt of each component for the current day and the previous 29 days is completed by sliding, and the value of Δt is used as the detection result for that day.

[0073] In step 5, the abnormal calculation method for removing the days of magnetic storms in the abnormal results is as follows:

[0074] Calculate the daily minimum value of the geomagnetic disturbance index. If the minimum value of the geomagnetic index on that day is less than -40 nT, it is considered that the anomaly is caused by a magnetic storm, and the outliers of each component on that day are set to 0. For the remaining days with outliers greater than 0, if the geomagnetic disturbance index on that day is greater than -40 nT, the outliers are set to 1 to convert them into anomaly indices. The purpose is to eliminate the anomalies caused by magnetic storms and eliminate the differences between the outliers of each component.

[0075] The calculation method of the station anomaly index matrix is as follows:

[0076] First, calculate the sum of the anomaly indices of the vertical component Z, horizontal component H, total amount F, and magnetic declination D of a single station on a single day as the anomaly index of the station on that day, that is, the maximum anomaly index on a single day is 4 and the minimum is 0. Then, use the haversine formula to calculate the distances between all longitude and latitude positions within the rectangular grid and the station. The calculation formula is as follows:

[0077]

[0078] where lat 1 and lat 2 represent the latitudes of two points respectively, lon 1 and lon 2 represent the longitudes of two points respectively, and R is the radius of the Earth. If the distance is less than or equal to 200 km, then add the anomaly index value of the station on that day to the value at the corresponding longitude and latitude matrix position. In this way, the single-day anomaly index matrix X t of a single station within the grid can be obtained. Then, add up all the single-day anomaly index matrices to get the cumulative anomaly index matrix of the single station The overall calculation algorithm for a single station in step 5 is shown as follows:

[0079]

[0080]

[0081] In step 6, repeat the process from step 2 to step 5 until all stations within the longitude and latitude range are calculated. Then, the cumulative anomaly index matrix of all stations in the region can be obtained. Then, sum up the anomaly index matrices of all stations to get the cumulative anomaly index matrix, as shown in the following formula:

[0082]

[0083] where n is the total number of stations, is the final cumulative anomaly index matrix. The longitude and latitude corresponding to the position with the largest anomaly index in the matrix is the most likely location of an earthquake, and a geomagnetic anomaly heat map of the longitude and latitude region is drawn as the output according to

[0084]

[0085]

[0086] ​It should be noted that the above embodiments are only the preferred embodiments of the present invention and do not limit the protection scope of the present invention. Any equivalent replacement or substitution made on the basis of the above technical solutions falls within the protection scope of the present invention.

Claims

1. A method for analyzing the spatio-temporal variation of seismic geomagnetic field for short-term earthquake monitoring, characterized in that, the method comprises the following steps: Step 1: Divide the longitude and latitude at equal intervals and generate a longitude coordinate array and a latitude coordinate array; Step 2: Preprocess the multi-component data of a single station within the grid to generate multiple segments of standard time series; Step 3: Input the time series into a prediction algorithm based on a long short-term neural network and an autoencoder to obtain a predicted time series; Step 4: Calculate the daily absolute mean error of the predicted time series and perform sliding interquartile range anomaly detection to obtain an outlier time series; Step 5: Remove the anomalies on the days of magnetic storms in the anomaly results and calculate the station anomaly index matrix; Step 6: Repeat the process of Step 2 to Step 5 until all stations within the longitude and latitude range are calculated. Then sum the anomaly index matrices of all stations to obtain a cumulative anomaly index matrix, and output the regional geomagnetic anomaly heat map corresponding to the longitude and latitude.

2. A method for analyzing the spatio-temporal variation of seismic geomagnetic field for short-term earthquake monitoring according to claim 1, characterized in that: in the said Step 1, the method for dividing the longitude and latitude at equal intervals and generating a longitude coordinate array and a latitude coordinate array is: Divide the meridian and parallel by 0.1 degrees. Through the upper and lower limits of the input meridian and parallel, a rectangular grid composed of multiple small squares with a side length of 0.1 degrees can be obtained. Then generate a longitude coordinate array and a latitude coordinate array corresponding to the number of longitude and latitude points in the grid according to the number of longitude and latitude points in the rectangular grid.

3. A method for analyzing the spatio-temporal variation of seismic geomagnetic field for short-term earthquake monitoring according to claim 1, characterized in that: in the said Step 2, the preprocessing method for a single component of a single station is: Take the time series data of the vertical component, horizontal component, total amount and magnetic declination of the station on the current day and the previous 60 days, and respectively downsample the time series data of each component. Divide a day into 96 time windows with a 15-minute time window, and then calculate the average value of the observation points within the window. Through downsampling, the geomagnetic data of a day is changed to 96 points, and then the data is normalized to the interval [0,1]. Set a sliding window with a length of 96, and slide on the normalized data with 1 point as the step length to obtain multiple segments of input time series, so as to obtain the input data of each component.

4. A method for analyzing the spatio-temporal variation of seismic geomagnetic field for short-term earthquake monitoring according to claim 1, characterized in that: in the said Step 3, the prediction algorithm model based on a long short-term neural network and an autoencoder is: The model consists of three parts: encoder, decoder, and fully connected layer network, namely LSTM-AE model. The encoder consists of an LSTM layer and a vector copy module. The input is encoded into a vector output through the LSTM network, and then the vector is copied through the copy module to achieve the same dimension as the input requirement. The decoder consists of an LSTM layer and a fully connected layer that are the same as the encoder. The LSTM layer and the fully connected layer are used to convert the encoded vector into a reconstructed vector with the same dimension as the input sequence. Finally, a linear fully connected layer is used to convert the decoder output vector into a predicted value. The model adopts a semi-supervised learning method. The LSTM-AE model is used to learn the time series when there is no earthquake to predict the value of the geomagnetic field data under normal conditions. Each component trains its LSTM-AE model separately. The calculation process of the forecast time series of a single component is: First, input the preprocessed data of each component into the trained LSTM-AE model in sequence, then obtain the predicted time series data for the current day and the previous 59 days, and finally, perform inverse normalization on the predicted data using the normalization parameters during model training to obtain the final predicted time series After the above process is executed for each component data, the time series of the vertical component can be obtained The time series of the horizontal component The time series of the total amount The time series of the magnetic declination 5. The method for analyzing the temporal and spatial variation of the seismic geomagnetic field for short-term and impending earthquake monitoring according to claim 1, Features: In step 4, the method for calculating the daily absolute average error of the forecast time series is: The absolute average error of each component within a single day is calculated in chronological order. The calculation formula is as follows: Among them is the average absolute error on the t-th day, represents the predicted value at the i-th moment on the t-th day, y i represents the observed value at the i-th moment on the t-th day, and then the time series of the daily difference of the vertical component can be obtained The time series of the daily difference of the horizontal component The time series of the daily difference of the total amount The time series of the daily difference of the magnetic declination Then, the sliding interquartile range method is used to detect each difference sequence to obtain the abnormal sequence of each component. The sliding interquartile range algorithm is as follows: The sliding window size is 30. The data within the window is arranged in ascending order to form an ascending array, and then the first quartile Q 1 and the third quartile Q 3 and the median Q 2 are obtained. Then, Q 3 is subtracted from Q 1 to obtain the IQR. The magnetic declination uses 2 times the IQR as the threshold for anomaly detection, with a 99.3% confidence probability of detecting anomalies. The remaining components use 2.5 times the IQR as the threshold for anomaly detection, with a 99.9% confidence probability of detecting anomalies. The upper and lower limits of the outliers are determined as shown in the following formula: Where up is the upper limit, low is the lower limit, and a is the IQR detection magnification. If the data t of one step after the sliding window is within the upper and lower limits, it is considered to be normal. If it exceeds the upper and lower limits, it is considered to be abnormal data, and the excess part is the abnormal value Δt. The formula is as follows: Then, with a step length of 1 day, the calculation of the abnormal value Δt of each component on the day and the 29 days before the day is completed, and the Δt value is used as the detection result of the day.

6. The method for analyzing the temporal and spatial variation of the seismic geomagnetic field for short-term and impending earthquake monitoring according to claim 1, Features: In step 5, the calculation method for removing the abnormal result of the day when the magnetic storm occurs is: Calculate the geomagnetic disturbance index, the daily minimum value. If the minimum value of the geomagnetic index on that day is less than -40nT, it is considered that the anomaly is caused by a magnetic storm, and the anomaly value of each component on that day is set to 0. For the other days when the anomaly value is greater than 0, if the geomagnetic disturbance index on that day is greater than -40nT, the anomaly value is set to 1 to convert it into an anomaly index. The purpose is to eliminate the anomalies caused by magnetic storms and eliminate the differences between the anomaly values ​​of each component. The calculation method of the station anomaly index matrix is: First, calculate the sum of the anomaly indexes of the vertical component Z, horizontal component H, total F and magnetic declination D of a single station on a single day as the anomaly index of the station on that day, that is, the maximum anomaly index of a single day is 4 and the minimum is 0. Then use the haversine formula to calculate the distance between all longitude and latitude positions in the rectangular grid and the station. The calculation formula is as follows: where lat 1 and lat 2 represent the latitudes of two points respectively, lon 1 and lon 2 represent the longitudes of two points respectively, R is the radius of the earth. If the distance is less than or equal to 200 km, then add the anomaly index value of the station on that day to the value at the corresponding position of the longitude and latitude matrix of the location, to obtain the single-station single-day anomaly index matrix X t , and then add all the single-day anomaly index matrices to obtain the cumulative anomaly index matrix of the single station The overall single-station calculation algorithm in step 5 is shown as follows:

7. The method for analyzing the temporal and spatial variation of the seismic geomagnetic field for short-term and impending earthquake monitoring according to claim 1, Features: In step 6, the process from step 2 to step 5 is repeated until all stations within the latitude and longitude range are calculated, and the cumulative anomaly index matrix of all stations in the region can be obtained. Then, the anomaly index matrices of all stations are summed to obtain the cumulative anomaly index matrix, as shown in the following formula: where n is the total number of stations, is the final cumulative anomaly index matrix. The longitude and latitude corresponding to the position with the largest anomaly index in the matrix are the most likely earthquake occurrence positions, and based on a geomagnetic anomaly heat map of the longitude and latitude region is drawn as the output.

Images