APMEG: Quadratic Time-Frequency Distribution Analysis of Energy Concentration Features for Unveiling Reliable Diagnostic Precursors in Global Major Earthquakes Short-Term Prediction

1. Introduction

In the recent years, the world has seen devastating earthquakes in various regions, namely the Turkey-Syria earthquake in 2023, Noto earthquake on New Year’s Day in 2024, and most recently, the Myanmar-Thailand earthquake.

Earthquakes are crucial for releasing energy accumulated in the Earth’s brittle outer layer, the lithosphere, through the movement of adjacent sections of the Earth’s crust driven by tectonic stresses. [1]. This energy is measured in Mw or Moment Magnitude Scale, which is the most popular metric used by seismologists worldwide to measure an earthquake in terms of energy released. Mw is derived from the seismic moment of an earthquake, calculated as the product of the fault displacement and the force required to cause the movement. This scale is logarithmic, where a one-unit increase in magnitude represents a tenfold increase in amplitude on a seismogram and approximately 31.6 times more energy release [2].

An early and accurate prediction of earthquake occurrences would lead to better disaster management and preparations for structural damage control, which reduces injuries and deaths caused by earthquakes, as well as economic losses [3,4,5]. As such, earthquake prediction is thought to be the ‘holy grail’ of seismology, as it is near impossible given the various factors and margin of error in time, place, and magnitude [6,7]. The Geological Society of United Kingdom in Earthquake Briefing, a policy document, mentions that currently, deterministically predicting the specific timing and location of earthquakes is not achievable; however, the identification of a ‘diagnostic precursor’ would be essential to enable such predictions [8]. Vladimir Keilis-Borok, the founding director of International Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, in his seminal paper entitled Earthquake Prediction, stated that “Earthquake prediction is necessary to undertake disaster preparedness measures, reducing the damage from earthquakes. This requires that the accuracy of prediction be known, but, contrary to common belief, a timely prediction of low accuracy may be useful” [9]. Perhaps due to these glimpses of hope and optimistic sentiments that there are now many detection and early warning systems for earthquakes using big data through a network of seismic monitoring stations around the world [10,11], including a network of smartphones by the participating population [12]. China, for example, launched currently the largest Earthquake Early Warning system (EEW) in the world, named the National System for Fast Seismic Intensity Report and Earthquake Early Warning, comprising about 2,000 seismic stations, 3,200 strong-motion stations, and 10,000 micro-electro-mechanical system (MEMs) stations [13]. Ironically, the earliest known seismometer used for earthquake detection and by extension, an EEW, to mankind was also invented in China circa 132 CE during the Han dynasty by the legendary astronomer Zhang Heng [14]. While EEWs are crucial in enhancing survivability, and reducing injuries by alerting the masses, its effectiveness is hindered by its alert time-window, which is only several seconds before the earthquake shock reaches the masses. Short-term prediction systems, meanwhile, aim to predict major earthquake events in a location with a longer time-window i.e., hours prior. However, reliable short-term prediction systems for earthquakes remain a challenge, receiving critical responses [15,16] due to the impact of such a prediction system existing. Based on the state-of-the-art literature reviews, the key to predicting seismic events lie in big data, signal processing, and machine learning, including neural networks [17,18,19,20]. Among approaches performed previously are analysis on changes in horizontal positioning of geodetic stations [21,22], as well as analysis on FM radio wave-based signals [18].

Disturbingly, most of the earth’s fault lines that cause earthquakes have yet to experience their largest possible earthquake along their individual fault lines [23]. The analyzed 2023 Turkey-Syria earthquake, 2018 Sulawesi earthquake, 2022 Luzon earthquake, 2016 Kaikōura earthquake, and 2024 Noto earthquake, are the most recent major earthquakes from these earthquake-prone countries that saw earthquakes of magnitudes 7 and greater, in which destruction of physical structures caused by the disaster are almost guaranteed [24]; the 2023 Turkey-Syria earthquake saw the largest earthquake of the region in historical period. As for the selected earthquake incidents, reports estimate the total casualty of over 56,000 deaths, and total structural damages and immediate production capacity losses are estimated upwards of US$101 billion where the latest Noto earthquake losses alone is estimated to cost US$17 billion [25,26,27,28,29,30,31].

