Continuous prediction method of earthquake early warning magnitude for high-speed railway based on support vector machine

Purpose – Usingthestrongmotiondataof K -netinJapan,thecontinuousmagnitudepredictionmethodbased on support vector machine (SVM) was studied. Design/methodology/approach – In the range of 0.5 – 10.0 s after the P -wave arrival, the prediction time window was established at an interval of 0.5 s. 12 P -wave characteristic parameterswere selectedas the model input parameters to construct the earthquake early warning (EEW) magnitude prediction model (SVM-HRM) for high-speed railway based on SVM. Findings – The magnitude prediction results of the SVM-HRM model were compared with the traditional magnitude prediction model and the high-speed railway EEW current norm. Results show that at the 3.0 s time window,themagnitudepredictionerroroftheSVM-HRMmodelisobviouslysmallerthanthatofthetraditional τ c method and P d method. The overestimation of small earthquakes is obviously improved, and the construction of the model is not affected by epicenter distance, so it has generalization performance. For earthquake events with themagnituderangeof3 – 5,thesinglestationrealizationrateoftheSVM-HRMmodelreaches95%at0.5safterthe arrivalof P -wave,whichisbetterthanthefirstalarmrealizationratenormrequiredby “ TheTestMethodofEEW andMonitoringSystemforHigh-SpeedRailway. ” Forearthquakeeventswithmagnitudesrangingfrom3to5,5to 7 and 7 to 8, the single station realization rate of the SVM-HRM model is at 0.5 s, 1.5 s and 0.5 s after the P -wave arrival, respectively, which is better than the realization rate norm of multiple stations. At the latest, 1.5 s after the P -wave arrival, the SVM-HRM model can issue the first earthquake alarm that meets the norm of magnitude prediction realization rate, which meets the accuracy and continuity requirements of high-speed railway EEW magnitude prediction.


Introduction
Earthquake is a sudden natural disaster that is extremely harmful to railway traffic safety. Even for an earthquake of small magnitude, its impact on subgrade, track and bridge may cause a grave accident that can endanger the life safety of passengers (Guo, 2016;Sun, Wang, Dai, Mu, & Liu, 2007;Wang, Dai, Xi, & Wei, 2018;Yang, 2018). Earthquake early warning system of high-speed railway is one of the effective means to reduce heavy losses caused by earthquake. If an alarm or emergency treatment can be implemented for high-speed railway a few tens or even seconds in advance before the destructive vibration of earthquake occurs, the probability of life and property losses of passengers can be significantly reduced (Song, 2013). With the rapid development of high-speed railway in China, it is extremely urgent to build an earthquake early warning system for high-speed railway and improve the earthquake early warning capacity. Meanwhile, relevant research will also provide supportive reference for the design and construction of Sichuan-Tibet Railway crossing multiple seismic fault zones.
Magnitude prediction is an important link in earthquake early warning of high-speed railway. All the release of earthquake early warning information, the judgment of earthquakeaffected range and the determination of emergency train control range depend on the results of magnitude prediction. For magnitude prediction, the linear proportional relation between the initial characteristic parameters of P-wave and magnitude is usually used, and a model for predicting magnitude is established based on this proportional relation. In the prediction model, the adopted characteristic parameters are mainly divided into two categories, namely, the periodic parameters (Huang, Lin, & Wu, 2015;Kanamori, 2005;Li & Zhang, 2018;Ma, 2008;Ziv, 2014) and the amplitude parameters (Wu & Zhao, 2006;Wang & Zhao, 2017;Zollo, Lancieri, & Nielsen, 2006). The commonly used magnitude prediction models based on the two types of parameters are τ c method (Kanamori, 2005) and P d method (Wu & Zhao, 2006).
In the earthquake early warning system based on dense network layout of earthquake monitoring station, the magnitude estimation can be carried out by means of weighted average or co-processing of multi-station information. Different from the earthquake industry, due to the linear characteristics of railway, the seismic monitoring stations of high-speed railway can only be arranged in a sparse and linear manner along the railway. Due to the limitation of the number and distribution of stations, the earthquake early warning of highspeed railway can only be subject to a single magnitude estimation mode, so it is necessary to establish a magnitude prediction model with higher accuracy.
For the earthquake early warning of existing high-speed railways, the fixed time window of 3.0 s after the earthquake P-wave arrival is usually used for magnitude prediction, but the commonly used magnitude prediction models, namely τ c method and P d method, only use the single parameter characteristics at the initial stage of earthquake wave, resulting in low generalization capacity of prediction model, great discreteness of prediction results, obvious overestimation of small earthquake and underestimation of large earthquake and the need to screen epicentral distance and signal-to-noise ratio when the linear proportional relation between characteristic parameters and magnitude is being established. With the rapid development of big data analysis technology, the application of machine learning method in the field of earthquake early warning has been gradually expanded. Some scholars and experts have applied machine learning method to the research of earthquake early warning magnitude prediction. For example, Ma (2008) has established the relationship between multiple parameters and magnitude based on artificial neural network, and Reddy and Nair (2013) have established the relationship between wavelet coefficient and magnitude with support vector machine (SVM), which provides a feasible reference for magnitude prediction with multi-parameter input and artificial intelligence methods.
In order to send accurate earthquake early warning and alarm as early as possible, meet the constantly updating requirements of earthquake early warning alarm for high-speed railway with the change of time and achieve the goal of "improving the accuracy of magnitude prediction RS within 3.0 s and increasing the continuity of magnitude prediction after 3.0 s," this paper selects the range of 0.5-10.0 s after the P-wave arrival and establishes a prediction time window at an interval of 0.5 s. Based on the classical SVM method in the field of artificial intelligence machine learning, 12 P-wave characteristic parameters in total, including such 4 categories as amplitude parameter, periodic parameter, energy parameter and derivative parameter, are taken as the input parameters of SVM to construct the SVM-based high-speed railway magnitude prediction model (SVM-HRM).The magnitude prediction results in 3.0 s time window are compared with τ c method and P d method, respectively. Moreover, according to relevant provisions of The Test Method of EEW and Monitoring System for High-Speed Railway in force (China Railway Corporation, 2018), the statistical realization rate of magnitude prediction is compared with the realization rate norm, aiming at providing reference for the construction and improvement of high-speed railway earthquake early warning system in China.

