LEAST SQUARE SUPPORT VECTOR MACHINE FOR DETECTION OF TECSEISMO- IONOSPHERIC ANOMALIES ASSOCIATED WITH THE POWERFUL NEPAL EARTHQUAKE (M<sub>w</sub> = 7.5) OF 25 APRIL 2015

Due to the irrepalable devastations of strong earthquakes, accurate anomaly detection in time series of different precursors for creating a trustworthy early warning system has brought new challenges. In this paper the predictability of Least Square Support Vector Machine (LSSVM) has been investigated by forecasting the GPS-TEC (Total Electron Content) variations around the time and location of Nepal earthquake. In 77 km NW of Kathmandu in Nepal (28.147 N, 84.708 E, depth=15.0 km) a powerful earthquake of Mw=7.8 took place at 06:11:26 UTC on April 25, 2015. For comparing purpose, other two methods including Median and ANN (Artificial Neural Network) have been implemented. All implemented algorithms indicate on striking TEC anomalies 2 days prior to the main shock. Results reveal that LSSVM method is promising for TEC sesimoionospheric anomalies detection.


INTRODUCTION
The pre-seismic disturbances occurring in lithosphere, atmosphere and ionosphere in absence of significant solar and geomagnetic perturbations are usually considered as earthquake precursors.
The ionospheric anomalies usually take place in the D, E and F layers, and may be observed 1 to 10 days prior to the earthquake and continue a few days after it. Many papers and reports have been published on satellite observations of the ionospheric plasma, the flux of charged particles, the DC electric field, electromagnetic waves and geomagnetic field associated with seismic activity (Parrot, 1995;Liu et al., 2004;Hayakawa and Molchanov, 2002;Pulinets and Boyarchuk, 2004;Akhoondzadeh, 2011). The preearthquake disturbances usually affect the TEC (Total Electron Content). In this study, GPS-TEC data have been downloaded via NASA Jet Propulsion Laboratory (JPL) website. Global Ionospheric Map (GIM) data constructed from a 5° × 2.5° (Longitude, Latitude) grid with a time resolution of 2 hours.
The accurate anomaly detection in time series of earthquake precursors is regarded as one of the most challenging tasks since the behaviors of precursors is complicated, dynamic and nonlinear. In addition, the affected factors in ionospheric precursors such as solar and geomagnetic indices intensify the complexity of the ionospheric precursors.
As a modified version of Support Vector Machine (SVM), Least Square SVM (LSSVM) was proposed by Suykens and Vandewalle in 1999. It retains principle of the Structural Risk Minimization (SRM) and has important improvement of calculating speed with traditional SVMs (Wang and Shang, 2014). It is because of changing inequality constraints into equations and takes a squared loss function. Therefore, LSSVM solves a system of equations instead of a quadratic programming problem. Due to the mentioned advantages of LSSVM, it has been implemented as a classification model (Wang and Shang, 2014). In this paper the predictability of LSSVM has been investigated by predicting the GPS-TEC variations around the time and location of Nepal earthquake.

