**1. Introduction**

Investigation of natural data is a meticulous, nontrivial task and is of great concern in the research of different processes and phenomena. Natural data processing and analysis are widely applied in different fields such as technical processes, geophysics, biology, economy, chemistry and so on. In the context of this problem, methods of object statistical analysis and anomaly detection are of special interest. Examples might include the problems of detection of outer space disturbances [1], floods, earthquakes [2], mud flows and landslides, tsunamis [3], changes in the geological environment [4] and other natural emergency situations. The problem of anomaly recognition is also very important in medicine. Timely and highly accurate diagnostics make it possible to detect and identify patient physical state disorders that in some cases may cost someone their life. Such methods should be fast and flexible in order to detect rapidly changing states of systems and objects. These changes characterize anomaly occurrences.

Different methods are used for natural data analysis. A wide spectrum of approaches to this problem is determined by the complexity of the data under investigation and includes deterministic [5] and stochastic [6] tools. At the present time, different combinations of such approaches are being actively developed [7]. As methods for natural data analysis developed, it became evident that classical approaches to modeling and analysis (AR, ARMA models [8], exponential smoothing [9], stochastic approximation [6], etc.) do not allow us to describe properly the complicated nonstationary structure of natural time series.

At the present time, hybrid approaches and the methods [7] allowing one to improve data analysis quality are developing rapidly to solve such problems. Wavelet transform is a flexible instrument and is widely used in the problems of data processing and analysis.

**Citation:** Mandrikova, O.; Polozov, Y.; Mandrikova, B. Natural Data Analysis Method Based on Wavelet Filtering and NARX Neural Networks. *Eng. Proc.* **2023**, *33*, 63. https://doi.org/10.3390/ engproc2023033063

Academic Editors: Askhat Diveev, Ivan Zelinka, Arutun Avetisyan and Alexander Ilin

Published: 16 August 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

A large library of wavelet filters and a wide choice of constructions for signal expansion provide the possibility to adapt this instrument for the data of different structures. For example, we developed the empirical wavelet filters [10], which are effective in image processing applications.

Now, neural network methods [11] are of great importance. They are known to be able to approximate complex non-linear dependencies in data and automate processes. However, to achieve high efficiency of the neural network apparatus, it is necessary to have high-quality training data and their representativeness. The article proposes a method based on combining the operations of wavelet filtering and NARX neural networks. NARX networks are an evolution of the ARIMA class of models. They approximate non-linear time series and may have different architectures. This allows you to flexibly adapt to the tasks of data analysis and forecasting. Wavelet filtering is used to reduce noise and simplify the data structure. This improves the efficiency of the neural network. Similar hybrid approaches to the analysis of complex data have already been used in other works [7,12]. One such work is the combination of ARIMA and discrete wavelet decomposition together with the LSTM neural network [7]. Another example is the work [12], which proposes a discrete wavelet decomposition with ARIMA and a neural network. Based on this method [12], the authors carry out a forecast of hydrological hazards (high water levels (floods), floods, rain floods, etc.). In this paper, the approach of combined application of wavelet filtering and NARX networks is used to improve the efficiency of autoregressive time series models. With the use of wavelet filtering, it is possible to reduce noise and simplify the structure of the initial data. The article shows that this leads to an increase in the performance of the neural network. The wavelet filtering process includes a combination of multiscale analysis [10] and threshold functions.

The article presents the wavelet filtering algorithm and the threshold estimation method. Thresholds were estimated based on statistical operations. The study of the time series of the critical frequency of the ionospheric layer F2 (foF2) was carried out. The structure of the recorded ionospheric data depends on many temporal, weather, positional, and other factors [13]. Time series of the ionosphere demonstrate a regular course and various anomalies in the form and duration. Anomalies can be observed under conditions of increased solar and geomagnetic activity, as well as during seismic activity [13].

It is known that there are cases of ionospheric anomaly negative impacts on highfrequency radio communication and navigation signals from GPS satellite. One example is an event in January 2005. Due to space weather anomalies, 26 flights of United Airlines were redirected to non-optimal routes in order to avoid catastrophic risks of communication losses. An event in October 2003 is also known when the American Wide Area Augmentation System remained inoperative for more than one day and did not navigate vertically because of ionospehric inhomogeneities [14]. As it was mentioned above, traditional approaches to the analyzed ionospheric parameters do not satisfy the requirements of anomaly detection accuracy. In this paper, we suggest a method that confirmed its efficiency in the task of detection of anomalies (ionospehric inhomogeneities) generated by magnetic storms. It was empirically proved that signal pre-processing based on wavelet filtering improves the accuracy of NARX neural network model for ionospheric parameter time series. The paper presents a comparison of the efficiency of the proposed method with direct application of NARX neural network.

