A novel time-frequency multilayer network for multivariate time series analysis

Unveiling complex dynamics of natural systems from a multivariate time series represents a research hotspot in a broad variety of areas. We develop a novel multilayer network analysis framework, i.e. multivariate time-frequency multilayer network (MTFM network), to peer into the complex system dynamics. Through mapping the system features into different frequency-based layers and inferring interactions (edges) among different channels (nodes), the MTFM network allows efficiently integrating time, frequency and spatial information hidden in a multivariate time series. We employ two dynamic systems to illustrate the effectiveness of the MTFM network. We first apply the MTFM network to analyze the 48-channel measurements from industrial oil–water flows and reveal the complex dynamics ruling the transition of different flow patterns. The MTFM network is then utilized to analyze 30-channel fatigue driving electroencephalogram signals. The results demonstrate that MTFM network enables to quantitatively characterize brain behavior associated with fatigue driving. Our MTFM network enriches the multivariate time series analysis theory and helps to better understand the complicated dynamical behaviors underlying complex systems.


Introduction
Exploring time-dependent complex system from observed multivariate time series represents a basic problem of great significance in broad-ranging fields [1][2][3]. Lots of solutions such as entropy analysis [4], recurrence plot [5], power spectral density [6], chaos analysis [7], etc have been applied to fulfill this challenging task. However, increasingly complex dynamics hidden in real-world systems seriously limit the understanding of their external system behaviors. Traditional time series analysis methods face enormous burden to address these problems and produce reliable results. In this regards, developing an effective multivariate time series fusion theory and then effectively characterizing the complex system dynamics represents a significant challenge of continuing interests.
During the last decade, complex network theory [8,9] has shown strong adaptability and anti-noise ability for characterizing the interaction among components of a variety of real-world systems [10][11][12][13][14]. Particularly, complex network analysis of time series [15][16][17][18][19][20][21], as an important branch, has been well developed and has gained a lot of attention. A review of complex network analysis of time series see [22]. However, in most cases, multivariate time series, as a record or mirror of real-world systems, exhibit obvious multiple characteristics. Multilayer network [23][24][25][26][27][28], possessing large number of nodes in multiple layers with different types of edges, allows providing a more intuitive and accurate characterization of multivariate time series. Plenty of successful applications have been achieved in diverse fields involving social network [29], multiphase flow [30], transportation system [31], biological systems [32], and information spreading [33,34].
The vast majority of natural systems show strong time-dependent multi-frequency characteristics. It is necessary to consider both time and frequency information during complex system analysis. As an improvement of the Fourier transform, wavelet analysis [35][36][37]  Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. mother wavelet. Here, a novel multilayer network analysis framework, i.e. multivariate time-frequency multilayer network (MTFM network), combining wavelet analysis and mutual information (MI), is proposed to explore the complicated dynamics underlying the observed multivariate time series. We first perform continuous wavelet transform (CWT) [38,39] on each channel of multivariate time series to conduct timefrequency representation (TFR). Then the key frequency (associated with the system dynamics) are screened out and divided into several segments. After inferring MI-based network from all the channels in each segment, MTFM network is derived, where each layer corresponds to one segment of the key frequency. MTFM network, integrating time, frequency and spatial information simultaneously, provides an effective mapping of the considered multivariate time series, which enables to characterize the complex but important system dynamics.
In this paper, two typical complex research objects, namely oil-water flows and driving fatigue in brain, are considered to verify the validity of MTFM network framework. Both of them have great need to reveal the complex behaviors from the observed multivariate time series. For the oil-water flows, two flow media are immiscible with each other and usually distribute themselves in different temporal and spatial structures, known as flow patterns. Different flow patterns possess different flow dynamics, which seriously constrains the understanding of the complex flow behaviors. Many previous studies have been done to understand oil-water flows. Hanafizadeh et al [40] studied the oil-water flow patterns in an inclined pipe through flow experiments. Castro et al [41] researched the interfacial waves of stratified flow patterns in viscous oil-water flows. Azizi et al [42] introduced the artificial neural network and predicted the water holdup in vertical and inclined oil-water flows. Lakshmanan et al [43] carried out systematic flow experiments and then constructed magnetic resonance images for the oil-water flows. Perera et al [44] studied the ECT performance to identify flow regimes in oilwater flows. Gao et al [30] designed a double-layer sensor and explored the flow behavior in oil-water flows. Tan et al [45] studied the factors which affect the flow pattern transition in horizontal oil-water flow, including viscosity of oil, pipe diameter, oil type, and pipe material. Zhao et al [46] presented a microwave dual frequency correction algorithm to measure the water content of oil-water flows. Despite the existing results, the complicated flow mechanisms ruling the transition and evolution of flow patterns are still elusive. On the other hand, driving fatigue arises from long-term continuous and repetitive driving operation, which would impair the drivers' capacity of perception, recognition and vehicle control. Nowadays, electroencephalogram (EEG) signals have been successfully used to quantitatively assess the driver mental state. Wang et al [47] employed wavelet entropy and pulse coupled neural network to detect potential danger during fatigue driving from EEG signals. Hajinoroozi et al [48] introduced the convolutional neural network to predict the driver's cognitive states from EEG signals. Huang et al [49] set up an online closed-loop fatigue detection and mitigation system based on the EEG signals. Fu et al [50] proposed a dynamic model based on Hidden Markov Model to estimate driving fatigue from EEG, electromyogram, and respiration signals. Zhao et al [51] investigated reconfiguration changes in functional networks of different EEG bands during fatigue driving. Min et al [52] detected the driving fatigue through multiple entropy fusion and got high accuracy. Nguyen et al [53] combined EEG and NIRS to predict driving drowsiness via studying different brain regions and EEG bands. Wang et al [54] used power spectrum density and sample entropy in the online detection of mental fatigue based on EEG signals. Though lots of methods have been developed, the brain mechanism related with driving fatigue still remains to be investigated. In the following sections, we take the multivariate time series from oil-water flows and fatigue driving experiments as two typical examples to demonstrate the effectiveness of our method.

