A forward earthquake prediction algorithm

Abstract

The present invention discloses a forward earthquake prediction algorithm, including: "AND operation" or "OR operation" for the following three conditions: Condition 1, the amplitude trend of the real-time waveform is constantly increasing; Condition 2, the amplitude jump at a specific frequency point of the real-time waveform; Condition 3, the real-time waveform is identified as a pre-earthquake waveform; When the "AND operation" is performed, if the above three conditions are all true, an earthquake warning is triggered; When the "OR operation" is performed, the operation mechanism is adjusted according to actual needs, and an earthquake warning is triggered when some of the conditions are true. Through the above method, the forward earthquake prediction algorithm described in the present invention can take into account the time domain characteristics and frequency domain changes of seismic waves, improve the accuracy of earthquake prediction, and reduce the problems of false alarms and missed alarms.

Description

Technical Field

[0001] The present invention relates to the technical field of earthquake prediction, and in particular to a forward earthquake prediction algorithm. Background Art

[0002] Earthquakes and other natural disasters caused by earthquakes usually cause serious casualties and property losses. In order to reduce the casualties and property losses caused by earthquakes, earthquake prediction is necessary to respond to the upcoming earthquake in advance.

[0003] Earthquake prediction is a common scientific problem throughout the world. Currently, the most widely used prediction method is to predict by looking at the expected time domain characteristics of the earthquake waveform. For example, if the amplitude of the waveform reaches a threshold, an earthquake is determined to have occurred, or the waveform characteristics of an earthquake are identified through machine learning model training. For example, the invention patent with application number 201910752213.2 discloses a method for identifying the vibration properties of earthquake waves based on machine learning to determine the category to which the earthquake waveform belongs.

[0004] Existing earthquake prediction algorithms usually only focus on one type of characteristics in the time domain, which can easily lead to false alarms or missed alarms. In addition, for prediction methods that focus on waveform amplitude, sudden changes in amplitude caused by environmental noise or human factors may trigger false alarms of earthquake warnings; for waveform models trained by machine learning, there may be misidentification or missed identification of earthquake waveforms, resulting in failure of early warnings, which needs to be improved. Summary of the invention

[0005] The main technical problem solved by the present invention is to provide a forward earthquake prediction algorithm, which predicts the time domain characteristics of the earthquake waveform and monitors the key frequency points in the frequency domain. It can take into account the time domain characteristics and frequency domain changes of the pre-seismic wave, improve the accuracy of earthquake prediction, and reduce false alarms and missed alarms.

[0006] In order to solve the above technical problems, a technical solution adopted by the present invention is to provide a forward earthquake prediction algorithm, including: an "AND operation" or "OR operation" for the following three conditions:

[0007] Condition 1: The amplitude trend of the real-time waveform is increasing;

[0008] Condition 2: The amplitude of the real-time waveform changes at a specific frequency point;

[0009] Condition 3: The real-time waveform is identified as the pre-earthquake waveform;

[0010] When performing “AND operation”, if all three conditions above are true, an earthquake warning will be triggered;

[0011] When performing the "OR operation", the operation mechanism is adjusted according to actual needs, and the earthquake warning is triggered when some of the conditions are true.

[0012] In a preferred embodiment of the present invention, the operation of condition 1 includes: predicting the amplitude trend of the real-time waveform and determining whether it is increasing continuously, and the specific steps are as follows:

[0013] A1. Create a stack of size 3, which is an array used to store the data segments to be processed;

[0014] A2. Receive real-time input data segment , perform the following operations on the most recently input data segment:

[0015] If the stack is not full, put the data segment at the top of the stack and update the top position of the stack;

[0016] If the stack is full, remove the data segment at the bottom of the stack first ,Will Put it on the top of the stack and update the position of the top of the stack;

[0017] A3. Calculation data segment , and The peak-to-peak value of , and ;

[0018] ,in For data segment All peak-to-peak values ​​of for The total number of peak-to-peak values;

