Seismic Intelligence: Machine Learning Models for Earthquake Magnitude Estimation

In the latest literature, different ML methods for earthquake magnitude prediction have been investigated and shows their own strengths and weaknesses. A remarkable work used Random Forest models on LA, Istanbul, and San Diego seismic data with high sensitivity rates of 91.65%, 98.53% and 97.97%. Of 16 tested ML and neural network models, Random Forest was the most successful. Nevertheless, as a few drawbacks still exist, it is hard to estimate the exact time and location of an earthquake and to apply in early warning systems [1].

In a different study, Artificial Neural Networks (ANNs) with Levenberg-Marquardt Backpropagation were applied for short-term forecasting (between 1–15 days) at regions covering Japan, Turkey, Greece and Indian Subcontinent. The model showed flexibility and low false alarm rates with the accuracy ranging from 78.36 to 87.65%. Nevertheless, it still relied on the past and its prediction power was not accurate in time [2].

Advanced models including Long Short-Term Memory (LSTM), Support Vector Machines (SVMs), Gradient Boosting, and hybrid stacking have also been used for prediction of magnitude and depth. These methods effectively capture temporal and geological characteristics, which help in preventing overfitting. However, they have the same limitations: computationally expensive and lack generalization to different seismic zones [3].

A comparative study assessed the performance between Random Forest Regressors and Neural Networks with heterogeneous seismic features. Whereas Random Forest performed well only with a variety of data, and Neural Network could recognize intricate trends but it has limitations on overfitting and in terms of computational time [4]. More recently Inverse Boosting Pruning Trees (IBPT) based on satellite infrared and hyperspectral to detect pre-earthquake anomalies were proposed. After evaluation with 1,371 events, IBPT surpassed six benchmarks (including Random Forest and SVM), and was promising for its short-term prediction. Nevertheless, the reliance on satellite information may hinder its scale and usage in resource-deprived areas [5].

In line with these works, Ding and Hu [6] showed how deep learning and computer vision could be employed in monitoring geological hazards, thus indicating that ML techniques can easily be adapted to the detection of hazard. JayaLakshmi et al. [7] employed transformer-based multimodal tweet classification for crisis management, highlighting the promise of ML in real-time events analysis. Ojeda et al. [8] improved seismic vulnerability assessment of confined brick masonry houses with the use of assemblage algorithms, cementing ML’s place within structural risk analysis. Also, Xing [9] developed a real-time monitoring and early warning system based on image analysis that is equivalent to the earthquake early warning application. Shantaf and Kutucu [10] also solved the issue of seismic sensor reliability by the use of ML in detecting signal distortion, Putra and Sunarno [11] presented a Monte Carlo aggregating technique to radon precursor-based earthquake prediction in Java, Indonesia.

Although these studies as shown in Table 1 highlight the promise of ML in seismic prediction, they mostly concentrate on single-model approaches with limitations either in terms of generalization, timing precision or data availability. This study fills these gaps by conducting a stacked ensemble method with CatBoost, LightGBM and Random Forest. By leveraging strengths of various models and selecting feature engineering systematically, the proposed method is expected to enhance accuracy and robustness of magnitude prediction for making better decision of seismic risk assessment and disaster preparedness.

Table 1. Comparison table (prior studies)

The magnitude of an earthquake is the measurement of energy release in an earthquake, and it is measured using a scale that can be correlated to other earthquakes.

4.1 Domain earthquake magnitude scales

The magnitude of an earthquake is a measure of the energy released during a seismic event and is an important parameter for predictive modeling. Mag Type (magType) is of particular interest for machine learning applications as source scales measure seismic energy differently, affecting model performance and transferability.

For larger earthquakes, the Moment Magnitude Scale (Mw) is commonly used due to its uniformity across fault types and depths [12, 13]. On the other hand, Local Magnitude (ML) or Body Wave Magnitude (Mb) may be used for smaller or shallow earthquakes and their use could add to statistical variability associated with regional calibration and station distance [14, 15]. Despite the fact that the Richter Scale was designed to be used for magnitude, these are usually determined using Mw for present seismic analysis [16]. In the context of machine learning, considering magType as a categorical feature is one way to allow models to control for such methodological differences. It is important to at least take it into account, because if not one might obtain biassed predictions, especially when the data contain mixed-scale measures. Careful encoding and normalization of magType is thus crucial to obtain accurate magnitudes.