Data and processing
Machine learning algorithm needs to be based on big data statistics. The data recorded by China's existing high-speed railway earthquake monitoring station is limited, while the strong earthquake observation station of the same type as the former is adopted for Japan's K-net strong earthquake observation network, and a large number of high-quality strong earthquake data have been recorded. Therefore, this paper mainly selects the data of Japan's K-net strong earthquake observation network (National Research Institute for Earth Science and Disaster Resilience (NIED), 2018) and the strong motion data based on the following principles (Song, Jiao, Li, & Hou, 2018).
The selected data were processed as follows.
(1) The methods proposed by Ma, Jin, Li, Chen, Liao, & Wei (2013) and Wang and Zhao (2016) were adopted to automatically detect the P-wave arrival and manually check the P-wave arrival based on the acceleration data.
(2) The acceleration data were integrated once to get velocity record, and velocity record was integrated once to obtain the displacement record. The integrated records were filtered by 0.075 Hz Butterworth high-pass filter to eliminate the low-frequency drift caused by the integration.
(3) Considering the numerical distortion caused by the great change of characteristic parameter values and the improvement of training efficiency of the model, the characteristic parameters were calculated at an interval of 0.5 s within the time range of 0.5-10.0 s after the P-wave arrival. The data were normalized in turn, and the normalization method could be expressed as follows Where b g is the normalized result of each characteristic parameter; x is the data corresponding to each characteristic parameter; x max and x min are the maximum and minimum values of corresponding P-wave characteristic parameters, respectively.
After processing, a total of 1,837 earthquakes were screened out, with 19,263 groups and 57,789 pieces of strong motion records, and 922 stations with strong motion data were recorded. The distribution of epicenters and recording stations obtained by screening is shown in Figure 1. In the Figure, the red hollow circles represent the position of epicenter, and the diameter of circle is directly proportional to the magnitude; green triangles represent stations where data are recorded.
The strong motion data obtained by screening were randomly divided into two groups that were not repeated, namely training data set and test data set, in which the training data set accounted for 80% of all data, with a total of 15,410 groups and 46,230 pieces of data; test data set accounted for 20% of all data, with 3,853 groups and 11,559 pieces of data. The relationship between magnitude, epicentral distance and number of records of screened data is shown in Figure 2, in which green dots represent the training data sets for establishing  Relationship between magnitude and epicentral distance and number of records RS SVM-HRM prediction model, and red dots represent test data sets for testing SVM-HRM prediction model.