As these seismic movements or waves can be measured and recorded, with the intensity of the waves noted as magnitude, the recorded waves are essentially a time-series function describing the occurrences of earthquakes and their respective magnitudes. In a research conducted in 2023 aimed to answer the question “does catalogues of smaller earthquakes contain information about future larger earthquakes?”, John B. Rundle concluded, “that catalogues contain significant information on predictability of future large earthquakes” [32], consistent with his other related research [33,34]. This sentiment is also shared by other researchers [35]. As an earthquake can also be described as the dissipation of accumulated energy from tectonic stress in the lithosphere [1], Time-Frequency Distributions (TFD) analysis, with its energy concentration feature to describe the variation of energy over time, may prove to be useful to understand earthquakes [36]. TFDs are used in signal processing to analyze and characterize signals in the time and frequency domains simultaneously and have been used to analyze earthquake motion measured using accelerograms [37,38,39]. In this research, we propose to analyze earthquake records using TFDs method. While analyzing seismic data using TFD in the form of Short-Time Fourier Transform (STFT) and Wavelet Transform (WT) has been performed previously to obtain geological and stratigraphic information [40] as well as for earthquake pattern recognition and detection [41,42], this research introduces the calculation of mean energy concentration, as well as utilization of Cohen’s class TFD specifically Wigner-Ville (WVD) and Choi-Williams (CWD) distributions, aiming to seek major earthquake precursory patterns in the analyzed results for the purpose of timely short-time earthquake prediction. This method of combining Cohen’s class TFDs and calculation of mean energy concentration has been shown to be able to successfully identify precursory patterns to arrhythmia and ischemia occurrences based on electrocardiogram (ECG) datasets with high accuracy [43,44].

2. Materials and Methods

2.1. Data Sources

Datasets of the Turkey-Syria (2023), Sulawesi (2018), Luzon (2022), Kaikōura (2016), and Noto (2024) major earthquakes were used for the analysis conducted in this research [45,46,47,48,49]. The chosen earthquakes have a similarity in that the major earthquakes reported have Mw magnitudes of 7 or greater, as they are destructive and thus are more important to predict, and it is found previously that such earthquakes have precursors [21,22]. The datasets selected to be analyzed span one month inclusive of the major earthquakes, since this research focuses on short-term prediction in the range of hours and days before the earthquake occurrences. It is noted that the datasets used in this research are sampled using geophonic sensors. For this research, we will only be focusing on the timestamp and magnitude data as the time-series function to undergo Time-Frequency Distributions (TFDs). The earthquakes that are selected to be part of the analysis are the latest major earthquakes of each region or country, and the minor earthquakes spanning 1 month before the respective major earthquakes. It is noted that the dataset obtained for the 5 countries are incomplete with irregular samples, thus a data interpolation method is required to obtain a regularly sampled time-series function for analysis. This is especially so for Noto, Kaikōura, and Luzon earthquake catalogs, as the threshold set for an earthquake to be recorded in their respective catalogues is Mw 3.0.

2.2. Data Interpolation

While the average sampling period can be calculated from observing the summary of each dataset (refer Table 1), all datasets appear to show irregular sampling, which is not meeting a crucial criterion for applying the key function for TFDs - Fast Fourier Transform (FFT) - i.e., the time-series data provided must be, or is assumed to be, regularly sampled [50]. Therefore, a data interpolation method is used to estimate and produce regularly sampled data based on each dataset. The interpolation method used is linear interpolation, which is done by estimating values between two adjacent points using a straight line. As earthquake occurrences are not a high frequency phenomenon, losing high frequency information due to linear data interpolation can be considered minimal impact to analysis result reliability, if not negligible. The linear interpolation formula used to determine the magnitude Mw of the interpolated data at a time T with the known discrete magnitude and time before T (Mw₀, T₀) and after T (Mw₁, T₁) is as follows [51]:

$Mw={\displaystyle \frac{{Mw}_0\left(T_1-T\right)+{Mw}_1\left(T-T_0\right)}{T_1-T_0}}.$ (1)

2.3. Time-Frequency Distributions (TFDs)

Time-frequency analysis is mainly achieved by formulating a function which represents the input signal in both time and frequency spectrum simultaneously, instead of either within time or frequency domain. This function or process can then be changed to suit the various use-cases that contain a time-varying spectrum [52]. There is no objective representation for signal energy changes over time and frequency, and this resulted in several proposals [53]. These proposals can be summarized in the following process in Figure 1.

2.4. Short-Time Fourier Transform (STFT)

The Short-Time Fourier Transform (STFT), also known as spectrogram distribution, is a signal processing method used to analyze multi-component signals. The key idea of this method is to apply FFT over a fixed width window moving over time, resulting in the change of spectrum of a signal over each frame.

2.5. Wigner-Ville Distribution (WVD)

One of the most significant and well-known representations of signal energy distribution is the quadratic class, or Cohen’s class [53]. However, Cohen’s class distribution is hampered by unwanted cross-terms between frequency components [55]. The Wigner-Ville distribution (WVD) belongs to this non-linear TFD category, with its calculation for a continuous signal on a specific frequency $\omega$ is given in (2) [56]. WVD covers most mathematical properties involved in analysis; real-valued, maintains time and shifting frequencies, and satisfies marginal features that generate an instantaneous power and power spectrum.

${WV}_{t,\omega }={\displaystyle \frac{1}{2\pi }}\int x\left(t+{\displaystyle \frac{\tau }{2}}\right)x^{\mathrm{*}}\left(t-{\displaystyle \frac{\tau }{2}}\right)e^{-j\omega \tau }d\tau$ (2)

Where x(t) is the interpolated magnitude signal in time-series with constant sampling rate, and $\mathsf{\omega }$ is a desired single frequency band.

2.6. Choi-Williams Distribution

One of the best methods from the Cohen’s Class of non-linear time-frequency for reducing these terms to the minimum is the Choi-Williams class kernel (CWD), with the distribution on a specific frequency $\omega$ of a discrete signal x(n) expressed as (3) [57].

${CWD}_{n,\omega }=\sum _{k={\displaystyle \frac{N}{2}}}^{{\displaystyle \frac{N}{2}}-1}{e}^{-i2\omega k}\sum _{m=-{\displaystyle \frac{M}{2}}}^{{\displaystyle \frac{M}{2}}}\sqrt{{\displaystyle \frac{\sigma }{4\pi k^2}}}{e}^{-{\displaystyle \frac{{\left(m-n\right)}^2}{4k^2}}}x\left(m+k\right)x^{\mathrm{*}}(m-k)$ (3)

2.7. Continuous Wavelet Transform (CWT)

Continuous Wavelet Transform is more dynamic to represent a signal without constant windowing as compared to STFT [58]. CWT utilizes a group of functions called wavelets, which are generated by shifting and scaling a mother wavelet function. The CWT of a signal x(t) can be defined as follows [58]:

$W_{a,b}={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x\left(t\right){{\psi }_{a,b}}^{\mathrm{*}}\left(t\right)dt$, (4)

where

${\psi }_{a,b}={\displaystyle \frac{1}{\sqrt{a}}}\psi \left({\displaystyle \frac{t-b}{a}}\right).$

2.8. Mean Energy Concentration

The mean energy concentration of TFDs refers to the distribution of energy across the time and/or frequency axes in the TFD representation and provides insights into the dominant time-frequency regions where the signal’s energy is concentrated. As earthquakes are essentially energy released by the lithosphere, mean energy concentration of TFDs would be an appropriate measurement to reflect the nature of earthquakes [59]. For this research, we will be focusing on the mean energy concentration on the time axis of the TFD representation of the earthquake datasets for any precursory behavior prior to the major earthquakes that are not apparent in the raw earthquake catalogues.

2.9. Analysis Procedure

The analysis procedure performed on the earthquake datasets providing the results for this research are shown in Figure 2.

3. Results

Each data is linearly interpolated to preserve the low frequency content. This allows the interpolated data to be analyzed using the time-frequency analysis methods of STFT, WVD, CWD, and CWT, as a constant sampling rate of a signal is a prerequisite for these discrete time-frequency analysis. It is noted that the datasets for these earthquakes have varying ranges of samples, and therefore, several parameters are adjusted, namely the resampling rate for interpolation of data, the window sizes to optimize both time and frequency resolution for STFT, and the time slots used for mean energy concentration calculation. The estimated time of occurrence for each major earthquake is noted and set as the origin time for better visualization and comprehension. Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 visualizes the comparisons for the interpolated data against the STFT, CWD, WVD, and CWT sets for each of the earthquakes respectively.

4. Discussion

From the TFD visualizations of the earthquakes, it can be said that the TFD methods were able to provide spectral information of the earthquakes, as well as the frequency variations over time, providing more insight into the earthquakes. Each of the mean energy concentration results were also able to correctly show the estimated time when the major earthquakes occurred, marking them with a peak.

On observation of the comparisons of the earthquake data, heatmap, and the mean energy concentration derived from the respective TFDs, it is seen that there are a few minor peaks visualized in the mean energy concentration graphs, which were initially barely, if not unnoticeable, in the earthquake signals. The elapsed time in seconds of the peaks since the first sample timestamp of each dataset were noted, and the estimated timestamps in UTC are calculated and tabulated in Table 2, Table 3, Table 4, Table 5 and Table 6.

It is shown that for all major earthquakes in question, there are at least 1 minor earthquake occurring between 11 to 66 hours prior. This timeframe is consistent with results from different approaches [18,21,22]. On the raw earthquake data, it is noted that lower magnitude level earthquakes are significantly more common compared to higher magnitude level earthquakes (Table 7). It is also noted that earthquakes of magnitudes 4.0 to 5.0 are consciously noticeable, and magnitudes 5.0 to 6.0 are noticeable but incidents or injuries caused by these magnitudes are not common [24].

Observing the percentage of minor earthquakes recorded against the percentage of major earthquakes recorded (Table 7), the ratio of recorded minor earthquakes to major earthquakes is staggeringly high. This indicates that the occurrences of minor earthquakes alone do not constitute a precursor to the occurrence of major earthquakes, and as such, observing the time-series data of earthquake signals is insufficient to predict major earthquakes.

In support of this research’s proposed method, the mean energy concentration of the TFD data of the earthquake signals provides possible hints towards the major earthquake occurrences, with peaks corresponding to some of the minor earthquakes. However, it is also noted that not all recorded minor earthquakes have a corresponding mean energy concentration peak, as seen from the overwhelming number of minor earthquakes with magnitudes 1 to 4, compared to only a few mean energy concentration peaks. In other words, the occurrence or lack of minor earthquakes within a timeframe is not the sole factor contributing to the mean energy concentration of that timeframe, but the mean energy concentration is indicative of a signal or data pattern which could indicate a major earthquake.

While the number of major earthquakes analyzed in this research is currently insufficient for rigorous statistical analysis, it can be said that there are no false negative occurrences thus far.

4.1. Endpoint Effect

It is worth noting that the earthquake catalog for Sulawesi earthquake spanned until 2018-11-29 23:54:50 UTC, which is relatively close to the main shock at around 2018-09-28 10:02:00 UTC. This may cause an error during TFD process known as endpoint effect, which causes distorted or inaccurate representations of the signal’s time-frequency content near the edges of the analyzed signal [59]. For this reason, a method to reduce said endpoint effect has been adapted [59].

4.2. TFD Comparisons

Among the visualized results (STFT, CWD, WVD, CWT), STFT visualization is the easiest to understand intuitively followed by CWT, and are more suitable for viewing, while CWD and WVD showed more well-defined mean energy concentration visualizations. In terms of accuracy, CWD shows a more accurate-to-reality graph of mean energy concentration, correctly placing all major peaks at 100% mean energy concentration, as opposed to STFT and WVD (see Figure 5c as opposed to Figure 5a,b), while CWT, although having some similar features in its mean energy concentration graph with STFT, CWD, and WVD such as the presence of minor peaks and major peaks, seems to be the most different from the rest. This can be attributed to the fact that CWT does not produce a time-frequency domain result of the input signal, but rather a time-scale output.

The fundamentals of earthquakes lie in the complex interactions of the Earth’s lithosphere at plate boundaries, and understanding these interactions is the key to studying earthquake mechanisms and potentially predicting their occurrence [60,61]. Due to the non-linearity nature of lithosphere dynamics causing non-stationary, time-varying, and multi-spectral earthquake motion, non-linear TFDs (WVD, CWD) prove to be more effective as compared to linear TFD (STFT), showing more well-defined mean energy concentration visualization. Of the Cohen’s class TFDs used, it can be concluded that CWD is more effective and accurate between the two. This conclusion is expected as CWD does have advantages over WVD in multi-component signal analysis [62]. A summary of the CWD TFDs results across all earthquake datasets are shown in Table 8.

4.3. Limitations and Further Research

This research employs only the magnitude and timestamp of the earthquake datasets in the investigated regions to form a time-series signal in the method of analysis, and therefore, can predict only an estimated time, magnitude, and the region of major earthquake (time, size, and location) [63]. Various other measured parameters have yet to be included in the calculation.

The time window between the precursor pattern to the major earthquake found in this research ranges between 11 to 66 hours, which perhaps is acceptable for a prediction timeframe, providing enough time for authorities to warn the public, however the authors believe that the accuracy of the time window can be further improved, with more complete data at a higher and regular sampling rate. The authors also believe that deep learning-based methods in combination with the proposed TFDs analysis would be able to provide a prediction process with enhanced efficiency and efficacy, possibly solving the issue of sparse data in our proposed process [64].

5. Conclusions

In this study, we have explored the utilization of Time-Frequency Distributions (TFDs) on seismic signals for earthquake prediction, demonstrating the effectiveness of TFD-based techniques and enabling the identification of patterns leading up to major earthquakes. Notably, non-linear quadratic Cohen’s class TFDs, especially Choi-Williams Distribution (CWD), prove to be more effective and accurate in analyzing earthquake records as they are non-stationary, multi-spectral, and time-varying in nature.

Through experimentation on real-world earthquakes, the results indicate that TFD-based techniques show great potential in short-term prediction of seismic events. The usage of TFDs provides valuable insights into the dynamics of earthquakes, while the identification of diagnostic precursor patterns enables improved prediction capabilities. This information can aid in the development of advanced prediction models and early warning systems, and consequently, proactive measures for disaster and potentially reduce human and economic losses caused by earthquakes. With a reliable prediction method, countries expecting a major earthquake in the upcoming years such as Japan [65] and Philippines [66] would be able to better prepare their citizens and crisis teams.

In a similar sentiment that smaller earthquakes contain insight to future larger earthquakes [32,33,34,35], this paper concludes that the analysis of earthquake records using TFDs opens new possibilities in earthquake short-term prediction, a possibility which many considered to be impossible previously. The advancements made in this study contribute to the growing body of knowledge in seismology.