4.2 Earthquake magnitude ranges

Earthquake Mahufacture sizes can vary enormously and are usually ranged into categories depending upon their size, the amount of energy released, and the potential for damage. The classification of seismic magnitudes with the resulting consequences are presented in Table 3 [17-19].

Table 3. Earthquake magnitude classification

Figure 3 illustrates the implementation workflow for the earthquake‑magnitude prediction project. Key stages comprise data acquisition, pre-processing, transformation, model testing, training, evaluation, stacking/ensemble learning, final deployment, prediction, and feedback.

tu_pian_35.png

Figure 3. Implementation diagram

6.1 Tools and libraries

The project was developed in python 3.12 and run in a Jupyter‑style environment on Google Colab, leveraging libraries for data preprocessing, model development, evaluation, and visualization.

6.2 Key libraries

Table 5 is about Key Libraries used in the study.

Table 5. Key libraries

6.3 EDA

Through EDA, we examined the dataset’s statistical properties, visualized feature interactions, and identified correlations and anomalies. This process revealed patterns that were crucial for feature selection and model design, ultimately shaping a more accurate and reliable predictive framework.

6.3.1 Distribution of earthquake magnitudes

Figure 4 illustrates the distribution of earthquake magnitudes. The x‑axis denotes magnitude, and the y‑axis denotes frequency. The predominance of low‑magnitude events (approximately 1–2) is evident from the tall bars on the left, while the overlaid curve represents a smoothed trend of the distribution.

tu_pian_36.jpg

Figure 4. Distribution of earthquake magnitudes

6.3.2 Distribution of earthquake depths

Figure 5 illustrates the distribution of earthquake depths. The x‑axis denotes depth (km), and the y‑axis denotes frequency. A predominance of shallow events (near‑surface) is evident from the tall bar near zero, while the overlaid curve indicates that deeper earthquakes are relatively rare.

tu_pian_40.jpg

Figure 5. Distribution of earthquake depths

6.3.3 Earthquake locations by magnitude with global outline

Figure 6 presents a global distribution of earthquake epicenters, indicated by dots. Dot color corresponds to earthquake magnitude, as defined in the adjacent color bar (yellow = high magnitude; dark purple = low magnitude). Longitude and latitude are shown on the x‑ and y‑axes, respectively. Earthquake activity is concentrated along tectonic plate boundaries, with notable clustering in the pacific “ring of fire” region.

tu_pian_41.jpg

Figure 6. Earthquake locations by magnitude with global outline

6.4 Global heatmap of earthquake locations

Figure 7 presents a global heatmap showcasing earthquake frequency, where dark red indicates high-density zones and light blue marks low-density areas. The map plots longitude and latitude along the x- and y-axes, respectively. Notably, seismic activity aligns closely with tectonic plate boundaries, with prominent clusters observed in the Pacific Ring of Fire, as well as regions near South America, Japan, and Indonesia.

tu_pian_42.jpg

Figure 7. Global heatmap of earthquake locations

6.4.1 Magnitude type vs. depth

Figure 8 illustrates the relationship between earthquake magnitude and depth in a scatter plot. Depth (km) is shown along the x‑axis and magnitude along the y‑axis. Colors represent different magnitude measurement types, as defined in the legend. Shallow earthquakes exhibit a broad magnitude range, while deep‑focus events are comparatively rare and generally exceed magnitude 4.

tu_pian_43.jpg

Figure 8. Magnitude type vs. depth

6.4.2 Distribution of magnitude types

Figure 9 displays the distribution of different earthquake magnitude measurement types. The x‑axis lists magnitude types (e.g., ml, md, mb), and the y‑axis indicates their count or frequency. The most common types are ml and md, with ml occurring most frequently. Other types such as mb, Mww, Mwr, and mB_LG are comparatively rare.

tu_pian_44.jpg

