Detecting Anomalies in Radon and Thoron Time Series Data Using Kernel and Wavelet Density Estimation Methods

1. Introduction

Many studies in the literature have recognized radon as a potential earthquake precursor. Scientists have studied radon extensively in the context of seismic forecasting. Various other techniques can be found in the literature for earthquake forecasting, such as, abnormal animal behavior [1], unusual geochemical activity [2], variations in seismic wave velocity [3,4] and alterations of several types of electromagnetic precursors e.g. see reviews [5,6]. Radon emission and monitoring have remained a popular technique in earthquake forecasting studies [5,6,7,8,9] and is considered among the precursors with the general model of Lithosphere Atmosphere Ionosphere Coupling (LAIC) e.g. see review [10]. Beside various studies reporting health impacts of elevated levels of radon and thoron e.g., reports in [11], many studies have reported that anomalous fluctuations in radon time-series data have been linked with seismic activity; however, these variations may also stem from additional factors, such as soil porosity, radionuclide composition, lithological characteristics, and meteorological influences, including rainfall, humidity, atmospheric pressure, and temperature [13,14]. Filtering out the effects of these perplexing factors, radon anomalies can serve as a reliable indicator for forecasting impending earthquakes.

Radon time series is analyzed for identification of anomalies within the dataset in order to get a better understanding of variations in radon and thoron concentrations and to establish a possible link between anomalies and forecasting of impending earthquakes e.g. see reviews [5,6,10]. It is very important to explore factors that cause changes in RTTS in order to find a link between these changes and seismic activity. Anomalies in RTTS may arise due to geophysical phenomena and environmental factors. Over the years, many advancements have been witnessed in refining existing and proposing new methods for the detection of anomalies in RTTS. Towards this, investigators have tried to account for large data volumes and the dynamic variability in data patterns. For complex dynamical systems, there is always a need to develop more advanced and intelligent-based anomaly detection techniques to address complicated processes that cause these anomalies [14].

Based upon specific applications, scientists have developed several algorithms, with certain advantages and disadvantages, to detect anomalies in time series. Some of the widely used methods to detect anomalies are clustering based approaches, Bayesian online change point detection (BOCD), auto-encoders, Statistical Process Control (SPC), moving average, seasonal decomposition of time series, Delegated Regressor Method [15] and Autoregressive Integrated Moving Average (ARIMA).

In the present study, anomaly in RTTS data was identified using Kernel Density Estimation (KDE) technique by evaluating the probability density function for the normal data points. The technique rests on the foundation that data points deviating significantly from the distribution are taken as anomaly. The non-parametric nature of KDE technique make it a good choice to detect anomaly. An additional advantage of KDE is that it generates density plot that display a clear graphical representation of data points thereby identifying potential outliers in original time series. KDE has ability to handle complex patterns and abrupt changes in time series making it a robust method for anomaly detection [16].

Kernel density estimation detects anomalies by estimating probability density function (PDE). It is also referred as Parzen’s window [17]. Kernel density estimation automatically estimates density of of the given dataset through automatic data processing. Laxhammar et al. [18] have proposed a Gaussian mixture model (GMM) and adaptive Kernel Density Estimator (KDE) for the purpose of anomaly detection. KDE characterizes the normal data accurately through examining the normality modeling. Anomaly detection results using these two methods shows there is no big difference between these two methods. Moreover, the performance of both approaches is deemed to be suboptimal.

Ramanna et al. [19] proposed a zone-free approach for better identification of areas with diffused and spread seismicity. This zone-free approach used the adaptive kernel density estimation (KDE) technique. The main merit behind the adaptive kernel technique outperforms the fixed kernel technique statistically is that the bandwidth of the kernel is spatially variable. This is based on how sparse or clustered the epicenters are. The fixed kernel approach performs poorly with multimodal and long-tail distributions, despite its general success in density estimation scenarios. Because the activity rate probability density surface is multimodal in nature, the adaptive kernel technique is useful in these circumstances and more pertinent to earthquake engineering.