[0019] ,

[0020] and Data segment Each peak value and its corresponding nearest valley value can be calculated by the same logic. and ;

[0021] A4. If , the output is True, otherwise the output is False;

[0022] A5. Receive the next data segment and repeat steps A2 to A4.

[0023] In a preferred embodiment of the present invention, the operation of condition 2 includes: in the frequency domain, monitoring the amplitude of the real-time waveform at a specific frequency point, and determining whether the amplitude has a jump. The specific steps are as follows:

[0024] B1. Create a stack of size 2, which is an array used to store two data segments to be processed;

[0025] B2. Receive real-time input data segment ,right and The Goertzel algorithm is used to calculate the amplitude of a specific frequency and store the calculation results separately. The calculation results of each frequency point are recorded as , The calculation result is recorded as ;

[0026] B3. If All elements in satisfy , the output is True, otherwise the output is False;

[0027] B4. Receive the next data segment and repeat steps 2-3.

[0028] In a preferred embodiment of the present invention, the Goertzel algorithm includes:

[0029] Using the properties of the discrete-time Fourier transform (DTFT), the DTFT transforms the discrete-time signal into Functions mapped to a continuous complex plane ,in is the frequency;

[0030] The DTFT definition is:

[0031] ,

[0032] Will Expressed as the difference of its prefix sum:

[0033] ,

[0034] Substituting the above formula into the definition of DTFT, we get:

[0035] ,

[0036] Expand the inner summation:

[0037] ,

[0038] In the second summation Replace with ,get:

[0039] ,

[0040] In the second summation Replace with ,get:

[0041] ,

[0042] In the second summation Replace with ,get:

[0043] ,Will Put it forward and get:

[0044] ,

[0045] The above formula can be regarded as a recursive formula, where Indicates the status of the filter:

[0046] ,

[0047] in, ;

[0048] Finally, the signal is obtained at a specific frequency The amplitude at :

[0049] ,

[0050] The calculation results of all frequency points are recorded as ;

[0051] Where e is a natural constant, n is the number of current data points, and j is a complex unit.

[0052] In a preferred embodiment of the present invention, the operation of condition three includes: using a machine learning algorithm to train a model, using the trained model to identify the characteristics of the real-time waveform, and determining whether it is a pre-earthquake waveform. The specific steps are as follows:

[0053] C1. Collect the detected pre-seismic waves as a data set and pre-process the data set, including denoising, filtering and normalization, to ensure data quality and consistency;

[0054] C2, extract the features of the data set, including time domain features and frequency domain features;

[0055] C3, divide the data set into training set and test set, and use cross-validation method to ensure the generalization ability and robustness of the model;

[0056] C4. Select a suitable machine learning model, train the selected machine learning model using the training set, and adjust the hyperparameters of the machine learning model to optimize performance;

[0057] C5. Use the training set to train the selected machine learning model and use the test set to evaluate the performance of the trained model;

[0058] C6. Adjust and optimize the machine learning model based on the evaluation results to improve its performance;

[0059] C7. After the machine learning model achieves satisfactory performance, it is input into the early warning system to identify the real-time waveform. If the pre-seismic wave is identified, the output is True, and if it is not identified, the output is False.

[0060] In a preferred embodiment of the present invention, the machine learning model adopts a support vector machine model. The support vector machine is a machine learning algorithm based on statistical learning theory, which can use n-dimensional vectors To represent the n quantized features of the object, and use To express The feature vector of an object; Indicates two different categories and uses To express The category to which the object belongs;

[0061] Support vector machine considers that the feature vector is in n-dimensional space points, if we can find a hyperplane that exists in this space to separate the points belonging to the two types of objects, we can use the eigenvector Establish the corresponding category relative to the position of the hyperplane The mapping:

[0062] ,

[0063] Among them, when hour, ,otherwise ;