LSSVM METHOD
SVM training is a time consuming process specially when analyzing huge dataset. For this purpose, LSSVM is proposed by Suykens and Vandewalle in 1999 to overcome these shortcomings.
Unlike the classical neural networks approach, SVM formulates the statistical learning problem as a quadratic programming with linear constraints, by the use of nonlinear kernels, high generalization ability, and sparseness of solution. However, for large-scale problems, the optimization process of SVM has high computational complexity, due to the high-dimensional matrix involved in the quadratic programming whose size is directly proportional to the training sample size (Zhou et al., 2011).
LS-SVM needs significantly less training effort than the standard SVM as a result of the model simplification.
The basic principle of SVM regression is to estimate the output variable y from Where w and b are the weight vector and bias term respectively (Zhou et al., 2011).
= . The first part of cost function regularizes weight sizes and penalizes large weights. Therefore, the weights tend to converge to similar values in that large weights cause excessive variance and hence deteriorate the generalization ability of LSSVM (Suykens and Vandewalle, 1999;Zhou et al., 2011).
The second part of Eq. (2) considers the regression error of all training data. The regularization parameter c controls the trade-off between the bias and variance of LS-SVM model. Note that the LSSVM model has equality constraints as shown in Eq. (3), rather than the inequality constraints with slack variables used in the standard SVM model. Moreover, as shown in Eq. (2), a squared loss function is considered in the objective function of LS-SVM model, while the standard SVM model has a linear combination of slack variables in its objective function. These two modifications simplify the quadratic optimization problem for the standard SVM to be linear for LSSVM (Suykens and Vandewalle, 1999;Zhou et al., 2011).
where i λ are the Lagrange multipliers. By the Karush-Kuhn-Tucker Theorem, the conditions of optimality are (Suykens and Vandewalle, 1999;Zhou et al., 2011) 0 Thus, b and λ can be solved from the following set of linear equations after eliminating w and e, is the kernel function (Suykens and Vandewalle, 1999;Zhou et al., 2011).
As a result, given vectors x and x i , the LS-SVM regression model for estimating y in Eq. (1) becomes  (7) In this study Gaussian kernel function has been used. ) where . denotes the 2-norm, and σ is a constant determining the width of Gaussian kernel.
To implement the LSSVM method, training and testing data were initially set respectively to 60% and 40% of all TEC data. The input patterns in the LSSVM method are, At each step, using the training data, the LSSVM method is implemented and the prediction error (PE) can be written as: Where, i x and i x are the actual value and the output from the LSSVM method, respectively.
Finally, the TEC value is predicted and then is compared to the true value in testing set. In the case of testing process, if the value of DX i (i.e. the difference between the actual value X i and the predicted value i X ) is outside the predefined bounds σ µ × ± k , ( µ and σ are the mean and the standard deviation of DX i values) the anomaly is detected.
In this study, geomagnetic and solar indices (i.e. K p , A p , D st and F10.7) were used to distinguish seismic anomalies from the other anomalies related to the geomagnetic and solar activities. Figure 1 illustrates the variations of K p , A p , D st and F10.7 indices, during the period of 01 March to 30 April 2015. An asterisk indicates the earthquake time. The X-axis represents the days relative to the earthquake day. The Y-axis represents the universal time coordinate.
Figure 2(a) shows TEC variations during the period of 01 March to 30 April 2015. To implement the Median method, Dx which will be called DTEC here, is calculated using Eq. (11).
Where x , high x , low x , M , IQR and Dx are the parameter value, higher bound, lower bound, median value, Interquartile range and differential of x, respectively. According to this, if the absolute value of Dx would be greater than k, ( k Dx > ), the behavior of the relevant parameter (x) is regarded as anomalous. Figure 2 Ap < and F10.7 <150, respectively, are jointly used using AND operator to construct the anomaly map (Figure 2(d)). The TEC value exceeds the higher bound ( . 1 ), 2 days before the earthquake time at 14:00, 16:00 and 18:00 UTC with the values of 18.68 %, 15.97% and 0.32% of the higher bound, respectively. It is seen that the other detected anomalies in Figure 2(c) have been masked by high geomagnetic activities (Figure 2(d)). In figure 2(d) increases (08.40% and 04.13%) in TEC are clearly observed at 16:00 and 18:00 UTC on earthquake date (Table 2).  To implement the Artificial Neural Network (ANN) method (Akhoondzadeh, 2013), training data were set to 60% of all data. Using the training data, the pattern vectors in feature space are constructed. In the case of testing process, if the difference value PE i between the actual value X i and the predicted value i X , is outside the pre-defined bounds σ µ × ± 0 . 2 , ( µ and σ are the mean and the standard deviation of PE i values) the anomaly is detected.
Red and green curves in figures 3(a) through (l) represent the observed and the predicted TEC values using the ANN method, respectively during the days selected as training and testing set. These figures indicate that the ANN method is a good estimator for non linear time series such as TEC variations. It can be seen that the ANN method efficiently predicts values during the time of testing.
Figures 4(a) through (l) represent the differences between the observed and the predicted TEC values during the testing data using the ANN method. Figure 5(a) is a representation of the differences values during the testing set. Figure 5 with the values of 8.14%, 2 days before the earthquake at 14:00 UTC, and also 4 days after the earthquake time at 22:00 UTC with the value of 38.49% (Table 2).
To implement the LSSVM method, training data were set to 60% of all data. Using the training data, the pattern vectors in feature space are constructed. In the case of testing process, if the difference value PE i between the actual value X i and the predicted value i X , is outside the pre-defined bounds σ µ × ± 0 . 2 , ( µ and σ are the mean and the standard deviation of PE i values) the anomaly is detected.
Red and green curves in figures 6(a) through (l) represent the observed and the predicted TEC values using the LSSVM method, respectively during the days selected as training and testing set.
Figures 7(a) through (l) represent the differences between the observed and the predicted TEC values during the testing data using the LSSVM method. Figure 8(a) is a representation of the differences values between the observed and predicted values during the testing set. Figure  8 . 3. a) through l) Variations of the observed (red curve) and predicted (green curve) TEC values obtained from ANN method during at different universal times. The X-axis represents the day relative to the Nepal earthquake day. ). The green horizontal line indicates the mean value ( µ ). The X-axis represents the day relative to the Nepal earthquake day. Fig. 5. a) Differences between the observed and the predicted values of TEC obtained from ANN method. b) DTEC variations. c) Detected anomalies using ANN method without considering the geomagnetic indices. d) Detected anomalies using ANN method with considering the geomagnetic indices. The X-axis represents the days relative to the Nepal earthquake day. The Y-axis represents the universal time coordinate. Fig. 6. a) through l) Variations of the observed (red curve) and predicted (green curve) TEC values obtained from LSSVM method during at different universal times. The X-axis represents the day relative to the Nepal earthquake day. ). The green horizontal line indicates the mean value ( µ ). The X-axis represents the day relative to the Nepal earthquake day. Fig. 8. a) Differences between the observed and the predicted values of TEC obtained from LSSVM method. b) DTEC variations. c) Detected anomalies using LSSVM method without considering the geomagnetic indices. d) Detected anomalies using LSSVM method with considering the geomagnetic indices. The X-axis represents the days relative to the Nepal earthquake day. The Yaxis represents the universal time coordinate. 20:00 UTC. Other post-seismic anomalies could be associated with the strong aftershocks (Table 1).

CONCLUSIONS
This paper attempts to acknowledge the ability of LSSVM as a good predictor to forecast the GPS-TEC variations around the time and location of Nepal earthquake. Also two classical and intelligent methods including Median and ANN have been implemented. Using all applied methods prominent TEC anomalies are observed 2 days prior to the main shock. Results reveal that LSSVM method is promising for TEC sesimo-ionospheric anomalies detection.