Zhang et al. [20] provided a density-based, unsupervised method for anomaly detection, determining a smooth and efficient anomaly detection method that can be applied to non-linear systems. They used the Gaussian kernel function to attain smoothness in the data. Additionally, adaptive kernel width is used to improve the discriminating power of the model. In high density regions, wide kernel widths are applied to smooth out the difference between normal samples, and in low-density regions, narrow kernel widths are applied to emphasize the abnormality of potentially anomalous samples

In order to mitigate noise from the input signal, Kordestani et al. [21] tested and presented a wavelet-based approach. A finite number of independent observations is used to estimate probability densities using a technique called wavelet-based density estimation (WBDE). It is not necessary to choose a global smoothing scale in advance, and the suggested method can adapt locally to the smoothness of the density depending on the provided discrete particle data. It also maintains the moments of the particle distribution function with a high degree of accuracy and has no restrictions on the system's dimensionality.

García-Treviño and Barria [22] have proposed online wavelet-based density estimation for non-stationary, rapidly changing streaming data (e.g., sensor data, computer network traffic, ATM transactions, web searches). Their methodology was constructed on a recursive formulation of the wavelet–based orthogonal estimator using a sliding window. Computer experimentation demonstrates that their proposed method works well, involving cases for non-stationary applications where the corresponding probability density function is changing with respect to time, and it also works for cases involving running average-based metrics. They have tested their algorithm for simulated and real-world environments and found that their proposed method with good adaptation capabilities that require a fixed amount of memory

In the current study, Kernel Density Estimation (KDE) was used to detect anomalies in radon and thoron time series data. Histogram-based Density Estimation (HBDE) was also employed. Wavelet-based decomposition was applied to the same dataset to discriminate between noise and denoised signals. Following the decomposition, wavelet-based density estimation (WBDE) was performed. The resulting densities were analyzed in relation to identified anomalies. The study spanned over a period of one year, collecting extensive data to support the analysis. This paper also demonstrates the potential of KDE for reliable anomaly detection in environmental radioactivity monitoring.

2. Materials and Methods

This section first presents the instrumental aspect, namely, the underlying geology, the monitoring apparatus and the related earthquake data. Then, the theoretical aspects are described namely the background of the KDE and WBDE methods. As aforementioned both methods are employed for detecting anomalies in radon and thoron time series. The main objective is to investigate the use of KDE and WBDE to identify patterns and anomalies that precede seismic activity. As will be presented both methods are promising and not extensively used in this specific subject. Moreover. as successful anomaly detection via KDE and WBDE will enhance our tools in forecasting earthquakes. The employed dataset spans over a year, comprising radon, thoron, and meteorological data, while nine earthquakes are included the study.

In the following subsections these will be presented in detail.

2.1. Instrumental Aspects

2.1.1. Geology of the Study Area

The data of the current study has been recorded in Muzaffarabad, the capital of the Pakistan-administered state of the Azad Jammu and Kashmir. Muzaffarabad is located at atitude 34.37002° and longitude 73.47082° at the confluence of two rivers, Neelum and Jhelum. The regional tectonic map of the westernmost part of the frontal Himalaya is shown in Figure 1.

Muzaffarabad district covers a total area of 1,642 km², while Muzaffarabad city itself spans approximately 17 km² [23]. Tectonically, Muzaffarabad is situated within the Hazara-Kashmir Syntaxis (HKS), a geologically significant and seismically active structure located in the western Himalayas. The HKS is formed by the folding of Himalayan thrust sheets in the northeastern region of Pakistan [24]. Several active faults are associated with the eastern limb of the HKS, including the Muzaffarabad Fault (MF), Riasi Fault (RF), Shahdara Fault (SF), Shaheed Gala Fault (SGF), Godri Badshah-Kotli Fault (GK), Pirgali Fault (PF), and Samwal Fault (SF).