[0064] for A vector in space, a normal vector to the hyperplane, and The dimensions are the same;

[0065] Express Perform transposition; is a scalar that determines the specific location of the hyperplane; The sign of reflects the position of the eigenvector relative to the hyperplane; the better the hyperplane is determined, and the corresponding The higher the proportion of agreement;

[0066] The search for the best hyperplane is transformed into solving the following conditional extremum problem:

[0067] ,

[0068] In order to use more data information, the restrictions are relaxed and the following conditional extreme value problem is obtained:

[0069] ,

[0070] st ,

[0071] In this extreme value problem, The value can be less than 1. Also adding weight It was then introduced into the process of solving the extreme value, and the SMO algorithm was used to quickly solve this problem and obtain the optimal separating hyperplane.

[0072] The beneficial effects of the present invention are: a forward earthquake prediction algorithm pointed out by the present invention adopts a multi-source collaborative architecture, predicts the trend of earthquake waveforms in the time domain, provides valuable reserved time for early warning and disaster avoidance before the earthquake arrives, and adopts a simplified Fourier transform in the frequency domain, which saves computing resources and time, improves the reliability and stability of the system in the long-term monitoring and operation of the early warning system, can take into account the time domain characteristics and frequency domain changes of earthquake waves, avoids the problems of missing pre-seismic waves, misidentifying background noise and low accuracy of pre-seismic wave identification, improves the accuracy of earthquake prediction, and reduces false alarms and missed alarms. DETAILED DESCRIPTION

[0073] The technical solutions in the embodiments of the present invention are described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0074] Taking C++ programming as an example, the embodiment of the present invention includes: Example 1

[0075] A forward earthquake prediction algorithm includes: an AND operation for the following three condition predictions:

[0076] Condition 1: The amplitude trend of the real-time waveform is increasing;

[0077] Condition 2: The amplitude of the real-time waveform changes at a specific frequency point;

[0078] Condition 3: The real-time waveform is identified as the pre-earthquake waveform;

[0079] If all three of the above conditions are true, an earthquake warning will be triggered;

[0080] The operation of condition one includes: predicting the amplitude trend of the real-time waveform and judging whether it is increasing continuously. The specific steps are as follows:

[0081] A1. Create a stack stack1 of size 3. The stack is an array used to store the three data segments to be processed. Initialize stack1 and set all positions in the stack to None, indicating that the stack is currently empty.

[0082] A2. Set an infinite loop while(ture), read the real-time incoming data segment newData in the loop, and push it into stack1;

[0083] Perform the following operations on newData:

[0084] If stack1 is not full, put newData on the top of stack1 and update the top position of stack1;

[0085] If stack1 is full, first remove the data segment at the bottom of stack1, then put newData at the top of stack1, and update the position of the top of stack1;

[0086] A3. Call CalDataVpp function for newData. CalDataVpp function is defined as follows:

[0087] CalDataVpp is used to process the data in stack1 and judge whether the processing result is true according to the set threshold value, which is recorded as ResultVpp. It accepts two parameters, stack1 and threshold value.

[0088] Pop out the three pieces of data stored in stack1 one by one, and call the calPeakToPeak function to calculate their respective average peak-to-peak values , , , where the calPeakToPeak function called in the CalDataVpp function is defined as follows:

[0089] The calPeakToPeak function accepts a vector type parameter data, which is used to calculate the average peak-to-peak value of data;

[0090] First, initialize maxVal and minVal to the first element of data, define the array Vpp to be empty, and then traverse data to find the peak value: for each data point, if it is larger than the previous and next data points, it is considered a peak value. ;

[0091] After finding the peak, continue to traverse the data to find the valley value corresponding to the peak: starting from the current peak position, traverse the data backward until a data point smaller than the next data point is found, which is considered to be the valley value corresponding to the peak. The corresponding valley value ;

[0092] calculate

[0093] ,

[0094] And store it in the array Vpp, calculate the average value of all elements in Vpp , and returns the calculation result.

