Natural Data Analysis Method Based on Wavelet Filtering and NARX Neural Networks †

: A method for analyzing natural data and detecting anomalies is proposed. The method is based on combining wavelet ﬁltering operations with the NARX neural network. The analysis of natural data and the detection of anomalies are of particular relevance in the problems of geophysical monitoring. An important requirement of these methods is their adaptability, accuracy and efﬁciency. Efﬁciency makes it possible to detect anomalies timely in order to prevent catastrophic natural phenomena. Wavelet ﬁltering operations include the application of a multi-scale analysis construction and threshold functions. The article proposes a wavelet ﬁltering algorithm and a method for estimating thresholds based on a stochastic approach. The operations of the method implementation are described. It is shown that the use of wavelet ﬁltering allows one to suppress noise, simpliﬁes the data structure and, as a result, allows one to obtain a more accurate NARX neural network model. The effectiveness of the method for detecting ionospheric anomalies during periods of magnetic storms is shown using the data of the critical frequency of the ionosphere as an example.


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.
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.

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.

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

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: -Current f 0 (t) and previous f 0 (t − 1), ..., f 0 (t − l x ) values of the input vector; -Output values during the previous time intervalsf 0 (t),f 0 (t − 1), ...,f 0 (t − l y ). The neural network input valuef 0 (t + 1) has the form: where F(·) is the neural network display function.   (Figure 1a) [11]. Analytical form for the NARX PA architecture is represented in expression (2). In NARX SPA, previous values f 0 (i) are applied to the input instead of the outputsf 0 (i), i = t, t − l y . In this case, the network output valuef 0 (t + 1) has the form 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: l x = l y = 2, l x = l y = 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.

Scheme of Method Realization
The scheme of method realization is illustrated in Figure 2. Anomalies are detected based on estimation of neural network errors There is an anomaly in the data if ε i > 2σ + ε mean , where σ is standard deviation of errors ε i during periods without anomalies, and ε mean is mean of errors ε i in periods without anomalies.

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 l x = l y = 2, l x = l y = 5. Tables 1 and 2 show the results of the obtained neural network operation. Root-meansquare (RMS) deviations of network errors were determined as 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 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 χ 2 1−α,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 l x = l y = 2, in accordance with the Ljung-Box test, network errors are uncorrelated. This confirms the adequacy of the obtained neural network model.  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 l x = l y = 2 are close to zero except for the anomalous period.
During the anomalous period, network errors significantly increase and the anomaly is clearly detected. When comparing the results of NARX PA and NARX SPA architectures, we find that for the delay lines l x = l y = 2, l x = l y = 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.

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.

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".