#### **2. Applied Data**

In this paper, we applied the ionospheric critical foF2 data [13], which have a one-hour sampling rage. The data have been recorded at the Paratunka site (Kamchatskiy kray, Russia, site coordinates: 53.0 N, 158.7 E) since 1969. Due to instrumentation failures and natural and man-made factors, there are gaps in foF2 data. We used a time series with data gaps no longer than one day. When there were gaps in the data of less than 25 h, they were filled in based on the median method.

#### **3. Method Description**

#### *3.1. Application of Wavelet Transform and Threshold Function*

Natural time series contain noise formed under a large variety of natural and manmade factors [15]. To improve the quality of natural data analysis based on neural networks, according to the papers [16], we applied noise suppression operations. They included the application of a multiscale analysis construction [17] and threshold function. The algorithm for noise suppression is the following:

1. Signal *f*0(*t*) wavelet decomposition into the components:

$$f\_0(t) = \sum\_{j=-1}^{-m} g\_j(t) + f\_{-m}(t),$$

where *<sup>f</sup>*−*m*(*t*) = <sup>∑</sup>*<sup>k</sup> <sup>c</sup>*−*m*,*kφ*−*m*,*k*(*t*) is the smoothed component, *<sup>c</sup>*−*m*,*<sup>k</sup>* = *f*0, *<sup>φ</sup>*−*m*,*<sup>k</sup>*, *<sup>φ</sup>*−*m*,*k*(*t*) = <sup>2</sup>−*m*/2*φ*(2−*mt* <sup>−</sup> *<sup>k</sup>*) is the scaling function, *gj*(*t*) = <sup>∑</sup>*<sup>k</sup> dj*,*k*Ψ*j*,*k*(*t*) are detailing components, *dj*,*<sup>k</sup>* <sup>=</sup> *f*0, <sup>Ψ</sup>*j*,*<sup>k</sup>*, <sup>Ψ</sup>*j*,*k*(*t*) = <sup>2</sup>*j*/2Ψ(2*<sup>j</sup> t* − *k*) is the wavelet, *j* is the decomposition level, the decomposition level *j* = 0 is assumed for the initial signal.

2. Application of the threshold fucntion to the component coefficients *gj*(*t*):