Model characteristic parameters
High-speed railway earthquake early warning can only be subject to the form of single magnitude estimation, which requires higher accuracy, so it is necessary to integrate multiple parameters when analyzing the magnitude estimation results. In this paper, the parameters were set based on the purpose of obtaining the optimal prediction results, and 12 characteristic parameters including such four categories as amplitude parameter, periodic parameter, energy parameter and derivative parameter were selected as the input parameters of SVM-HRM prediction model, in which the amplitude parameter, energy parameter and derivative parameter were normalized and uniformly corrected to the reference focal distance of 10 km (Peng, Yang, Xue, Chen, & Zhu, 2013;Yamada, Heaton, & Beck, 2007;Zollo, Lancieri, & Nielsen, 2006). Each characteristic parameter is defined as follows.

Amplitude parameter
It includes peak displacement P d (Wu & Zhao, 2006), peak velocity P v and peak acceleration P a , and the calculation formulas are as follows Where 0 is the arrival time of P-wave; τ 0 is the time window length after the P-wave arrival; s(t) is the vertical displacement time history; v(t) is the vertical velocity time history; a(t) is the vertical acceleration time history.

Periodic parameter
It includes characteristic period τ c (Kanamori, 2005), peak ratio T va (Ma, 2008) and structural parameter P P (Huang, Lin, & Wu, 2015). The calculation formulas are as follows Where r is the intermediate variable obtained according to Parseval's theorem.
3.3 Energy parameter It includes velocity square integral I V2 (Festa, Zollo, & Lancieri, 2008), cumulative absolute velocity C AV (Reed & Kassawara, 1990;Song, 2013) and cumulative energy change rate D I (Nakamura, 2003) . The calculation formulas are as follows Where a 3(t) is the acceleration synthesized in three components.

Derivative parameter
It includes vertical cumulative absolute displacement C ad , vertical cumulative absolute velocity C av and vertical cumulative absolute acceleration C aa . The calculation formulas are as follows

Prediction model algorithm
SVM is a machine learning method based on statistical learning in the field of artificial intelligence, which can use multiple parameters for mode classification and nonlinear regression analysis (Saunders, Stitson, Weston, Holloway, Bottou, Scholkopf, & Smola, 2002). In this paper, the training data set obtained from the above screening is used to establish the SVM model algorithm based on the obtained fitting function. The key of the algorithm is to establish Gaussian kernel parameters and training parameters.

Support vector machine algorithm
Definition of parameters: f(X) is the predicted magnitude; W is the weight vector of each characteristic parameter; X is the vector composed of characteristic parameters; b is the intercept; n is the number of training set data, and i is the ith data; y i is the magnitude corresponding to the characteristic parameter; E is the tolerance error. The linear regression function f (X) 5 W T · X þ b is used to fit 12 parameters in 4 categories calculated based on the above training data set. Assuming that after fitting, all sample data can be expressed by linear function number f(X) in the range of [ÀE, E] (Li, 2011), and its mathematical expression is On this basis, the SVM model is established through the following steps. 4.1.1 Problem transformation. Optimize the distance between data points and linear regression functions, and transform the problem of calculating distance into an extremum problem to obtain Where kW k is the module of the weight vector.
4.1.2 Introduction of penalty parameters. Different from the least square fitting method, the SVM algorithm allows a certain range of fitting error, so the penalty parameter C can be introduced to indicate the penalty degree of the sample when the allowable error is exceeded, and the objective function (extremum calculation formula) of Formula (16) is transformed into In which, Where S i , S* i are slack variables, which are also recorded as loss functions, that is, training errors on the upper and lower sides of the regression function.
The relationship between linear regression function and loss function of SVM is shown in Figure 3. The solid line in the figure represents the linear regression function f (X) 5 W T $X þ b, and the dashed line range [ÀE, E] is called the tolerance range, representing the hyperplane of the SVM. If the error of sample points (hollow circles) falling within the dashed lines is negligible, the sample points (red solid circles) falling on these two dashed lines are recorded as support vectors; the sample points (hollow circles) falling outside the range of dashed line, that is, exceeding the range of [ÀE, E], can be recorded as loss function S. 4.1.3 Transformation of objective function. In order to solve the extremum problem in step 2, the above objective function is transformed into Lagrangian function, so as to solve the constrained extremum problem. Thus, Formula (17) can be transformed into L À W ; b; S; S * ; α; β; γ Where α, β and γ are Lagrangian factors.