[0095] A4. Calculation

[0096] ;

[0097] ;

[0098] judge Is it true and return the result;

[0099] The operation of condition 2 includes: in the frequency domain, monitoring the amplitude of the real-time waveform at a specific frequency point, and determining whether the amplitude has a jump. The specific steps are as follows:

[0100] B1. In the main function, create a stack stack2 with a size of 2 to store the two data segments to be processed;

[0101] B2, enters an infinite loop, reads the real-time data newData in the loop and pushes it into stack2, and calls the function CalDataGoe to process newData and the previous data segment formerData;

[0102] The function CalDataGoe is used to process two pieces of data stored in the stack stack2, and determine whether the conditions are met, and the returned result is recorded as ResultGoe; the function accepts three parameters: stack stack2, sampling rate fs and specific frequency point F;

[0103] Pop out the two data segments in stack2, call the Goertzel function, and calculate the amplitudes of the two data segments at specific frequency points respectively;

[0104] The Goertzel function is used to calculate the amplitude at a specific frequency point. It accepts three parameters, the data segment data, the sampling rate fs and the specific frequency point F. In this embodiment, the specific frequency points are set to 1 Hz, 2 Hz, 3 Hz, 4 Hz, 5 Hz, 6 Hz, 7 Hz, 8 Hz, 9 Hz and 10 Hz.

[0105] First calculate the Goertzel coefficient

[0106] ,

[0107] initialization

[0108] ,

[0109] For each data Do iterative calculations:

[0110] ,

[0111] ,

[0112] ,

[0113] Get the amplitude A of data at frequency F:

[0114] ,

[0115] ,

[0116] ,

[0117] Returns the calculation result ;

[0118] B3. Determine whether the amplitude of newData is greater than formerData and return the result;

[0119] The operation of condition three includes: using machine learning algorithms to train the model, using the trained model to identify the characteristics of the real-time waveform, and determining whether it is a pre-earthquake waveform. The specific steps are as follows:

[0120] C1. Collect the detected pre-seismic waves as a data set and pre-process the data set, including denoising, filtering and normalization, to ensure data quality and consistency;

[0121] C2, extract the features of the data set, including time domain features and frequency domain features;

[0122] C3, divide the data set into training set and test set, and use cross-validation method to ensure the generalization ability and robustness of the model;

[0123] C4. Select a suitable machine learning model, use the training set to train the selected machine learning model, and adjust the hyperparameters of the machine learning model to optimize performance. In this embodiment, the machine learning model adopts a support vector machine model;

[0124] In the main function of the model training code, input the training set samples, which is a two-dimensional vector that stores all samples and the features of each sample;

[0125] Define a one-dimensional vector labels as sample classification, which corresponds one-to-one to the samples in samples. Pre-seismic waves belong to category 1, and other samples belong to category -1.

[0126] Define the svm_problem structure prob to store training data, prob.1 represents the number of training samples, and prob.y is a pointer to a label array for storing sample labels;

[0127] Prepare feature vectors, each of which is an array of svm_node structures; use a loop to create a feature vector for each sample and set the index and value of the feature, with the index starting at 1. Then store the pointer of each feature vector in the nodePtrs vector, and concatenate the data of all feature vectors into the nodes vector. Finally, point the x member of the svm_problem structure to the starting position of the feature vector array.

[0128] Set the parameters of SVM. Set param.svm_type to C_SVC, which means using C_SVC, i.e. multi-category support vector classification; set param.kernel_type to RBF, which means using radial basis function as kernel function.

[0129] C5. Use the training set to train the selected machine learning model and use the test set to evaluate the performance of the trained model;

[0130] Call the svm_train function to train the SVM model and return the trained model, which is stored in the model file svm_model_file.model, that is, a new support vector machine model is trained;

[0131] C6. Adjust and optimize the machine learning model based on the evaluation results to improve its performance;