The Muzaffarabad Fault is also known as Himalayan Frontal Thrust (HFT), Tanda Fault (TF) and Balakot-Bagh Fault (BBF) [25]. The Muzaffarabad Fault creates the massive active rupture zone, which ranged from one to three kilometers. Landslides, gouges, fractures, and cracks are all common in this rupture zone. The morphology along Muzaffarabad Fault is governed by warps or folded scarps, fault-related fold scarps. The local pop-ups, triangular zones, pressure ridges, and stretching movements are evident across this fault.

The Muzaffarabad Fault was recognized to be active prior to the devastating earthquake that struck on October 8, 2005, with the magnitude of Mw=7.6, which occurred at 8:52 am PST (3:50:40 UTC). The epicenter for this earthquake was situated 18 kilometers north-northeast of Muzaffarabad, the capital of Azad Jammu and Kashmir, with a focal depth of 26 km [27].

2.1.2. Data Acquisition

The radon and thoron radioactive gases were measured using the SARAD RTM 1688-2 active radon thoron monitor [28]. By using the above active monitor, soil radon, soil thoron, and the meteorological parameters temperature, pressure, and humidity, have been continuously monitored from March 2017 up to May 2018. To install the RTM 1688-2 for in-situ measurements, a borehole 1 m deep from the soil’s surface was dug. There, the packer probe of RTM 1688-2 was brought down, properly sealed from the surrounding air so that the pumped gas to be only from the nearby soil. Measuring period was 40 min. In this manner, RTM 1688-2 collected 36 samples per day yielding to a total of 15692 data samples were taken during the study period.

2.1.3. Earthquake Related Data

The earthquake data was obtained from the National Seismic Monitoring Centre of Pakistan [29] and refer to the parallel period of measurements Section 2.1.2 plus-minus one month to allow for pre- and post-earthquake data. In oder to select earthquakes from this period, the earthquake preparation radius of Dobrovolsky et al. [30] was calculated:

$$R_D={10}^{0.43M}$$

(1)

The rear should not that equation (1) is still in extensive use in the related literature. The hypothesis in [30] is that, prior to seismic events, a circle of radius $R_D$ exists around the epicenter, within which, major elastic crustal deformation occurs. In this consensus, it is essential that radon and thoron monitoring stations are installed within the radius $R_D$ in order to ensure (according to [30]) that the station’s measurements are influenced mainly by the seismic activity. This means that to account for (mainly) seismic activity, the distance of the radon-thoron measurement station from the earthquake epicenter $R_E$ must be within the earthquake preparation radius, namely the Dobrovolsky’s radius must greater than the epicentral one, i.e., $R_D$≥$R_E$. Under this criterion, nine earthquakes are identified during the earthquakes’ period. These earthquakes are presented in Table 1 together with the Dobrovolky’s and epenentral radii.

2.2. Theoretical Aspects

2.2.1. Kernel Density Estimation

KDE has been introduced in statistical and geostatistical studies, however nowadays is utilised in physics, public health, astronomy, agriculture, social economy [31] and seismic research [32]. It is considered among most popular non-parametric techniques for the estimation of the of probability density function (PDF) of a time series or a variable of a dataset. According to Zhang et al. [20] KDE can be used in anomaly detection in non-linear systems by detecting the density variations in KDE heatmap, which visualize these anomalies A low probability density in the KDE heatmap is associated with an anomaly in the data series [20].

Mathematical expression for the kernel density can be estimated by $Y(n\times m)$ considering it as a matrix $\left\{y_1,y_2,y_{3,\dots \dots \dots \dots \dots \dots \dots \dots ..}{y}_n\right\}$ of identically distributed samples. These samples are taken from a probability kernel density,$b\left(y\right)$, from the m-dimensional Euclidean space given as [18]:

$$b\left(y\right)={\displaystyle \frac{1}{m}}{\sum }_{i=1}^n{a}^{-m}k\left({\displaystyle \frac{y-y_i}{a}}\right)$$

(2)

where $k$ is the density function and $a$ is the width smoothing parameter. The coefficients ${\displaystyle \frac{1}{m}}$ and $a^{-m}$ are used to normalize the density estimator $k$. Probability density functions that are frequently for $k$ are the Gauss, Laplace, Uniform, Tricube and Epanechnikov ones. The following properties are necessary for the generation of a smoothing kernel [33]:

$$\int k\left(y\right)dy=1$$

(3)

$$\int yk(y))dy=0$$

(4)

$$\int y^{2}k\left(y\right)dy\succ 0$$

(5)

KDE has been employed in detection of anomalies in time series. This is because it offers a nuanced understanding of the anomaly detection by encoding, decoding and evaluating post-training, events against the estimated density, providing hence a comprehensive notion of normality [34]. This is achieved by comparing the density of each sample with the one neighborhood through calculation of the average density of all nearby samples [35]. Anomaly detection via KDE can be accessed in geochemistry [36], pattern recognition [37], droughts [38], magnetic fields [39] and other disciplines.

2.2.2. Histogram-Based Density Estimation

Estimating the densities of any time series data can be also carried out using histogram-based density estimation. It can be easily set and its implementation is of low complexity in time. HBDE distribution density, being discontinuous, restrict it to be good choice for the purpose of anomaly detection [35]. However, HBDE can be effectively used for estimating distribution of data when the actual distribution is not known. Estimated density due to HBDE technique is sensitive to the choice of bin width. While considering limitations of histogram-based density estimation technique, other density estimation models are used for the detection of anomalies.

2.2.3. Wavelet-Based Density Estimation

Wavelets are common mathematical tools used for the analysis and computation of non-stationary time series. Wavelets can be locally adjusted so as to follow the smoothness of any function due to their multi-scale nature [40]. Due to this fundamental characteristic, wavelets have been used in a variety of applications. Their most important application is the denoising of non-stationary signals [41]. The key achievement of wavelet thresholding to reduce noise depends upon the fact that the extension of input signals in a wavelet basis is usually sparse, i.e.,several features of the input signal are effectively summarized by a small fraction of wavelet coefficients [40,412].

Wavelet-based density estimators are actually orthogonal series estimators, that is a broad family of non-parametric techniques [42]. The two most widely used orthogonal series are the Fourier and Hermite functions. In comparison to other orthogonal functions for WBDE, the above functions provide greater flexibility in smoothness and convergence [22]. Initially, in. the so called, conventional wavelet density estimation only the preliminary density was approximated. With the evolution of the knowledge in the field, a recursive data modification was achieved which takes into account the data in calculations and also the exponential depreciation of previous data. Through this updating process the coefficients in orthogonal series estimators are represented by the expected value of each data point is projected across the orthogonal series [43].

Wavelet based estimators have not been employed extensively so far and, as a result, the related papers are scarce during the last decade. Jiang et al. [38] employed also wavelet based estimators for droughts.WBDE has been utilised to estimate the uncertainty and bias of hydro-climatological changes [44], to quantify anomalies for Mobile Information Carriers [45] and to study precipitation extremes and large scale climate anomalies [46].

3. Results and Discussion