Solving Lagrangian function. Derive Formula (18) and get
Where α i and β i are not 0 simultaneously, and α i > 0; K(Xi, X) is the kernel function of SVM.
The linear regression function of SVM can be obtained as follows through the above four steps

Gaussian kernel parameters of support vector machine
The most important difference between the SVM and linear regression is that the former maps sample data to high-dimensional space through kernel function operation, that is, based on the linear combination of nonlinear functions, and its network structure is shown in There are many forms of kernel functions. In this paper, the kernel function refers to the Gaussian kernel function, which has a wide application range and strong differential ability and can be used to better extract the local features of sample data (Steinwart, Hush, & Scovel, 2009). The calculation method is as follows Where, λ is the Gaussian kernel parameter.

Training parameters of support vector machine
Appropriate training parameters shall be selected for a reasonable SVM model, so as to minimize the error between predicted results and true values. The error-related parameters in SVM include penalty parameter C, tolerance error E and Gaussian kernel parameter λ. The calculation in this paper is based on the empirical calculation given by Ma (2002, 2004), which is as follows:

Continuous prediction of magnitude
Where u is the average value of the output results of training set data; γ is the standard deviation of the output results of training set data; n is the number of data in the training set; η is the standard deviation of the error between the predicted value and the true value obtained by training; q is the range of input characteristic parameter values; m is the number of input characteristic parameters. According to the calculation method of training parameters given in Formulas (23)-(25), the penalty parameter C, tolerance error E and Gaussian kernel parameter λ under different time windows were determined after six crossing validations, and then the SVM-HRM prediction model was established under different prediction time windows.

Comparison of magnitude prediction results of SVM-HRM prediction model
Based on the test data set, the magnitude prediction results of SVM-HRM prediction model are compared as follows.
(1) Under the time window of 3.0 s, the magnitude prediction results of SVM-HRM prediction model were compared with the traditional magnitude prediction models of τ c method and P d method.
(2) According to relevant provisions of The Test Method of EEW and Monitoring System for High-Speed Railway in force, the single-station realization rate of magnitude prediction of SVM-HRM prediction model was calculated and compared with the magnitude prediction realization rate norm required in The Test Method of EEW and Monitoring System for High-Speed Railway.