Figure 9. Distribution of magnitude types

6.4.3 Magnitude vs. magnitude error

Figure 10 displays a scatter plot illustrating the relationship between earthquake magnitude and measurement error. The x-axis represents magnitude, while the y-axis shows the corresponding error. Lower-magnitude earthquakes (near zero) exhibit greater variability in error, whereas higher magnitudes (above 4) show smaller, more consistent errors indicating improved measurement accuracy with increasing magnitude.

tu_pian_45.jpg

Figure 10. Magnitude vs. magnitude error

6.4.4 Earthquake distribution by state (with other places)

Figure 11 shows that California leads in earthquake occurrences (42.8%), followed by Hawaii (9.9%) and Alaska (7.0%). Less activity is noted in other states, with “Other Places” making up 17.4%. The data highlights a strong concentration of seismic events in specific regions, especially California and Hawaii.

tu_pian_46.png

Figure 11. Earthquake distribution by state (with other places)

6.5 Data preprocessing

To ensure the dataset's relevance, integrity, and suitability for machine learning modeling, a comprehensive data pre-processing pipeline was implemented. Initially, data filtering was performed by retaining only rows labeled as ‘earthquake’ in the type column, excluding other seismic events such as quarry blasts, explosions, and sonic booms. Missing data were addressed by identifying columns with NaN values specifically nst, gap, dmin, horizontalError, magError, and magNst - and removing rows with incomplete entries, reducing the dataset from 7,884 to 6,146 records. Duplicate rows were checked using the df.duplicated() function, confirming no redundancies. Redundant columns such as net and location source were dropped to streamline the dataset. Outlier removal was conducted by excluding records with negative earthquake magnitudes, which lack physical meaning in standard seismic interpretation.

For categorical transformation, the mag_type feature was one-hot encoded into binary columns representing each magnitude type (e.g., magType_mb, magType_ml), enabling effective model input. Finally, feature standardization was applied using Z-score normalization via StandardScaler, ensuring consistent feature scaling across training and test sets to enhance model convergence and performance.

6.5.1 Algorithmic breakdown and discussion

Table 6 provides a detailed overview of the algorithmic breakdown along with relevant discussions.

Table 6. Algorithmic breakdown and discussion

6.6 Correlation analysis

Pearson's correlation coefficient measures the strength and direction of the linear relationship between two continuous variables such as the mag feature and other numerical attributes. Its value ranges from -1 (indicating a perfect negative correlation) to +1 (indicating a perfect positive correlation), with 0 representing no linear association. The coefficient is calculated using the following formula:

r = Σ((xᵢ - mean(x)) * (yᵢ - mean(y))) / √(Σ(xᵢ - mean(x))² * Σ(yᵢ - mean(y))²)

where,

xᵢ, yᵢ: individual data points for the two variables respectively.

x, y: Two variables of dataset.

R: Pearson correlation coefficient.

The resulting matrix has the correlation coefficients between all pairs of numerical features.

The correlation matrix (Figure 12) provides insight into which variables are strongly associated with earthquake magnitude and which exhibit minimal or negative correlations, potentially indicating limited or no predictive value with respect to the target variable. The accompanying heatmap visualizes these relationships among the dataset’s numerical features. Strong positive or negative correlations are depicted in dark red or dark blue, respectively, while weaker correlations appear in lighter shades.

tu_pian_47.jpg

Figure 12. Heatmap image (confusion matrix)

Values near 1 indicate strong positive correlations, while those near -1 indicate strong negative correlations. The diagonal represents perfect correlations, as each variable is correlated with itself. This visualization enables rapid assessment of feature relationships, supporting effective feature selection for modeling.

6.7 Standardization and scaling

Feature scaling was performed to ensure uniformity and prevent dominance of features with larger numerical ranges. Table 7 represents the overview of regression methods for earthquake magnitude prediction.

Table 7. Overview of regression methods for earthquake magnitude prediction

6.8 Cross validation