Time series data of radon, thoron, temperature, pressure, and humidity parameters have been collected continuously for 436 days. According to the in situ experience and the publications of the team [e.g references in 15 and 41], this time period is adequate for the study of radon and thoron disturbances in association with earthquakes. Table 2 shows the values from the KDE and WBDE methods applied to the radon and thorοn time series during period of study. It can be observed that values of the RTS of the earthquakes 6 and 9 tend to be higher than the values for the other earthquakes. WBDE has a tendency to higher value only for earthquake 6. KDE values for TTS tend to be higher for earthquakes 4 and 5, while from WBDE for earthquake 1 and 2. This differentiation can be attributed to a different sensitivity of each method in respect to the time series and earthquake occurrence. To study the density value tendencies further, the average density values are calculated and then the cut-off limit of average+2standard deviation is arbitrary set. The reasoning is that the density values that are above the cut-off limits are unusually high. These cut-off points are given in the last row of Table 2. Only the densities of the radon time series of earthquake 9 from KDE values and of earthquake 7 from the WBDE values are increased. It is obvious that the pure density values of Table 2 provide non-systematic pre-earthquake estimations and therefore the values are only indicative.

It is obvious fro the above viewpoint, that other types of calculations can be advantageous in analysing the radon and thoron data in relation the the studied earthquakes. Towards this the time evolution of the KDE heat map could provide adequate information (Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6). The reader should note here that according to the colour bar of each Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6, where the density is high, hot colors occur, such as reds and yellows. Taking into account the colour bar values. High density corresponds to profile values, i.e. baseline levels, whereas low density localities are associated to anomalous values and cooler colours such as blues or greens. Finally, time intervals with a significant lack of color or very low density are associated with unusual or irregular measurements.

Under the above consensus, radon (Figure 2), thoron (Figure 3) and atmospheric (Figure 4, Figure 5 and Figure 6) time series shows unusual KDE density values at certain time samples. To interpret the time axis, the reader may recall form Section 2.1.2 that each sample correspond to a 40 min measuring interval. In this view, radon time series show anomalous KDE densities on 12 March 2017 i.e. 11 days prior to the occurrence of the M_W=4.3 (21 March 2017, EQ1) earthquake and 13 days prior to the M_W=2.5 (23 March 2017, EQ2) earthquake. The radon also shows anomalous kernel density in density on 15 August 2017, i.e., 12 days before the occurrence of earthquakes of magnitude M_W=4.8 (27 August 2017, EQ3). The anomaly in KDE density was also observed in radon on 17 September 2017, i.e., 7 days before the occurrence of earthquakes of magnitude M_W=4.6 (23 September 2017- EQ4). Another KDE density anomaly was observed in radon time series on 19 January 2018, i.e., 15 days before the occurrence of earthquakes of magnitude M_W=0.8 (3 February 2018.EQ6). The reader should however note here that the above results are valid in the logic of what KDE heat map provides. The corresponding evidence hence are valid under these constrains and they may be subjected to methodology bias if different mathematical approaches are selected.

Thoron also shows anomalous nature in density (Figure 3) on 17 August 15 2017, i.e., 12 days prior to the occurrence of earthquakes of magnitude M_W=4.8 (27 August 27 2017, EQ3). KDE density density was also observed in thoron on 17 September 2017, i.e., 7 days before the occurrence of earthquakes of magnitude M_W=4.6 (23 September 23 2017, EQ4). Another KDE density anomaly was also observed in thoron on 28 February 2018, i.e., 15 days before the occurrence of earthquake of magnitude M_W=0.8 (03 February 2018, EQ6).

From the results presented so far, it is evident that the KDE method manages to identify anomalies in radon and thoron time series in link with the earthquakes of Table 1. This gain of KDE is noteworthy especially if the reader accounts that according to Frehner at el. [34] Among the techniques specifically developed for anomaly detection, there is no clear winner. It is evident from the recent review of Nikolopoulos et al. [5] on earthquake precursors, that there is no single certain rule to link earthquakes and anomalies but several papers report such links for even mild earthquakes. In the sense of this review, the present results fall within the tendencies of both recent and old papers. Moreover, the recent review of Pulinets et al. [10] the Lithosphere Atmosphere Ionosphere Coupling (LAIC) adopts anomalies from all the parts. The reader should note that the first author of reference [10] is the maim founder of the LAIC model which is used extensively nowadays for the explanation of pre-earthquake anomalies. The present results referring on anomalies from data in soil, i.e., from the lithosphere, is both within the LAIC model and the international literature. The review of Conti et al. [6] also verifies this aspect. Hence it can be supported that both radon and thoron anomalies identified through KDE are supported by evidence. Moreover, the fact that KDE has been applied to several types of anomalies (see references of Section 2.2.1), it is also evident that this method is indeed robust, modern and succeeds in delineating anomalies. The method managed to identify five from the nine earthquakes of Table 1.