Comparison with traditional magnitude prediction model
The predicted magnitude of SVM-HRM prediction model is defined as the predicted value, and the cataloged magnitude of earthquake events is defined as the true value. If the difference between the predicted value and the true value is error ω i , the standard deviation σ of ω i is Where N is the number of data in the test set, and k is the kth data. Under the time window of 3.0 s after the arrival of P-wave, the magnitude prediction results of SVM-HRM prediction model and traditional magnitude prediction models of τ c method and P d method were compared, and the standard deviation σ value of the three models was calculated according to Formula (26). The results are as shown in Figure 5. In the figure, the black solid line represents the 1:1 linear proportional relationship between the predicted value and the true value, and the red dashed line represents the 1x standard deviation σ of the error between the predicted value and the true value.
According to the calculation and Figure 5, the magnitude prediction error standard deviation of τ c method, P d method and SVM-HRM prediction model is 1.64, 0.43 and 0.30 magnitude units, respectively, that is, the error standard deviation of SVM-HRM prediction model is much smaller than that of τ c method, and also smaller than that of P d method.
It should be noted that τ c method significantly overestimates the magnitude prediction results of earthquakes below magnitude 5, that is, there is a phenomenon of "small earthquake overestimation," because τ c method needs to screen the epicentral distance and signal-to-noise ratio of earthquake data before a suitable prediction model is constructed.
There is no such step for SVM-HRM prediction model constructed in this paper, which greatly improves the universality of the model. At the same time, through the comparison of the RS Figure 5. Comparison of magnitude prediction results between τ c method, P d method and SVM-HRM Continuous prediction of magnitude magnitude prediction results of P d method, it can be clearly found that the "small earthquake overestimation" of SVM-HRM prediction model has also been improved.
Obviously, through the comparison of τ c method with P d method, the accuracy of magnitude prediction obtained by SVM-HRM prediction model is significantly improved.
5.2 Comparison with the realization rate norm of high-speed railway magnitude prediction The change of magnitude prediction error standard deviation of SVM-HRM prediction model based on the training data set and that based on the test data set with the prediction time window is shown in Figure 6. In the figure, the abscissa is the time window, and 0 s is the arrival time of P-wave. It can be seen from Figure 6 that the two curves tend to coincide, and the maximum difference of error standard deviation corresponding to the same time window does not exceed 0.02, which indicates that the SVM-HRM prediction model constructed in all prediction time windows has strong generalization ability after the P-wave arrival, that is, the SVM-HRM prediction model has strong adaptability under fresh data samples, and the trained network can also give appropriate output of data other than the training data set with the same rule; with the increase of time window, the error standard deviation of magnitude prediction decreases significantly, which indicates that SVM-HRM prediction model has the continuity of magnitude prediction, and with the gradual increase of time window after the P-wave arrival, the accuracy of magnitude prediction increases significantly.
The Test Method of EEW and Monitoring System for High-Speed Railway stipulates and puts forward requirements for the predicted magnitude realization rate in the first alarm (the Figure 6. Change of error standard deviation of SVM-HRM magnitude prediction with time window RS alarm issued for the first time) of earthquake early warning (China Railway Corporation, 2018), and takes the realization rate as an index parameter to measure the accuracy of magnitude prediction of high-speed railway earthquake early warning system. The realization rate is defined according to the percentage of the number of predicted magnitude errors with their absolute values less than or equal to 1, and the calculation formula is Where r is the realization rate; h is the number of predicted magnitude error ω i ≤ 1 in the test data; H is the total number of test data. See Table 1 for the realization rate norm of magnitude prediction specified in The Test Method of EEW and Monitoring System for High-Speed Railway. In this paper, the change rule of realization rate with prediction time window was studied, aiming at analyzing how long the SVM-HRM prediction model can give the accurate and reliable first alarmed magnitude after the P-wave arrival.
It should be pointed out that in high-speed railway earthquake early warning, the magnitude prediction results of multiple stations need to be weighted and averaged according to the prediction results of single station. Therefore, in the comparative analysis of the realization rate of multiple stations, as long as the statistical results of the realization rate of single station in different magnitude ranges are better than the realization rate norm of multiple stations, the realization rate of multiple stations is deemed as better than the norm. Therefore, it is only necessary to compare the realization rate of single station in different magnitude ranges with the corresponding realization rate norm of magnitude range in The Test Method of EEW and Monitoring System for High-Speed Railway.
Under different time windows, the calculation results of single station realization rate of magnitude prediction by SVM-HRM prediction model and the comparison with the realization rate norm specified in The Test Method of EEW and Monitoring System for High-Speed Railway are shown in Figure 7. The dashed lines in the figure are the single station (red dashed line) and multiple stations (blue dashed line) realization rate norm bureau specified in The Test Method of EEW and Monitoring System for High-Speed Railway in force (China Railway Corporation, 2018).
It can be seen from Figure 7a that for earthquake events with magnitudes ranging from 3 to 8, the SVM-HRM prediction model obtains that the realization rate of the predicted magnitude of a single station reaches 95% at 0.5 s after the P-wave arrival, which is better than the realization rate norm of a single station specified in The Test Method of EEW and Monitoring System for High-Speed Railway. With the increase of prediction time window, the realization rate gradually increases, indicating that the accuracy of magnitude prediction continues to increase; when the time window reaches 2.0 s, the magnitude prediction realization rate of SVM-HRM prediction model is close to 100%.