$$T(d\_{j,k}) = \begin{cases} 0, \text{ if } |d\_{j,k}| \le T\_j \\ d\_{j,k'} \text{ if } |d\_{j,k}| > T\_j \end{cases} \tag{1}$$

where *Tj* = *<sup>t</sup>*1<sup>−</sup> *<sup>α</sup>* <sup>2</sup> ,*N*−1*σ*ˆ*j*, *<sup>t</sup>α*,*<sup>N</sup>* are *<sup>α</sup>*-quantiles of Student's distribution, *<sup>σ</sup>*ˆ*<sup>j</sup>* is the sample standard deviation, the decomposition levels *j* = −1, −*m*.

3. Signal wavelet recovery:

$$
\tilde{f}\_0(t) = \sum\_{j,k} T(d\_{j,k}) \Psi\_{j,k}(t) + f\_{-m}(t).
$$

Determination of thresholds *Tj* in operation (1) was based on the paper results [17] and it was assumed that absolute values of the coefficients |*dj*,*k*| of the components *gj*(*t*) were close to zero out of the neighborhood containing signal structural characteristics. In this case, the probability (*α* ≈ 0.99) that the coefficients |*dj*,*k*| with respect to the argument *k* are within the interval (*μ* − 3*σ*; *μ* + 3*σ*) is high. Here, *μ* ≈ 0 is the mathematical expectation of the value |*dj*,*k*|, *σ* is the standard deviation (three-sigma rule [18]). Then, for normally distributed |*dj*,*k*|, we can estimate the thresholds *Tj* with defined confidence level *<sup>α</sup>* as *Tj* = *<sup>t</sup>*1<sup>−</sup> *<sup>α</sup>* <sup>2</sup> ,*N*−1*σ*ˆ*j*, where *<sup>t</sup>α*,*<sup>N</sup>* are the *<sup>α</sup>*-quantiles of the Student's distribution [18].

#### *3.2. Application of NARX Neural Networks*

After noise suppression operations, data are analyzed based on NARX neural networks [11]. We applied feedback NARX networks [11], the architectures of which are illustrated in Figure 1.

The value vector applied to the hidden layer neurons consists of the following components:



The neural network input value ˆ *f*0(*t* + 1) has the form:

$$\hat{f}\_0(t+1) = F[\tilde{f}\_0(t), \tilde{f}\_0(t-1), \dots, \tilde{f}\_0(t-l\_x), \hat{f}\_0(t), \hat{f}\_0(t-1), \dots, \hat{f}\_0(t-l\_y)],\tag{2}$$

where *F*(·) is the neural network display function.

**Figure 1.** Architecture of a generalized recurrent network. (**a**) NARX PA architecture; (**b**) NARX SPA architecture.

Figure 1 shows the NARX network with Series-Parallel Architecture (NARX SPA) (Figure 1b) and NARX network with Parallel Architecture (NARX PA) (Figure 1a) [11]. Analytical form for the NARX PA architecture is represented in expression (2). In NARX SPA, previous values *f* - <sup>0</sup>(*i*) are applied to the input instead of the outputs ˆ *f*0(*i*), *i* = *t*, *t* − *ly*. In this case, the network output value ˆ *f*0(*t* + 1) has the form

$$
\hat{f}\_0(t+1) = F[\overline{f}\_0(t), \overline{f}\_0(t-1), \dots, \overline{f}\_0(t-l\_x), \overline{f}\_0(t), \overline{f}\_0(t-1), \dots, \overline{f}\_0(t-l\_y)].
$$

It is known from [11], that NARX SPA are applied when data of a previous value vector are available. NARX PA networks are applied when only modeling results can be used, for example, when making a forecast with a large step ahead.

The network memory parameters were used in the experiments: *lx* = *ly* = 2, *lx* = *ly* = 5. The Bayesian Regularization algorithm with Levenberg–Marquardt optimization was used to train the networks [11]. That allows us to minimize the linear combination of error squares and weights. After that, their appropriate combination is formed in such a way that the network has the highest degree of generalization.

#### *3.3. Scheme of Method Realization*

The scheme of method realization is illustrated in Figure 2. Anomalies are detected based on estimation of neural network errors

$$\varepsilon\_{i} = \sum\_{i=i-s}^{i+s} |\hat{f}\_{0}(i) - \tilde{f}\_{0}(i)|. \tag{3}$$

There is an anomaly in the data if *ε<sup>i</sup>* > 2*σ* + *εmean*, where *σ* is standard deviation of errors *ε<sup>i</sup>* during periods without anomalies, and *εmean* is mean of errors *ε<sup>i</sup>* in periods without anomalies.

**Figure 2.** Scheme of method realization.

#### **4. Results of Method Application for Ionospheric Data**

Regular foF2 data depend on the time of year and the level of solar activity. Therefore, separate training samples were formed (to build a neural network) for different seasons and levels of solar activity. Two periods of solar activity were studied in the work: high and low. The level of solar activity was estimated based on the average monthly values of the solar index f10.7. We used the ionospheric data recorded during a calm geomagnetic state for the period 1969–2015. When constructing each neural network, training sampling contained about 2000 counts. According to the suggested method, neural networks were trained using the data obtained after wavelet filtering. When making wavelet filtering operations, orthonormal Daubeches third-order wavelets were applied [10]. They were determined by minimizing the data approximation errors. To estimate the method, neural networks were also trained via applying the foF2 initial data. We constructed 16 neural networks using the NARX PA and NARX SPA architectures, which had the memory *lx* = *ly* = 2, *lx* = *ly* = 5.

Tables 1 and 2 show the results of the obtained neural network operation. Root-meansquare (RMS) deviations of network errors were determined as

$$S = \sqrt{\frac{1}{N-1} \sum\_{i=1}^{N} (e\_i - \bar{\varepsilon})^2} \sqrt{}$$

where *e*¯ = <sup>1</sup> *<sup>N</sup>* <sup>∑</sup>*<sup>N</sup> <sup>i</sup>*=<sup>1</sup> *ei*, *ei* <sup>=</sup> <sup>ˆ</sup> *f*0(*i*) − *f* - <sup>0</sup>(*i*).

**Table 1.** Standard deviations of NARX neural network errors with delays (2) and (5). *Tj* is the wavelet filtering.



**Table 2.** Results of operation of neural networks (Summer, low solar activity). *Tj* is the wavelet filtering.

An analysis of the results from Table 1 shows that NARX SPA architecture describes more adequately the foF2 data time series. When using a large number of delay lines, the results become more comparable in favor of NARX SP (Table 1). A comparison of the results in Table 2 shows that application of wavelet filtering operation allows us to significantly increase the neural network operation quality that confirms the suggested method efficiency.

To test the adequacy of neural network models, the autocorrelation of network errors was evaluated based on the Ljung-Box test *Q* = *M*(*M* + 2) ∑*<sup>L</sup> s*=1 *p*2 *s <sup>T</sup>*−*s*, where *<sup>M</sup>* is the observation number, *ps* is the autocorrelation of the *s*-th order, and *L* is the checked lag number. If *Q* > *χ*<sup>2</sup> <sup>1</sup>−*α*,*L*, where *<sup>χ</sup>*<sup>2</sup> <sup>1</sup>−*α*,*<sup>L</sup>* are the distribution quantiles chi-square with *L* freedom degree, the presence of autocorrelation of the *L*-th order in the time series is recognized.

An analysis of the results in Table 3 shows that Ljung–Box test values for the networks constructed without wavelet filtering significantly exceed the critical values *χ*<sup>2</sup> <sup>1</sup>−*α*,*L*. The results indicate a significant correlation of errors of these neural networks and, as a consequence, insufficient quality of foF2 data approximation. When the network memory is increased, data approximation quality improves. The networks, trained applying data wavelet filtering operations, show the best results. For a small memory of networks *lx* = *ly* = 2, in accordance with the Ljung–Box test, network errors are uncorrelated. This confirms the adequacy of the obtained neural network model.


**Table 3.** Estimate of error autocorrelation of neural networks (Summer, low solar activity). *Tj* is the wavelet filtering.

Figure 3 shows the results of neural network operation during a magnetic storm, which occurred on 16 July 2017. The red dashed line in Figure 3 marks the magnetic storm beginning. According to Space Weather site data (http://ipg.geospace.ru, accessed on 1 March 2023), the magnetic storm had a mixed nature, arrival of an accelerated flux from a coronal hole and coronal mass ejection. During strong geomagnetic disturbances, the DST index (DST index of geomagnetic activity) reached the minimum of -72 nTl (Figure 3f,l). Analysis of foF2 initial data (Figure 3a,g) and comparison with the moving median (in Figure 3a,g, the median is a green dashed line) show the distortion of data regular time variation on July 16. That is determined by ionospehric disturbance occurrences. Anomalous changes in foF2 data on July 16 confirm the results of the neural networks, as a significant error increase is observed (Figure 3b–e,h–k). However, as the analysis shows, the application of wavelet filtering improves the network performance quality. The errors of the network with small memory *lx* = *ly* = 2 are close to zero except for the anomalous period.

During the anomalous period, network errors significantly increase and the anomaly is clearly detected.

**Figure 3.** Result of data processing foF2. (**a**,**g**) foF2 data (black) and foF2 median (green); (**b**,**c**) NARX SPA errors with memory 2 and 5, respectively; (**d**,**e**) NARX PA errors with memory 2 and 5, respectively; (**h**,**i**) NARX SPA errors with wavelet filtering and memory 2 and 5, respectively; (**j**,**k**) NARX PA errors with wavelet filtering and memory 2 and 5, respectively; (**f**,**l**) DST. Red dashed line is the magnetic storm beginning.

When comparing the results of NARX PA and NARX SPA architectures, we find that for the delay lines *lx* = *ly* = 2, *lx* = *ly* = 5, NARX SPA architecture detects anomalous changes more clearly both when applying initial data to the network input and when applying the data after wavelet filtering to the input. A comparison of the results of NARX SPA neural networks with moving median confirms the effectiveness of the method. The moving median has errors in the period after the storm on July 18, 2017, which are absent in the neural network model.

#### **5. Conclusions**

The application of wavelet filtering together with a NARX neural network is effective for the task of detection of ionospehric anomalies. Direct application of NARX neural networks does not allow one to approximate ionospheric time variation due to the presence of long-term dependencies. Wavelet filtering operation and application of stochastic thresholds significantly decrease the error level even for small memory networks. Wavelet filtering suppresses noise, simplifies the data structure and, as a consequence, makes it possible to obtain a more accurate NARX neural network model. The application of the Ljung–Box test showed no correlation of network errors and confirmed the adequacy of the obtained neural network model.

Comparison of NARX SPA and NARX PA showed the advantage of the NARX SPA architecture, which makes it possible to obtain a more adequate model for describing the ionospheric time series. Using the example of a strong magnetic storm, the possibility of the method for detecting ionospheric anomalies is shown. Comparison of the NARX SPA network with the moving median, widely used to analyze ionospheric data, showed the efficiency of the proposed approach.

With an extension of the statistics and a more detailed study of ionospheric data during anomalous periods, the authors plan to continue research in this direction.

**Author Contributions:** Conceptualization, O.M. and Y.P.; methodology, O.M., Y.P. and B.M.; software, Y.P.; validation, O.M. and Y.P.; formal analysis, O.M. and Y.P.; investigation, O.M. and Y.P.; resources, Y.P.; data curation, Y.P.; writing—original draft preparation, Y.P. and B.M.; writing—review and editing, O.M.; visualization, O.M. and Y.P.; project administration, O.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was carried out as a part of the implementation of the State Task AAAA-A21- 121011290003-0. The work was carried out by means of the Common Use Center "North-Eastern Heliogeophysical Center".

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