Methods
Multivariate time series hold rich system dynamical and behavioral features. Our methods allow to construct a multilayer network from multivariate time series. The schematic diagram of our method is shown in figure 1.
, MTFM network is derived via the following steps. First of all, CWT [41] is first performed on each channel via different scale a and continuous translation b. Compared with Fourier transform and short-time Fourier transform, CWT has good time-frequency localization characteristics, which allows providing a more precise TFR. After the selection of mother wavelet t , y ( ) the analysis process can be expressed mathematically as follows: y ( ) is the scaled and translated wavelet.
y ( ) Mother wavelet t y ( ) must satisfy three allowable conditions, including: And in frequency domain, equation (1) could be interpreted as: X w ( ) and * y w ( )present the Fourier transforms of x t ( ) and t , * y ( ) respectively. Wavelet coefficients in WT exactly reflect the energy intensity and distribution associated with scale a and time t. And higher coefficients indicate better matching between signal and wavelets. In particular, scale a enables to adjust the frequency through the following relationship: where F c means the center frequency of mother wavelet t . y ( ) F a is the pseudo-frequency corresponding to the scale a. The change of scale a leads to the squash and stretch of wavelet, which allows abstracting and mapping the abundant characteristics of time series from the perspective of frequency. Real-world systems always run their complex dynamics and behaviors in a fixed frequency range, called key frequency. And this phenomenon presents obvious system specificity. For multivariate time series x t k p , 1, 2, ..., , channel-based CWT allows constructing a series of wavelet TFR and further positioning the key frequency associated with the studied system. Key frequency (i.e. key scales) are then split into L frequency segments. We average the wavelet coefficients in each segment among scales to get a new coefficient sequence, such as the red lines and yellow lines in TFR (see figure 1). So for each channel of the p-channel time series, its wavelet TFR is transformed into L coefficient sequences (corresponding to L frequency segments). For one frequency segment, we set each channel as a node. The edges between the nodes are derived via the normalized MI (NMI) between coefficient sequences. MI captures the shared information between two variables. Based on the entropy theory, MI between two variables (e.g.
where H(X) and H(Y) denote the entropies of variables X and Y, respectively. And H X Y , ( )is the mutual entropy between variables X and Y. Then Taking X as an example, amplitude interval of X is first split into n l 1 3 = ⌊ ⌋ / parts with same width d. And the jth-interval is:   where p C C , j x k y ( )means the joint probability distribution of X and Y. Different frequency segments allow inferring different MI-based networks, respectively. All of the above can enable us to construct MTFM network with L layers and p nodes in each layer.

Characterizing the complex flow dynamics in oil-water flows
We use MTFM network to characterize the complex flow dynamics in oil-water flows. A distributed conductance fluid sensor, named as High-speed Cycle Motivation Conductance sensor (HCMC sensor), is welldesigned for capturing the multichannel spatial flow information. Three typical flow patterns are studied here, i.e. oil-in-water slug flow (D OS/W), oil-in-water bubble flow (D O/W) and very fine dispersed oil-in-water bubble flow (VFD O/W). Tap-water and No. 3 white oil are set as the flow media. The flow process is controlled via two essential flow parameters, namely total flow velocity (V m ) and water-cut (K w ) (i.e. the volume fraction of water). Different combinations of V m and K w allow producing different flow situations, called as flow conditions. We adjust V m and K w to get lots of flow conditions, classified into the above described three flow patterns. And for each flow condition, water and oil are pumped into a small-diameter pipe respectively to conduct oil-water flow experiment. After systematic flow experiments, a series of 48-channel measurements (time series) are obtained via the HCMC sensor, which reflect different spatial flow states and contain abundant information of the flow structures and behaviors. V m is in the range of 0.0368-0.2579 m s −1 . And K w is set as 60%, 70%, 82%, 95%, 96% and 97%, respectively. For each 48-channel time series, MTFM network is inferred. In this part, key frequency of the oil-water flows are exactly split into L=12 segments. Consequently, we totally have 12 layers in each MTFM network G L , 1, 2, ..., . a = a For each layer G α of the MTFM network, we choose a threshold I to retain the edges with higher weight (i.e. higher NMI between nodes). ⟷ is the shortest weighted path between nodes i and j. f is a map from weight to length. Spectral radius R of network G corresponds to the maximal absolute value of the eigenvalues.
where i l denotes the ith eigenvalue of G and V means the node set. D OS/W and D O/W possess exactly different flow structures and behaviors. When K w is not very high, such as 60%, and V m is low, D OS/W flow appears, where oil bubbles collide with each other and gather into large oil slugs. With the continuous increase of V m , the oil slugs become unstable. When V m reaches a critical value, the oil slugs die away and the D O/W flow occurs, in which oil phase exists in the form of dispersed oil bubbles with diverse sizes. We calculate the weighted global efficiency and spectral radius for all the inferred MTFM networks, corresponding to these two flow patterns, respectively. The results are presented in figures 2(a)-(c). As can be seen, weighted global efficiency and spectral radius present relatively larger values for D OS/W flow. Both of them gradually decrease to follow the transition of flow patterns from D OS/W flow to D O/W flow. The flow process becomes more random during such a transition.
In addition, VFD O/W flow emerges when oil bubbles are stroked into very fine oil droplets. This arises from the enhanced turbulence kinetic energy induced by the higher K w (e.g. 95%) and relatively larger V m (e.g.

Detecting the driving fatigue from multichannel EEG via MTFM network
We utilize MTFM network to detect the driving fatigue from EEG recordings. Via recruiting volunteers, we conducted the simulation driving experiments. A driving simulator, a Neuroscan data acquisition device, a computer, and a screen are used to set up experimental system. We totally recruited eight right-handed volunteers (5 males and 3 females) to conduct the experiments. Formal simulated driving experiment was conducted in the afternoon (14:00-17:00). And each subject simulated for about 90 min. EEG signals is recorded Figure 2. Distributions of spectral radius and weighted global efficiency E g w for different flow conditions with water-cut (a) at sampling rate of 1000 Hz. And the electrodes are arranged in accordance to the standard international 10/20 system. In the first 5 min of the simulation experiment, the subjects are awake. And we define the EEG signals in this period as alert state. Then the EEG signals acquired during the last 5 min are named fatigue state. Left mastoid is set as the reference electrode. For each subject, the continuous EEG signals are segmented into 5 s non-overlapped epochs, and then we obtain 60 alert epochs and 60 fatigue epochs. The raw signals are down-   [51,53]. We then overlap all 13 layers into a weighted projection network P i j , via equation (14). For each generated P , i j , a wide range of sparsity (10%S35%) with an interval of 1% is utilized. This way for determining the connections has been widely used in the brain network construction [55,56]. Sparsity (S) is defined as the ratio of the number of existing edges to the maximum possible edge number in the network. The integrals (over the sparsity range) of two network measures, namely weighted global efficiency E g w (see equation (10)) and graph energy E, are then calculated to quantitatively characterize the network topology. The graph energy E of network G is defined as: where i l denotes the ith eigenvalue of G. We display all the calculated results in the figures 3(a)-(h) corresponding to 8 subjects, respectively. MTFM network analysis shows good performance in distinguishing alert state and fatigue state during the driving. P-values of graph energy E (weighted global efficiency E g w ) between alert state and fatigue state are all much smaller than 0.001 (a benchmark for determining the significant difference from the p-value), indicating that our MTFM network method allows characterizing difference between alert state and fatigue state during driving. For each subject, weighted global efficiency presents growth trends, suggesting that a greater synchronization of neural assemblies is achieved as the brain state changes from alert to mental fatigue, which agree with the findings in [51,57]. We combine weighted global efficiency and graph energy to construct a two-dimensional feature vectors. Then we employ SVM and 10 fold cross-validation (10 times) to classify the alert state and fatigue state for different subjects. The average classification accuracy for all subjects reaches 98.51%. All these results indicate that our MTFM network method enables to explore the changes of brain from alert state to fatigue state during the driving process.

Conclusion
Complex systems contain complex dynamics, and complex dynamics spawn complex systems. These two are inseparable. In this paper, aiming at revealing the complex dynamics from the time-dependent complex system, MTFM network framework is developed to serve as a solution. This framework fully exploits the respective advantages of wavelet transform and MI. And it can effectively reveal complex dynamics from the multivariate time series via integrating time, frequency and spatial information. The successful applications of MTFM network in both fluid dynamics and neuroscience have been demonstrated. Particularly, gradually decreased weighted global efficiency and spectral radius can effectively mirror the complex dynamics ruling the evolution of flow patterns and the transitions among them. While the weighted global efficiency with growth trend presents the changes of brain as driving fatigue occurs. These findings render our method particularly efficient for exploring the system dynamical behaviors from the observed multivariate time series. Considering the generality and flexibility, we hope that MTFM network analytical framework could make contributions in cases where great challenges remain faced for grasping complex system dynamics from observed multivariate time series, such as smart power grids, intelligent transport, and global climate.