[0132] C7. After the machine learning model achieves satisfactory performance, it is input into the early warning system to identify the real-time waveform. If the pre-seismic wave is identified, the output is True, and if it is not identified, the output is False. Specifically:

[0133] In the early warning system code, define the extFeatures function to convert real-time data into feature vectors. This function accepts a std::vector <double>The parameter of type vibSignal represents real-time data. In the function, by traversing vibSignal, each data value is converted into a svm_node structure and added to the features vector. The svm_node structure contains the index and value of the feature. The feature index starts from 1, and a termination node is added to the end of the feature vector, adding an svm_node structure with an index of -1.

[0134] In the main function of the early warning system, call the svm_load_model function to load the trained model file svm_model_file.model and assign it to the structure model;

[0135] In an infinite loop, the extFeatures function is called to extract the features of real-time data and convert the real-time data into a feature vector;

[0136] Call the svm_predict function to use the trained model model to predict the extracted feature vector and get the predicted label predictLabel;

[0137] The prediction result is recorded as ResultSVM. If predictLabel=1, ResultSVM is True;

[0138] Earthquake early warning condition judgment, whether the above three conditions are all true:

[0139] Declare a bool type warning flag variable alarm in the main function and set it to false;

[0140] Add a judgment statement in the infinite loop. If "ResultVpp&&ResultGoe&&ResultSVM" is true, the alarm is set to true, and the loop is jumped out to trigger the early warning program.

[0141] Example 2

[0142] On the basis of Example 1, the user adjusts the operation mechanism according to actual needs, takes some results in ResultVpp, ResultGoe and ResultSVM as true, and triggers the early warning program.

[0143] In summary, the forward earthquake prediction algorithm pointed out by the present invention adopts a multi-source collaborative architecture, which improves the accuracy of earthquake prediction and can reduce the problems of false alarms and missed alarms.

[0144] The above descriptions are merely embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention specification, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.< / double>

Claims