Since KDE method is successful in finding radon and thoron anomalies as well as time series anomalies in general, it may be used also for the delineation of anomalies of the atmospheric parameters related to radon and thoron data series. According to this sense, Figure 4, Figure 5 and Figure 6 present the variations of the atmospheric parameters with the KDE method. The reader should note here, that traditionally, radon and thoron series are presented in parallel to atmospheric data series. This is because the atmospheric parameters may evoke anomalies in radon and thoron due to the physics of the lithosphere (underlying geology). Classical papers explain that and the reader may find several such references in the references of the reviews [5,6,10]. It should be noted though, that this is a not necessary condition, meaning that it is not obligatory or priory that some kind radon or thoron anomaly is due to atmospheric parameters disturbances. That further means, that even if some indications may co-exist with atmospheric parameters and radon or thoron anomalies, it is not obligatory to associate the disturbances to variations of the atmospheric parameters ad not to seismic events. After these evidence, the reader may observe that the KDE heat map of the atmospheric parameters id not identical to the one for radon and is completely different from the one of thoron time series. In an alternate interpretation, radon time series may has KDE heat maps of all atmospheric series and no single one. On the contrary, thoron time series has no related heat map. According to this consensus, the claims presented regarding the potential links on the earthquakes of Table 1 and radon time series are valid whereas the claims from thoron time series are more strict. The link between the earthquakes of Table 1 and the radon and thoron time series, is reinforced by the findings of the HBDE method (Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11).

Indeed the HBDE estimation of radon time series are very different from the ones of temperature (Figure 9) and relative humidity (Figure 10). The HBDE of air pressure time series (as with the KDE heat maps) have some similarities with the radon time series implying also a bias due to atmospheric pressure in radon time series. Regarding thoron time series there is again no similarities. Therefore tit can be supported that thoron time series manage to delineate better the link to earthquakes of Table 1, but the claims for the links with radon time series are also valid.

Figure 12 through Figure 16 present the wavelet coefficients of radon, thoron, and meteorological time series using the Bior (Biorthogonal) wavelet. The reader should note here that the black or coloured spots indicate wavelet coefficients with higher amplitudes. These spots may represent the important information in the data at the various scales and locations. Black dots may indicate strong coefficients, making them more prominent. A black spot in the wavelet coefficient plot usually denotes areas with substantial amplitudes of the wavelet coefficients. Areas with abrupt changes, discontinuities, or localized features are represented by black dots. They may show whether a feature is present at that specific resolution and location. In this sense and as can be observed from Figure 12, the levels 1-3 show that some noise is present in the radon time series, since some maximum data points for the levels (1-3) are below or above the blue dotted lines indicating the corresponding thresholds. On the other hand according to level 4, noise is not present in radon time series. This is an alternative implication that radon is co-affected by other factors (as already mentioned mainly atmospheric pressure) but this bias is not systematic and can not be separated from the effect of the lithosphere’s disturbances. In other words, the the wavelet coefficients with the Bior wavelet, imply only some noise present in the radon series and therefore the reported claims between potential links of radon time series and the five earthquakes of Table 1 are more evident. This is reinforced by the fact that when smaller coefficients are set to zero for the radon time series the noise and black spots are removed.

Figure 12. Wavelet coefficients of thorium time seies using Bior wavelet.

Figure 12. Wavelet coefficients of thorium time seies using Bior wavelet.

Figure 14. Wavelet coefficients of air temperature using Bior wavelet.