The study used 5-fold cross-validation to assess model generalizability by splitting data into five parts, training on four, and validating on one ensuring all data points were used. Preprocessing steps like standardization were applied within each fold via pipelines to prevent data leakage. This method provided reliable performance estimates, reduced overfitting, and enabled consistent evaluation using metrics such as MSE, RMSE, NRMSE, and R².

6.9 Training and testing

Models were trained on a normalized dataset using 80% of the data, with preprocessing ensuring consistency. A 5-fold cross-validation approach guided performance evaluation and hyperparameter tuning, emphasizing key metrics and ensemble model selection. The remaining 20% was reserved for testing, with prediction accuracy validated through error metrics and Predicted-vs-Actual plots.

The code defines evaluate_model_cv, which takes a model, training data (x_train, y_train), and a cross‑validation parameter (default = 5). It evaluates models using two scorers MSE (lower is better) and r² (higher is better) via cross_val_score, returning their mean values. Models are wrapped in a pipeline for consistent application and potential preprocessing. A loop iterates through models, creates pipelines, calls the function, and stores results (MSE and r²) in a dictionary, enabling systematic comparison of model performance on the same dataset.

6.10 Ensemble learning

Ensemble learning techniques were applied to combine the strengths of multiple models (Table 8), improving predictive accuracy and generalization. The stacking model implemented in the notebook incorporated three base learners and a meta‑regressor to aggregate their predictions effectively.

This line creates a stacking regressor that combines three machine learning models CatBoost Regressor (with verbose output disabled), LightGBM Regressor, and Random Forest Regressor. In the stacking approach, each base model generates its own predictions, which are then combined by a meta‑model to produce the final output.

Table 8. Rationale for selected regression models

The stacking ensemble model (Figure 13), integrating CatBoost, LightGBM, and Random Forest, delivered the most accurate predictions with the lowest MSE. In contrast, the standalone Lasso regression model showed the highest error, reflecting its relatively limited predictive performance.

tu_pian_48.jpg

Figure 13. MSE access models

The stacking ensemble (Figure 14) learning model (CatBoost, LightGBM, and Random Forest) achieved the highest R² value, reflecting superior predictive performance, while the Lasso regression model recorded the lowest R² value, highlighting its weaker predictive accuracy.

tu_pian_49.jpg

Figure 14. MSE access models

A low RMSE (Figure 15) on the left side of the plot reflects greater predictive accuracy and a better data fit for models such as Ensemble Learning and CatBoost, whereas higher RMSE values seen in models like Lasso indicate reduced effectiveness in minimizing errors.

tu_pian_50.jpg

Figure 15. RSME models

The line plot illustrates (Figure 16) the NRMSE across various machine learning models, ranging from approximately 4% for the high‑performing Ensemble Learning model to over 17% for the low‑performing Lasso model. This variation underscores marked differences in predictive accuracy and suggests potential trade‑offs between model complexity and error minimization.

tu_pian_51.jpg

Figure 16. NRMSE models

The heatmap highlights substantial differences (Figure 17) in MSE across the evaluated machine learning models, with Ensemble Learning and CatBoost achieving the lowest MSE values, while Lasso records the highest, indicating pronounced variation in predictive performance.

tu_pian_52.jpg

Figure 17. The heatmap (MSE in various machine learning models)

The heatmap (Figure 18) illustrates the variation in mean R² values across the models, with Ensemble Learning and CatBoost demonstrating the highest model fit and explaining the greatest proportion of variance in the target variable. In contrast, Lasso is an outlier, displaying a near‑zero R² score.

tu_pian_53.jpg

Figure 18. The heatmap (R-squared values in various machine learning models)

The heatmap presents (Figure 19) RMSE values for the different machine learning models, showing that Ensemble Learning achieves the lowest RMSE, whereas Lasso produces the highest, reflecting notable differences in predictive accuracy.

tu_pian_54.jpg

Figure 19. RMSE across various machine learning models

tu_pian_55.jpg

Figure 20. NRMSE models

The heatmap compares prediction accuracy (Figure 20) across machine learning models, showing that ensemble‑based methods (CatBoost + LightGBM + Random Forest Regressor) with an NRMSE of 4.6% consistently outperform linear regression models such as Lasso and ElasticNet.