1. A forward earthquake prediction algorithm, characterized in that: include: The AND or OR operation of the following three conditions predictions: Condition 1: The amplitude trend of the real-time waveform is increasing; Condition 2: The amplitude of the real-time waveform changes at a specific frequency point; Condition 3: The real-time waveform is identified as the pre-earthquake waveform; When performing "AND operation", if all three conditions above are true, an earthquake warning will be triggered; When performing "OR operation", the operation mechanism is adjusted according to actual needs, and earthquake warning is triggered when some of the conditions are true; The operation of condition one includes: predicting the amplitude trend of the real-time waveform and judging whether it is increasing continuously. The specific steps are as follows: A1. Create an array of size 3 to store the data segments to be processed; A2. Receive real-time input data segment D N , perform the following operations on the most recently input data segment: If the array is not full, put the data segment at the top of the array and update the top position of the array; If the array is full, first remove the data segment D at the bottom of the array N-2 , D N Put it at the top of the array and update the position of the top of the array; A3. Calculate data segment D N-2 , D N-1 and D N The peak-to-peak value of N-2 、Vpp N-1 and Vpp N ; Where Vpp N (i) is data segment D N The peak-to-peak value of T is D N The total number of peak-to-peak values; Vpp N (i)=D Npeak (i)-D Nvalley (i), D Npeak (i) and D Nvalley (i) respectively for data segment D N Each peak value and its corresponding nearest valley value, Vpp can be calculated by the same logic N-2 and Vpp N-1 ; A4. If Vpp N >Vpp N-1 >Vpp N-2 , the output is True, otherwise the output is False; A5, receive the next data segment and repeat steps A2 to A4; The operation of condition 2 includes: in the frequency domain, monitoring the amplitude of the real-time waveform at a specific frequency point, and judging whether the amplitude jumps. Assuming the specific frequency point is fi, the specific steps are as follows: B1. Create a stack of size 2, which is an array used to store two data segments to be processed; B2. Receive real-time input data segment D N , for D N and D N-1 The Goertzel algorithm is used to calculate the amplitude of a specific frequency and store the calculation results separately. N The frequency calculation result is recorded as D N-1 The calculation result is recorded as B3, if G N All elements in The output is True, and False in other cases; B4, receive the next data segment and repeat steps B2 to B3; The Goertzel algorithm consists of: Using the properties of discrete-time Fourier transform, DTFT maps the discrete-time signal x(n) to a function X(e) on a continuous complex plane by recursively calculating the amplitude of the signal at a specific frequency. jω ), where ω is the frequency; The DTFT definition is: Express x(n) as the difference of its prefix sums: Substituting the above formula into the definition of DTFT, we get: Expand the inner summation: Replacing n in the second summation with n+1 yields: Replacing k in the second summation with k+1 yields: Replacing n in the second summation with n+1 yields: E -jωn Put it forward and get: The above formula can be regarded as a recursive formula, where Q n Indicates the status of the filter: Q n =e -jωN / 2 Q n-1 +x(n)-e -jωN(n-1) / N x(n-N), Among them, Q -1 =Q -2 =0; Finally, the amplitude A of the signal at a specific frequency f is obtained: A=Q N-1 2 +Q N-2 2 -2cos(2πf)Q N-1 Q N-2 , The calculated results of all frequency points are recorded as G N ; Where e is a natural constant, n is the number of current data points, and j is a complex unit; The operation of condition three includes: using machine learning algorithms to train the model, using the trained model to identify the characteristics of the real-time waveform, and determining whether it is a pre-earthquake waveform. The specific steps are as follows: C1. Collect the detected pre-seismic waves as a data set and pre-process the data set, including denoising, filtering and normalization, to ensure data quality and consistency; C2, extract the features of the data set, including time domain features and frequency domain features; C3, divide the data set into training set and test set, and use cross-validation method to ensure the generalization ability and robustness of the model; C4. Select a machine learning model, train the selected machine learning model using the training set, and adjust the hyperparameters of the machine learning model to optimize performance; C5. Use the training set to train the selected machine learning model and use the test set to evaluate the performance of the trained model; C6. Adjust and optimize the machine learning model based on the evaluation results to improve its performance; C7. After the machine learning model achieves satisfactory performance, it is input into the early warning system to identify the real-time waveform. If the pre-seismic wave is identified, the output is True; if not, the output is False. The machine learning model adopts a support vector machine model. The support vector machine is a machine learning algorithm based on statistical learning theory. It uses a g-dimensional vector r to represent the g quantized features of the object and uses {x i |i=1,2,…m} to represent the feature vectors of m objects; use y∈{-1,1} to represent two different categories, and use {y i |i=1,2,…m} to indicate the category to which m objects belong; The support vector machine considers the feature vector to be m points in the g-dimensional space. If a hyperplane existing in this space can be found to separate the points belonging to the two types of objects, then a mapping to the corresponding category y is established based on the position of the feature vector x relative to the hyperplane: k α，δ (x)=p(w T ·x+δ), Among them, when z=w T When x+δ>0, p(z)=1, otherwise p(z)=-1; α is a vector in the X space, a normal vector of the hyperplane, and has the same dimension as x; T represents the transposition of α; δ is a scalar that determines the specific position of the hyperplane; α T The sign of x+δ reflects the position of the eigenvector relative to the hyperplane. The search for the best hyperplane is transformed into solving the following conditional extremum problem: In order to use more data information, the restrictions are relaxed and the following conditional extreme value problem is obtained: st y i ·(a T ·x i +δ)≥1-ξ i ,x i ≥0, i=1,…, m, In this extreme value problem, y i (α T ·x i +δ) is less than 1, and the degree to which this is less than 1 is ξ i After adding the weight C, it is also introduced into the process of solving the extreme value. The SMO algorithm is used to quickly solve this problem and obtain the optimal separating hyperplane.

Images