Magnitude range/ magnitude
The realization rate about the deviation not more than 1 between the first alarmed predicted magnitude and the actual magnitude Single station 3-8 Not less than 50% Multiple stations 3-5 Not less than 30% 5-7 Not less than 90% 7-8 Not less than 60% It can be seen from Figure 7b that for earthquake events with magnitudes ranging from 3 to 5, the realization rate of SVM-HRM prediction model at 0.5 s after the P-wave arrival reaches 98%, which is better than the multi-station realization rate norm specified in The Test Method of EEW and Monitoring System for High-Speed Railway. With the increase of prediction time window, the realization rate gradually increases, indicating that the accuracy of magnitude prediction continues to increase; when the time window reaches 1.0 s, the magnitude prediction realization rate of SVM-HRM prediction model is close to 100%. It can be seen from Figure 7c that for earthquake events with magnitudes ranging from 5 to 7, the realization rate of SVM-HRM prediction model reaches 92% at 1.5 s after the P-wave arrival, which is better than the multi-station realization rate norm specified in The Test Method of EEW and Monitoring System for High-Speed Railway. With the increase of prediction time window, the realization rate gradually increases, indicating that the accuracy of magnitude prediction continues to increase; when the time window reaches 4.5 s, the magnitude prediction realization rate of SVM-HRM prediction model is close to 100%.
It can be seen from Figure 7d that for earthquake events with magnitudes ranging from 7 to 8, the realization rate of SVM-HRM prediction model at 0.5 s after P-wave arrival is 67%, which Comparison between single-station realization rate of SVM-HRM prediction model and the current realization rate norm under different time windows RS is better than the multi-station realization rate norm specified in The Test Method of EEW and Monitoring System for High-Speed Railway. With the increase of prediction time window, the realization rate gradually increases, indicating that the accuracy of magnitude prediction continues to increase; when the time window reaches 2.5 s, the magnitude prediction realization rate of SVM-HRM prediction model reaches 95% and enters the platform stage. When the time window reaches 6.0 s, the realization rate continues to increase, and when the time window reaches 7.0 s, the magnitude prediction realization rate reaches 100%.
Through the above analysis, it can be found that the accuracy and continuity of magnitude prediction of SVM-HRM prediction model established in this paper have been greatly improved, which meet the requirements of relevant provisions on magnitude prediction in The Test Method of EEW and Monitoring System for High-Speed Railway, and can be used for magnitude prediction of high-speed railway earthquake early warning system.

Conclusion
(1) Under the 3.0 s time window after the P-wave arrival, compared with τ c method and P d method for traditional earthquake early warning magnitude prediction, the magnitude prediction error of SVM-HRM prediction model is obviously reduced; the phenomenon of overestimation of small earthquakes is substantially improved; and the accuracy is greatly increased.
(2) The training data set of SVM-HRM prediction model cannot be screened for epicentral distance and signal-to-noise ratio, which indicates that the model has strong universality; the error standard deviation of magnitude prediction results between the training data set and the test data set of the model tends to be consistent, and decreases significantly with the increase of time window, which indicates that the model has strong generalization performance and the continuity of magnitude prediction is greatly improved.
(3) After comparing the magnitude prediction results of SVM-HRM prediction model with the first alarmed magnitude prediction realization rate norm required in The Test Method of EEW and Monitoring System for High-Speed Railway, it can be found that the single-station magnitude prediction realization rate of SVM-HRM prediction model reaches 95% at 0.5 s after the P-wave arrival, which is better than the first alarmed magnitude prediction realization rate norm required by The Test Method of EEW and Monitoring System for High-Speed Railway. For the earthquake events with the magnitudes ranging from 3 to 5, 5 to 7 and 7 to 8, the single-station magnitude prediction realization rate is at 0.5 s, 1.5 s and 0.5 s after the P-wave arrival, respectively, which is better than the multi-station realization rate norm. This indicates that with the SVM-HRM prediction model constructed in this paper, the first alarm of earthquake early warning can be issued at the latest 1.5 s after the P-wave arrival according to the magnitude prediction realization rate norm in The Test Method of EEW and Monitoring System for High-Speed Railway, and the accuracy is greatly improved.