Figure 14. Wavelet coefficients of air temperature using Bior wavelet.

Figure 15. Wavelet coefficients of air pressure using Bior wavelet.

Figure 15. Wavelet coefficients of air pressure using Bior wavelet.

Figure 16. Wavelet coefficients of air humidity using Bior wavelet.

Figure 16. Wavelet coefficients of air humidity using Bior wavelet.

Up to now, via KDE heat maps, the HBDE method and the four levels of the the wavelet coefficients with the Bior wavelet, it became evident that thoron time series is not affected by the studied meteorological factors. Radon time series are surely not strongly affected and only some bias may be present, mainly, by the atmospheric pressure. To check further if such a bias may affect reason time series Perarson’s correlation matrix was produced. Table 3 presents the main linear correlations between atmospheric-meteorological parameters and radon or thoron time series. As is observed some very low positive correlation of radon is observed with pressure and humidity, (coefficients of 0.121 and 0.432, respectively), whilst radon is negligible negatively correlated with temperature (-0.141). On the other hand, the correlation matrix shows greater effects of the meteorological parameters to thoron time series since the corresponding Parsons’s coefficient are higher. In specific, thoron time series is negligibly positively correlated with temperature, (Pearson’s coefficient of 0.579) and a negligible negative correlation with air pressure and air humidity (Paerson’s coefficients of -0.510 and -0.211, respectively). It can be supported hence, that the potential bias of atmospheric parameters in radon time series (where only the atmospheric pressure shows some increased bias) is not justified by the linear correlation between pressure and radon. Therefore this reinforces the claims that the links between radon time series variations and the earthquakes of Table 1 are the most probable. This is also valid for thoron time series, where the corresponding series is not affected by atmospheric parameter from the other methods, whereas the bias identified from Person’s linear correlation is more prominent (still low though). It can be supported hence from all the evidence presented so far, that the link between the variations of radon and thoron time series and the earthquakes of Table 1 are the most probable.

Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21 present the WBDE plots of radon, thoron, and meteorological time series. As can be observed from Figure 17, low-density areas are observed on 12 March 2017, i.e., 11 days before two earthquakes of magnitudes M_W=4.3 (21 March 2017, EQ1) and M_W=2.5 (23 March 2017, EQ2). According to Figure 18, anomalous WBDE are observed on 12 March 2017, i.e., 11 days before two earthquakes of magnitudes M_W=4.3 (21 March 2017, EQ1) and M_W=2.5 (23 March 2017, EQ2). Peaks, valleys and flat areas are observed in the WBDE plots of the meteorological time series, indicating the anomalies in temperature, pressure, and humidity time series data.

4. Conclusions

This study employed the wavelet transform on radon and thoron time seies with Biorthogonal wavelets to separate noise and extract meaningful signals. Kernel density estimation wavelet-based density estimation methods were then applied to detect anomalies in the radon and thoron time series. The outcomes of this study confirm the presence of anomalies in both radon and thoron time series data. Moreover, the links of these anomalies to seismic events that occurred during the study period are also reported. The KDE applied to the radon and thorough time series identified nine significant anomalies. The anomalies in radon time series data were detected on 12 March 2017, 15 August 2017, 17 September 2017 and on 19 January 2018. These anomalies were identified several days before the corresponding seismic events. Anomalies in thoron time series were observed on 12 March 2017, 15 August 2017, 17 September 2017, and on 28 February 2018, all preceding the seismic events by some days. Five out of nine seismic events with anomalies in both radon and thoron time-series data have been accurately were identified by the Kernel Density Estimation technique. WBDE identified four radon anomalies and two thoron anomalies that were associated to the corresponding seismic activities. Future research should focus on refining these techniques and expanding the dataset to enhance the accuracy and reliability of earthquake prediction based on environmental radioactivity measurements. The installation of more radon thoron stations is always a target in the seismic prone areas of Pakistan.