Nonlinear and Non-Stationary Detection for Measured Dynamic Signal from Bridge Structure Based on Adaptive Decomposition and Multiscale Recurrence Analysis

To test the nonlinearity and non-stationarity of measured dynamic signals from a bridge structure with high-level noise and dense modal characteristics, a method that combines the adaptive signal decomposition with the recurrence analysis is proposed to solve the difficulty of testing nonlinearity and non-stationarity of bridge structure signals. A novel white noise assistance and cluster analysis are introduced to the ensemble empirical mode decomposition to alleviate mode-mixing issues and generate single-mode intrinsic mode functions. Combining the hypothesis-testing scheme of nonstationary and nonlinear synchronization and surrogate techniques, a data-driven recurrence quantification analysis method is proposed and a novel recurrence quantification measure pairs are set up. To demonstrate the efficacy of the proposed methodology, complex signals, which are collected from a carefully instrumented model of a cable-stayed bridge, are utilized as the basis for comparing with traditional nonlinear and non-stationary test methods. Results show that the proposed multiscale recurrence method is feasible and effective for applications to a nonlinear and non-stationary test for real complex civil structures.


Introduction
With the fast advancement of microprocessors, miniature sensors, wired and wireless digital transmission networks, and contemporary monitoring devices are able to collect an enormous amount of signals exhibited from bridge structures. However, it is well known that complex bridge structures show high-level nonlinear and non-stationary behaviors due to a variety of causes such as continuous ambient vibrations, vehicle-bridge coupling effect, cumulative crack growth, and degradation of columns, joints, and beams [1]. This brings considerable challenges for civil engineers to inspect the time-varying performance of complex bridge structures. It is of great significance to consider the characteristics of bridge structure signals for detecting the inherent non-stationary and nonlinear nature of signals.
Conventional stationary analysis methodologies and linear statistical approaches tend to have limitations to capture the nonlinear and nonstationary variations in the signals [2]. For example, the fast Fourier transformation cannot provide the temporal localization of the frequency components and assumes that spectral components are stationary. Furthermore, linear statistical methods, e.g., regression models or the analysis of variance, have certain difficulties to capture the nonlinearity 6.
Decompose r k (t) + ε k E k w i (t) , i = 1, · · · , N, until their first EMD mode and define the (k+1)-th mode using the equation below.
8. Go to step 6 for the next k. 9.
Step 6 to 8 are performed until the obtained residue is no longer feasible to be decomposed (the residue does not have at least two extrema). The final residue is expressed using the equation below.
where K is the total number of modes. 10. Let us define that each IMF contains two parts (i.e., the useful characteristic signal and noise signal), and set I MF j = s j (t) + n j (t). s j (t), which stands for the characteristic signal component, and n j (t) stands for the noise component. The cross-correlation function between different IMFs is shown below.
R ij (τ) = R s i s j (τ) + R s i n j (τ) + R s j n i (τ) + R n i n j (τ) (i = j) 11. It is known that corresponding IMFs of different series of noise have no correlation with each other. Then, we can set R s i n j (τ) = R s j n i (τ) = R n i n j (τ) = 0, and Equation (11) can be expressed by the formula below. R ij (τ) = R s i s j (τ) (i = j) (12) Appl. Sci. 2019, 9, 1302 5 of 21 12. Normalized cross correlation coefficient r ij (τ) is introduced to quantify the correlation between two different IMFs, which can be expressed by the equation below.
r ij (τ) = R s i s j (τ) In general, when r ij (τ) ≥ 0.8, the two IMFs are considered to have a strong correlation and the two IMFs are merged [30]. Therefore, the given signal x(t) can be expressed by the formula below.
where M is the total number of modes after cluster analysis, and M ≤ K.
Observe that ε i can be a different value at each decomposition stage. Concerning the amplitude of the added white noise, Huang [31] suggested using small amplitude values for data dominated by high frequency signals, and vice versa. Following the IEEMD method, we fixed the same SNR for all stages. This value might depend on the level of noise contained in the given signal.
The flow chart of the EEMD method and the IEEMD method is shown in Figure 1.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 5 of 20 Observe that can be a different value at each decomposition stage. Concerning the amplitude of the added white noise, Huang [31] suggested using small amplitude values for data dominated by high frequency signals, and vice versa. Following the IEEMD method, we fixed the same SNR for all stages. This value might depend on the level of noise contained in the given signal.
.  Figure 1 has shown that the crucial point for EEMD to lessen the mode-mixing problem is that the number of IMF of different signal decomposed by EMD must be equal. However, it is difficult to ensure the same IMF number when the signal with different series of noise is decomposed by EMD in each period. For the proposed IEEMD method, the mode-mixing problem that the oscillations of very disparate amplitude appear in a mode is improved by taking advantage of the uniformity of the white noise probability density function. The other mode-mixing problem where oscillations or frequencies with the same time scale are assigned to different estimated modes, is improved by the cluster analysis.

Testing and Verification
In this section, a synthetic signal will be used to test and verify the advantage of the IEEMD method when compared with the EEMD method. The noisy signal from the bridge structures is  Figure 1 has shown that the crucial point for EEMD to lessen the mode-mixing problem is that the number of IMF of different signal x i (t) decomposed by EMD must be equal. However, it is difficult to ensure the same IMF number when the signal x i (t) with different series of noise is decomposed by EMD in each period. For the proposed IEEMD method, the mode-mixing problem that the oscillations of very disparate amplitude appear in a mode is improved by taking advantage of the uniformity of the white noise probability density function. The other mode-mixing problem where oscillations or frequencies with the same time scale are assigned to different estimated modes, is improved by the cluster analysis.

Testing and Verification
In this section, a synthetic signal will be used to test and verify the advantage of the IEEMD method when compared with the EEMD method. The noisy signal from the bridge structures is simulated as follows.
where two cosine signals with center frequency 1 Hz and 5 Hz stand for the "true" signals, which contain the structural dynamical information. 10rand models the monitoring noise with a signal-to-noise ration [32] 2.66, and a peak signal to noise ration [32] 15.95. The sampling frequency is 256 Hz and the sampling period is 10 s. The decomposition results are shown in Figure 2.  Figure 2 has clearly shown that the number of IMFs decomposed by EEMD is 11, which is two more than the IEEMD. By using EEMD, the frequency of IMF4 is 1 Hz. However, the IMF5 and IMF6 both have 5 Hz frequency components, which indicates that different IMFs contain the same frequency components, and the mode mixing problem still exists. By contrast, by using IEEMD presented in this paper, it can be seen clearly that IMF5 and IMF6 are the main frequency components. The frequency of IMF5 is 1 Hz, and the frequency of IMF6 is 5 Hz. The IEEMD method effectively alleviates the mode-mixing issues.
The decomposition error is compared by superimposing all the decomposed IMFs and comparing them with the original signal. The reconstruction error of the two methods is shown in Figure 3. From Figure 3, it can be seen clearly that the error of IEEMD is in the order of 10 −13 and that EEMD is in the order of 1. Clearly, the decomposition accuracy of IEEMD is much higher than that of the original EEMD method.  Figure 2 has clearly shown that the number of IMFs decomposed by EEMD is 11, which is two more than the IEEMD. By using EEMD, the frequency of IMF4 is 1 Hz. However, the IMF5 and IMF6 both have 5 Hz frequency components, which indicates that different IMFs contain the same frequency components, and the mode mixing problem still exists. By contrast, by using IEEMD presented in this paper, it can be seen clearly that IMF5 and IMF6 are the main frequency components. The frequency of IMF5 is 1 Hz, and the frequency of IMF6 is 5 Hz. The IEEMD method effectively alleviates the mode-mixing issues.
The decomposition error is compared by superimposing all the decomposed IMFs and comparing them with the original signal. The reconstruction error of the two methods is shown in Figure 3. From Figure 3, it can be seen clearly that the error of IEEMD is in the order of 10 −13 and that EEMD is in the order of 1. Clearly, the decomposition accuracy of IEEMD is much higher than that of the original EEMD method.
presented in this paper, it can be seen clearly that IMF5 and IMF6 are the main frequency components. The frequency of IMF5 is 1 Hz, and the frequency of IMF6 is 5 Hz. The IEEMD method effectively alleviates the mode-mixing issues.
The decomposition error is compared by superimposing all the decomposed IMFs and comparing them with the original signal. The reconstruction error of the two methods is shown in Figure 3. From Figure 3, it can be seen clearly that the error of IEEMD is in the order of 10 −13 and that EEMD is in the order of 1. Clearly, the decomposition accuracy of IEEMD is much higher than that of the original EEMD method.

Brief Description on Recurrence Plot (RP) and RQA
Recurrence plot (RP) measure recurrences of a trajectory x i ∈ R m in its m-dimensional phase space by the recurrence matrix R, which is defined as the equation below [33].
where N is the number of measured points, ε is a recurrence threshold distance, Θ(·) is the Heaviside function, and · is a norm. Then, R i,j (ε) = 1 if the states of the system at time i and that at time j are in an ε-neighborhood, and R i,j (ε) = 0 otherwise. The phase space trajectory x i is reconstructed from a time series {u i } N i=1 by time delay embedding [33].
where m is the embedding dimension and τ is the time delay. The RP is obtained by plotting the recurrence matrix R in its binary entries: a black dot if R i,j = 1 (recurrence) and a white dot if R i,j = 0 (no recurrence). Thus, by the definition of R, the RP has a black main diagonal line and is symmetric with respect to the main diagonal [33].
Since the initial purpose of RP was to visualize trajectories in the phase space, the interpretation of RP is qualitative and subjective. Therefore, it lacks quantitative standards [1]. Hence, a quantification of the obtained structures of RP is necessary for a more objective investigation of the considered system. To this end, the recurrence quantification analysis (RQA) was proposed as a quantitative description of RP by Marwan et al. [34], Zbilut et al. [35], and Webber et al. [36] whose diagonal and vertical lines are the base [37]. Furthermore, the RQA can provide information about different characteristics of the RP and can be used to detect transitions in dynamical systems. The RQA measures, which are adopted in this paper, are as follows: Recurrence rate, RR, is a measure of the density of recurrence points in the RP [37].
where R i,j (ε) is the recurrence matrix defined in Equation (13) and the main diagonal (i = j) is usually excluded. The RR can be used to reveal changes in dynamic systems. Determinism, DET, is the ratio of recurrence points that form diagonal structures of at least length l min to all recurrence points [37]. where P(l) is the histogram of diagonal lines of length l given the threshold ε. Processes with uncorrelated or weakly correlated, stochastic or chaotic behavior cause none or very short diagonals, whereas deterministic processes cause longer diagonals and less single, isolated recurrence points. Therefore, the DET provides a measure for determinism or predictability of the system. Average diagonal line length, L, is the average time that two segments of the trajectory are close to each other [37].
This can be interpreted as the mean prediction time.
Longest diagonal line, L max , is related to the exponential divergence of the phase space trajectory [37]. (21) where N l = ∑ l≥l min P(l) is the total number of diagonal lines. The faster the trajectory segments diverge, the lower the L max will be. Analogously to L max , the maximal length of the vertical lines, V max , in the RP can be regarded as [37] the following equation. (22) where N v is the absolute number of vertical lines.
Entropy, ENTR, refers to the Shannon entropy of the probability p(l) = P(l)/N l to find a diagonal line of exactly length l in the RP [37].
ENTR reflects the complexity of the RP with respect to the diagonal lines. Therefore, for uncorrelated noise, the value of ENTR is rather small, which indicates its low complexity.
Laminarity, LAM, is the ratio between the recurrence points forming the vertical structures and the entire set of recurrence points [37].
Its definition is analogous to the DET in Equation (19), which represents the occurrence of laminar states in the system without describing the length of these laminar phases.
Analogously to ENTR, a new RQA measure, the Entropy of vertical lines (VENTR) is defined in this paper. where is the histogram of vertical lines of length v given the threshold ε.

Data-Driven RQA
The synchronization analysis is usually used to detect the transitions of different complex systems. The basic types of synchronization in coupled complex systems may be broadly divided into three categories [38]: completely synchronization, generalized synchronization, and phase synchronization. In this section, the relationship between recurrences and phase synchronization of complex systems will be adopted. If two systems synchronize in the same way, their recurrences in phase space are not independent from each other. Hence, by comparing the recurrences of each single system with RQA, some statistical indications about their synchronization type and degree will be set up and tested.

Improved RQA Method
RQA is very sensitive in detecting any dynamic change. However, it can be easily influenced by the value of the setup parameters. One crucial parameter of an RP is the Euclidean distance threshold ε. Too small value of ε may result in quantification of much noise and a very low recurrence rate, which leads to the loss of recurrence information of the underlying system. While, on the other hand, too large values of ε may produce false recurrences or saturate the measurement, which leads to many artefacts [1]. Inspired by the DVV method [39] and the DVRQA method [40], an improved RQA method is proposed here to solve the difficulty of adaptive determination of the Euclidean distance threshold. The basic process is as follows.

1.
Set the initial value of recurrence threshold ε 0 and target recurrence rate vector RR tar , the value of ε 0 can be set up based on the equation [39]: where v k,j is a delay matrix consisting of embedding the dimension and delay time.

2.
Set ε i = ε i−1 and construct a recurrence matrix RP i based on Equation (16). Obtain the RR i value based on Equation (18).

3.
If RR i and RR tar,i satisfy the set tolerance, record the target threshold value ε tar,i . If not, the solution process is performed iteratively until the allowable error is satisfied. 4.
Repeat step 2 and step 3 until the threshold vector ε tar corresponding to RR tar is obtained.

5.
According to every value of ε tar , the RQA is performed to obtain the RQA measure vector RQAs.
The flow chart of the previously mentioned Improved-RQA method is shown in Figure 4.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 20 4. Repeat step 2 and step 3 until the threshold vector corresponding to is obtained. 5. According to every value of , the RQA is performed to obtain the RQA measure vector . The flow chart of the previously mentioned Improved-RQA method is shown in Figure 4.
No，Reset ε i

Selection of Other RQA Parameters and Surrogate Techniques
Generally, the mutual information (MI) function [41] and false nearest neighbors (FNN) [41] method are common and robust tools to, respectively, decide the appropriate time delay τ and embedding dimension m. In this paper, the reasonable value for τ is decided based on the first local minimum of the time-delayed mutual information and the selection of the minimum embedding dimension m, which is based on the false nearest neighbor algorithm proposed in Reference [42].
Theiler et al. [43] introduced the concept of "surrogate data," which has been extensively used in the context of statistical non-stationary and non-linearity testing. A surrogate time series is generated as a realization of the null hypothesis under study. Thus, given an original signal, realizations of this data must be generated by modifying only the desired characteristic of the signal that is being tested while the rest of the characteristics remain the same. In this paper, the classic Fourier-based technique [44] and the iAAFT (Iterative Amplitude Adjusted Fourier Transform, iAAFT) technique [44] for the nonstationary and nonlinearity generation of surrogates that will be used later in a statistical analysis based on RQA.

The Statistical RQA Measures
The statistical RQA measures proposed in this case are based on the computation of an ensemble of surrogates, which are representative stationary and linear realizations of the null hypothesis. A statistical measure is computed for the original given time series and the surrogates. If the statistic is significantly different from the values obtained for the surrogate set, the null hypothesis can be rejected. The details of the analysis process is as follows.

Selection of Other RQA Parameters and Surrogate Techniques
Generally, the mutual information (MI) function [41] and false nearest neighbors (FNN) [41] method are common and robust tools to, respectively, decide the appropriate time delay τ and embedding dimension m. In this paper, the reasonable value for τ is decided based on the first local minimum of the time-delayed mutual information and the selection of the minimum embedding dimension m, which is based on the false nearest neighbor algorithm proposed in Reference [42].
Theiler et al. [43] introduced the concept of "surrogate data," which has been extensively used in the context of statistical non-stationary and non-linearity testing. A surrogate time series is generated as a realization of the null hypothesis under study. Thus, given an original signal, realizations of this data must be generated by modifying only the desired characteristic of the signal that is being tested while the rest of the characteristics remain the same. In this paper, the classic Fourier-based technique [44] and the iAAFT (Iterative Amplitude Adjusted Fourier Transform, iAAFT) technique [44] for the nonstationary and nonlinearity generation of surrogates that will be used later in a statistical analysis based on RQA.

The Statistical RQA Measures
The statistical RQA measures proposed in this case are based on the computation of an ensemble of surrogates, which are representative stationary and linear realizations of the null hypothesis. A statistical measure is computed for the original given time series and the surrogates. If the statistic is significantly different from the values obtained for the surrogate set, the null hypothesis can be rejected. The details of the analysis process is as follows.

1.
The stationary surrogate technique and linear surrogate technique are adopted separately to generate N groups of stationary and linear surrogate signals.

2.
Improved-RQA is carried out for each group of stationary and linear signals and the RQA measure matrix corresponding to the target recurrence rate are obtained. Each column of the RQA measure matrix represents a set of index vectors of the surrogate signals.

3.
Calculate the mean and the variance of the RQA measures to obtain the non-stationarity testing and nonlinearity testing confidence interval

4.
Improved-RQA analysis is carried out for the original given signal. Considering the strong noise test environment of the bridge structure, the noise influence factor is introduced to avoid the strong noise submerging the signal and cover up of the real recursive topology of the signal. The mathematical expression of the RQA measures is as follows.
RQAs ori = a noise ·RQAs improved−RQA (26) where a noise is the noise influence factor. a noise = 1 if the given signal is pure. Otherwise, a noise is selected according to the SNR level.

5.
If each element RQAs ori,i of the given signal is within the confidence interval, the null hypothesis will be accepted. The nonstationary element Ξ S,i = 0 and the nonlinear element Ξ L,i = 0. Otherwise, Ξ S,i = 1, Ξ L,i = 1. Let the total size of the RQAs be N (i = 1, 2, · · · , N), the number of the nonstationary elements and nonlinear elements, which are within the confidence interval, be n s ,n l . When the ratios of n s /N,n l /N are greater than the threshold η, the null hypothesis is accepted. Otherwise, the null hypothesis is rejected. Let the binary non-stationarity statistical RQA measure be SNS and nonlinearity statistical RQA measure index be SNL. SNS and SNL can be calculated as follows. 6.
Let the median value of RQAs be RQAm, then the element pair (RQAm,SNS) and (RQAm,SNL) will be set up as the statistical nonstationary and nonlinear test measures.
The flow chart of the data-driven RQA algorithm is shown in Figure 5.
6. Let the median value of RQAs be RQAm, then the element pair (RQAm,SNS) and (RQAm,SNL) will be set up as the statistical nonstationary and nonlinear test measures. The flow chart of the data-driven RQA algorithm is shown in Figure 5.  Figure 5. A data-driven RQA framework for nonstationary and nonlinear signal testing.

Test Description
A 1:20 scale model cable-stayed bridge subjected to dynamic tests will be taken as the case study for the validation of the proposed algorithm in this section. The model bridge, which is shown in The superstructure of the model bridge is constructed of a steel π-shaped girder with box crosssection on both sides of the girder. The width of the girder is 1.1 m and the height of the girder is 0.081 m. The superstructure is divided into six sections, which are made separately and connected by connecting plates which are numbered inside a circle in the figure e.g.①, ②, etc. The superstructure is supported by 60 pairs (or 120 in total) of stay-cables, which are fabricated from wire-ropes of 10 mm diameter and are anchored to the top of the tower. The tower is constructed of two rectangular (0.21 m × 0.14 m) hollow steel columns with wall thickness of 8 mm, which stands 5.95 meters in height above the footing. Three cross-beams connect the columns at each tower. The two upper crossbeams are at spacing of 0.34 m center-to-center at the top, and the lower cross-beam is at spacing of 1.397 m center-to-center at the columns. The end spans of the superstructure are supported by four Figure 5. A data-driven RQA framework for nonstationary and nonlinear signal testing.

Test Description
A 1:20 scale model cable-stayed bridge subjected to dynamic tests will be taken as the case study for the validation of the proposed algorithm in this section. The model bridge, which is shown in Figure 6, may consist of five spans based on the arrangement of the location of towers and piers.   The superstructure of the model bridge is constructed of a steel π-shaped girder with box cross-section on both sides of the girder. The width of the girder is 1.1 m and the height of the girder is 0.081 m. The superstructure is divided into six sections, which are made separately and connected by connecting plates which are numbered inside a circle in the figure, e.g., 1 , 2 , etc. The superstructure is supported by 60 pairs (or 120 in total) of stay-cables, which are fabricated from wire-ropes of 10 mm diameter and are anchored to the top of the tower. The tower is constructed of two rectangular (0.21 m × 0.14 m) hollow steel columns with wall thickness of 8 mm, which stands 5.95 m in height above the footing. Three cross-beams connect the columns at each tower. The two upper cross-beams are at spacing of 0.34 m center-to-center at the top, and the lower cross-beam is at spacing of 1.397 m center-to-center at the columns. The end spans of the superstructure are supported by four piers. The towers and superstructure are constructed of steel plate, with a yield strength of f y = 235 MPa. The piers are constructed of concrete, with a compressive strength of f c = 30 MPa.
Dynamic properties of the model cable-stayed bridge are characterized by uni-directional piezoelectric accelerometers, which have a sensitivity of 2500 mv/g and an operating frequency in the range of 0.05 to 300 Hz. The arrangement of accelerometers is shown in Figure 6, where 13 solid triangles have been used to denote the location of the vertical accelerometers mounted on the girders. Pulse excitation at the location of the middle span is used to characterize the dynamic properties of the bridge. The damped acceleration attenuation responses of the girder are measured by 13 vertical accelerometers. Different damage locations of the girder is modeled by loosening the bottom bolts of the connecting plate, as shown in Figure 7b. The cable damage is modeled by uninstalling the stay cable, as shown in Figure 7c. In this paper, the responses measured from the girder corresponding to the damage of the connecting plate 2 are used to verify the proposed algorithm.  Dynamic properties of the model cable-stayed bridge are characterized by uni-directional piezoelectric accelerometers, which have a sensitivity of 2500 mv/g and an operating frequency in the range of 0.05 to 300 Hz. The arrangement of accelerometers is shown in Figure 6, where 13 solid triangles have been used to denote the location of the vertical accelerometers mounted on the girders. Pulse excitation at the location of the middle span is used to characterize the dynamic properties of the bridge. The damped acceleration attenuation responses of the girder are measured by 13 vertical accelerometers. Different damage locations of the girder is modeled by loosening the bottom bolts of

Results and Discussions
The response measured from the number 6 accelerometer under damage conditions of the connecting plate 2 is used as an example to show the analysis process. The acceleration waveform in the time domain is shown in Figure 8.  Figure 7b. The cable damage is modeled by uninstalling the stay cable, as shown in Figure 7c. In this paper, the responses measured from the girder corresponding to the damage of the connecting plate ② are used to verify the proposed algorithm.

Results and Discussions
The response measured from the number 6 accelerometer under damage conditions of the connecting plate ② is used as an example to show the analysis process. The acceleration waveform in the time domain is shown in Figure 8. The IMFs of the given signal can be obtained by using the IEEMD method proposed in section 2.2, which is shown in Figure 9.  The IMFs of the given signal can be obtained by using the IEEMD method proposed in Section 2.2, which is shown in Figure 9. The IMFs of the given signal can be obtained by using the IEEMD method proposed in section 2.2, which is shown in Figure 9. The stationary and linear surrogate signals of each IMF can be generated by using the method in Section 3.2.2, which are shown in Figure 10 and Figure 11. The stationary and linear surrogate signals of each IMF can be generated by using the method in Section 3.2.2, which are shown in Figures 10 and 11. The IMFs of the given signal can be obtained by using the IEEMD method proposed in section 2.2, which is shown in Figure 9. The stationary and linear surrogate signals of each IMF can be generated by using the method in Section 3.2.2, which are shown in Figure 10 and Figure 11.  The appropriate time delay and embedding dimension of each IMF can be obtained by using the method explained in section 3.2.2. The mutual information, false nearest neighbors, and phase-space figures of IMF1 are shown in Figure 12. From the figure, we can see that the first local minimum of the time-delayed mutual information is 3∆t (∆t is the Sampling interval) and the false nearest neighbors equal zero when the embedding dimension is 7. The appropriate time delay and   IMF1  IMF2  IMF3  IMF4   IMF5  IMF6  IMF7  IMF8   IMF9 IMF10 IMF11 Figure 11. The linear surrogate signals of each IMF.
The appropriate time delay and embedding dimension of each IMF can be obtained by using the method explained in Section 3.2.2. The mutual information, false nearest neighbors, and phase-space figures of IMF1 are shown in Figure 12. From the figure, we can see that the first local minimum of the time-delayed mutual information is 3∆t (∆t is the Sampling interval) and the false nearest neighbors equal zero when the embedding dimension is 7. The appropriate time delay and embedding dimension of other IMFs is similar to the analysis process of the IMF1. The appropriate time delay and embedding dimension parameters of all IMFs are listed in Table 1. The appropriate time delay and embedding dimension of each IMF can be obtained by using the method explained in section 3.2.2. The mutual information, false nearest neighbors, and phase-space figures of IMF1 are shown in Figure 12. From the figure, we can see that the first local minimum of the time-delayed mutual information is 3∆t (∆t is the Sampling interval) and the false nearest neighbors equal zero when the embedding dimension is 7. The appropriate time delay and embedding dimension of other IMFs is similar to the analysis process of the IMF1. The   Parameters IMF1 IMF2 IMF3 IMF4 IMF5 IMF6 IMF7 IMF8 IMF9 IMF10 IMF11     3  3  5  12  16  18  28  35  45  40  23   m  8  9  6  5  3  2  2  2  2  2  2 Eckmann [33] and Marwan [1] both have stated that the length of the diagonal lines is related to the largest positive Lyapunov exponent, which is usually taken as an indication that the system is chaotic and nonlinear. In addition, the length of vertical lines allows for the investigation of intermittency, even for rather short and non-stationary data series [1]. According to this point, in this paper, the diagonal line-based RQA measures (e.g. Lmax, ENTR, DET) are used as indicators of the nonlinear level. The vertical line-based RQA measures (e.g. Vmax, VENTR, LAM) are used as the indicators of the nonstationary level. The statistical six indicators of each IMF can be calculated by using the method described in section 3.2.3. Due to the length of the paper, we only select the statistical RQA measures of the IMF1 and IMF7 to compare the results, which are shown in Figure  13.   Parameters IMF1  IMF2  IMF3  IMF4  IMF5  IMF6  IMF7  IMF8  IMF9 IMF10 IMF11   3  3  5  12  16  18  28  35  45  40  23  m  8  9  6  5  3  2  2  2  2  2  2 Eckmann [33] and Marwan [1] both have stated that the length of the diagonal lines is related to the largest positive Lyapunov exponent, which is usually taken as an indication that the system is chaotic and nonlinear. In addition, the length of vertical lines allows for the investigation of intermittency, even for rather short and non-stationary data series [1]. According to this point, in this paper, the diagonal line-based RQA measures (e.g., Lmax, ENTR, DET) are used as indicators of the nonlinear level. The vertical line-based RQA measures (e.g., Vmax, VENTR, LAM) are used as the indicators of the nonstationary level. The statistical six indicators of each IMF can be calculated by using the method described in Section 3.2.3. Due to the length of the paper, we only select the statistical RQA measures of the IMF1 and IMF7 to compare the results, which are shown in Figure 13.
From Figures 13 and 14, it can be clearly seen that the nonstationary indicators (Vmax, VENTR, LAM) and the nonlinear indicators (Lmax, ENTR, DET) of IMF1 are all in the confidence interval, which indicates that the stationary and linear hypothesis should be accepted. However, the nonstationary indicators and nonlinear indicators of IMF7 all fall outside the confidence interval, which indicates that the stationary and linear hypothesis should be rejected. The conclusion in this case is consistent with the results shown in Figure 9 in which IMF1 mostly contains the stationary and linear white noise, and IMF7 mostly contains the non-stationary and nonlinear structural dynamic information.
To further prove the effectiveness of the above nonstationary indicators, the time-frequency analysis of IMF1 and IMF7 is performed by the continuous wavelet transform algorithm [45], which is shown in Figure 15. The figure has shown that the IMF1 has a relatively uniform frequency distribution with time. Unlike the former, the frequency range 1-2 Hz of IMF7 is the major frequency components and the main frequency in IMF7 is steadily decreasing with time. The time-frequency analysis results of IMF1 and IMF7 indicate that the IMF1 is stationary and IMF7 is nonstationary, which are consistent with the results of the proposed nonstationary test indicators in this paper. Appl. Sci. 2019, 9, x FOR PEER REVIEW 14 of 20  From Figure 13 and Figure 14, it can be clearly seen that the nonstationary indicators (Vmax, VENTR, LAM) and the nonlinear indicators (Lmax, ENTR, DET) of IMF1 are all in the confidence interval, which indicates that the stationary and linear hypothesis should be accepted. However, the nonstationary indicators and nonlinear indicators of IMF7 all fall outside the confidence interval, which indicates that the stationary and linear hypothesis should be rejected. The conclusion in this case is consistent with the results shown in Figure 9 in which IMF1 mostly contains the stationary and linear white noise, and IMF7 mostly contains the non-stationary and nonlinear structural dynamic information.
To further prove the effectiveness of the above nonstationary indicators, the time-frequency analysis of IMF1 and IMF7 is performed by the continuous wavelet transform algorithm [45], which is shown in Figure 15. The figure has shown that the IMF1 has a relatively uniform frequency distribution with time. Unlike the former, the frequency range 1-2 Hz of IMF7 is the major frequency components and the main frequency in IMF7 is steadily decreasing with time. The time-frequency  From Figure 13 and Figure 14, it can be clearly seen that the nonstationary indicators (Vmax, VENTR, LAM) and the nonlinear indicators (Lmax, ENTR, DET) of IMF1 are all in the confidence interval, which indicates that the stationary and linear hypothesis should be accepted. However, the nonstationary indicators and nonlinear indicators of IMF7 all fall outside the confidence interval, which indicates that the stationary and linear hypothesis should be rejected. The conclusion in this case is consistent with the results shown in Figure 9 in which IMF1 mostly contains the stationary and linear white noise, and IMF7 mostly contains the non-stationary and nonlinear structural dynamic information.
To further prove the effectiveness of the above nonstationary indicators, the time-frequency analysis of IMF1 and IMF7 is performed by the continuous wavelet transform algorithm [45], which is shown in Figure 15. The figure has shown that the IMF1 has a relatively uniform frequency distribution with time. Unlike the former, the frequency range 1-2 Hz of IMF7 is the major frequency components and the main frequency in IMF7 is steadily decreasing with time. The time-frequency The bispectrum-based method [46] is used to test for linearity of IMF1 and IMF7, which can further prove the effectiveness of the above nonlinear indicators. The bi-spectrum-based test results are shown in Figure 16. The figure has shown that the estimated (full line) and theoretical (dashed line) values of R, the interquartile range of the estimated bi-coherence values, are very close, which indicates that the test of IMF1 show strong evidence of linearity. Contrary to the former, the estimated and theoretical values of R of IMF7 are clearly different and the estimated value is large, which indicates that the test of IMF7 show strong evidence of nonlinearity. The results are consistent with that of the proposed nonlinear test indicators in this paper. analysis results of IMF1 and IMF7 indicate that the IMF1 is stationary and IMF7 is nonstationary, which are consistent with the results of the proposed nonstationary test indicators in this paper. The bispectrum-based method [46] is used to test for linearity of IMF1 and IMF7, which can further prove the effectiveness of the above nonlinear indicators. The bi-spectrum-based test results are shown in Figure 16. The figure has shown that the estimated (full line) and theoretical (dashed line) values of R, the interquartile range of the estimated bi-coherence values, are very close, which indicates that the test of IMF1 show strong evidence of linearity. Contrary to the former, the estimated and theoretical values of R of IMF7 are clearly different and the estimated value is large, which indicates that the test of IMF7 show strong evidence of nonlinearity. The results are consistent with that of the proposed nonlinear test indicators in this paper. The analysis process has repeated for other IMFs, and the statistical nonstationary and nonlinear RQA measures pairs of all IMFs have been listed in Table 2. In the table, = 1, ⋯ ,6 refers to the median value of the statistical RQA measures. = 1, ⋯ ,6 refers to the binary statistical RQA measures index.   The bispectrum-based method [46] is used to test for linearity of IMF1 and IMF7, which can further prove the effectiveness of the above nonlinear indicators. The bi-spectrum-based test results are shown in Figure 16. The analysis process has repeated for other IMFs, and the statistical nonstationary and nonlinear RQA measures pairs of all IMFs have been listed in Table 2. In the table, = 1, ⋯ ,6 refers to the median value of the statistical RQA measures. = 1, ⋯ ,6 refers to the binary statistical RQA measures index.  The analysis process has repeated for other IMFs, and the statistical nonstationary and nonlinear RQA measures pairs of all IMFs have been listed in Table 2. In the table, Z i (i = 1, · · · , 6) refers to the median value of the statistical RQA measures. D i (i = 1, · · · , 6) refers to the binary statistical RQA measures index. From Table 2, it can be clearly seen that the binary statistical RQA measures of IMF1~IMF2 are zero and that of IMF3~IMF11 are one, which indicates that IMF1~IMF2 are linear and stationary and IMF3~IMF11 are nonlinear and non-stationary. To prove the correctness of the conclusion based on binary statistical RQA measures, the energy [47] of each IMF is calculated, which is shown in Figure 17. From Table 2, it can be clearly seen that the binary statistical RQA measures of IMF1~IMF2 are zero and that of IMF3~IMF11 are one, which indicates that IMF1~IMF2 are linear and stationary and IMF3~IMF11 are nonlinear and non-stationary. To prove the correctness of the conclusion based on binary statistical RQA measures, the energy [47] of each IMF is calculated, which is shown in Figure  17. According to References [24,47], the first group of IMFs whose energies are monotonically decreasing can be viewed as white noises. From Figure 17, it can be clearly seen that the first groups of IMFs, which can be regarded as noises are IMF1~IMF2. In addition, it has become a consensus that white Gaussian noise is stationary and linear, which is identical with the conclusion based on the binary statistical RQA measures index proposed in this paper.
The RQA measures median value and statistical RQA measure index of each IMF are shown in Figure 18. The left Y-axis refers to the median value of each RQA measure and the right Y-axis refers to the statistical RQA index. The Y-axis headings have the same meaning as Table 2.  With the increase of IMF, the level of non-linearity and non-stationarity increases, which indicates that apparent dynamics of the considered data are mainly controlled by its medium-frequency and low-frequency components. The conclusion happens to be consistent with that of Reference [19]. It Figure 17. The IMFs energy of the given signal.
According to References [24,47], the first group of IMFs whose energies are monotonically decreasing can be viewed as white noises. From Figure 17, it can be clearly seen that the first groups of IMFs, which can be regarded as noises are IMF1~IMF2. In addition, it has become a consensus that white Gaussian noise is stationary and linear, which is identical with the conclusion based on the binary statistical RQA measures index proposed in this paper.
The RQA measures median value and statistical RQA measure index of each IMF are shown in Figure 18. The left Y-axis refers to the median value of each RQA measure and the right Y-axis refers to the statistical RQA index. The Y-axis headings have the same meaning as Table 2. From Table 2, it can be clearly seen that the binary statistical RQA measures of IMF1~IMF2 are zero and that of IMF3~IMF11 are one, which indicates that IMF1~IMF2 are linear and stationary and IMF3~IMF11 are nonlinear and non-stationary. To prove the correctness of the conclusion based on binary statistical RQA measures, the energy [47] of each IMF is calculated, which is shown in Figure  17. According to References [24,47], the first group of IMFs whose energies are monotonically decreasing can be viewed as white noises. From Figure 17, it can be clearly seen that the first groups of IMFs, which can be regarded as noises are IMF1~IMF2. In addition, it has become a consensus that white Gaussian noise is stationary and linear, which is identical with the conclusion based on the binary statistical RQA measures index proposed in this paper.
The RQA measures median value and statistical RQA measure index of each IMF are shown in Figure 18. The left Y-axis refers to the median value of each RQA measure and the right Y-axis refers to the statistical RQA index. The Y-axis headings have the same meaning as Table 2.  Figure 18 shows that the vibration signal monitored from the model bridge has nonstationary and nonlinear components. IMF1 and IMF2 are stationary and linear and IMF3~IMF11 are nonstationary and nonlinear. The median value of statistical RQA measures indicate the nonstationary and nonlinear level of IMF3~IMF11. The nonstationary level and nonlinear level of IMF3~IMF4 are weak, and the nonstationary level and nonlinear level of IMF5~IMF11 are strong. With the increase of IMF, the level of non-linearity and non-stationarity increases, which indicates that apparent dynamics of the considered data are mainly controlled by its medium-frequency and low-frequency components. The conclusion happens to be consistent with that of Reference [19]. It  Figure 18 shows that the vibration signal monitored from the model bridge has nonstationary and nonlinear components. IMF1 and IMF2 are stationary and linear and IMF3~IMF11 are nonstationary and nonlinear. The median value of statistical RQA measures indicate the nonstationary and nonlinear level of IMF3~IMF11. The nonstationary level and nonlinear level of IMF3~IMF4 are weak, and the nonstationary level and nonlinear level of IMF5~IMF11 are strong. With the increase of IMF, the level of non-linearity and non-stationarity increases, which indicates that apparent dynamics of the considered data are mainly controlled by its medium-frequency and low-frequency components. The conclusion happens to be consistent with that of Reference [19]. It is worth noting that the previously mentioned law cannot be obtained only according to the single RQA measure.
On the basis of the nonlinear and non-stationary test of the model bridge vibration signal, the time-varying dynamical characteristics of the model bridge are studied further. The identification of the time-varying frequency damping ratio has been conducted by using the stochastic subspace identification method with a sliding window technique [48], and the results are shown in Figure 19.
is worth noting that the previously mentioned law cannot be obtained only according to the single RQA measure.
On the basis of the nonlinear and non-stationary test of the model bridge vibration signal, the time-varying dynamical characteristics of the model bridge are studied further. The identification of the time-varying frequency damping ratio has been conducted by using the stochastic subspace  Figure 19 has shown that the first two-order frequency and damping ratio of the model bridge are gradually rising with a zigzag line, which indicates obvious time-varying characteristics, and the model bridge system is nonstationary and nonlinear. It is remarkable that the time-varying dynamic information will be lost if the non-stationarity and nonlinearity of the signal is not considered.

Conclusions
Nonlinear and non-stationary detection of complex signals is usually a prerequisite for further developing time-varying system identification methods. In this paper, the framework of the combination of signal decomposition methods with statistical RP and RQA is adopted to solve the difficulty of original RP's application in bridge engineering. The following concluding remarks can be drawn from this study.
1. By introducing a new white noise assistance way by considering the uniformity of white noise probability density function, and cluster analysis to original EEMD, an improved EEMD method is proposed to alleviate the mode mixing issues of the EEMD method. The decomposition and reconstruction of a numerical simulation signal has proven the superiority of IEEMD.
2. By combining the surrogate techniques and hypothesis-testing scheme of nonstationary and nonlinear synchronization, a data-driven RQA method is performed on the single-mode signal components and the true indicator pairs, which contain the statistical RQA measures median and a binary RQA index are set up. The results of the conventional time-frequency analysis method, nonlinear test method, and energy analysis are separately used to prove the effectiveness of the nonstationary, nonlinear RQA measures and binary RQA index. Based on the effective implementation of a model cable-stayed bridge, the frequency-evolving RQA of the given signal indicates that the apparent dynamics are mainly controlled by its medium-frequency and lowfrequency components.   Figure 19 has shown that the first two-order frequency and damping ratio of the model bridge are gradually rising with a zigzag line, which indicates obvious time-varying characteristics, and the model bridge system is nonstationary and nonlinear. It is remarkable that the time-varying dynamic information will be lost if the non-stationarity and nonlinearity of the signal is not considered.

Conclusions
Nonlinear and non-stationary detection of complex signals is usually a prerequisite for further developing time-varying system identification methods. In this paper, the framework of the combination of signal decomposition methods with statistical RP and RQA is adopted to solve the difficulty of original RP's application in bridge engineering. The following concluding remarks can be drawn from this study.

1.
By introducing a new white noise assistance way by considering the uniformity of white noise probability density function, and cluster analysis to original EEMD, an improved EEMD method is proposed to alleviate the mode mixing issues of the EEMD method. The decomposition and reconstruction of a numerical simulation signal has proven the superiority of IEEMD.

2.
By combining the surrogate techniques and hypothesis-testing scheme of nonstationary and nonlinear synchronization, a data-driven RQA method is performed on the single-mode signal components and the true indicator pairs, which contain the statistical RQA measures median and a binary RQA index are set up. The results of the conventional time-frequency analysis method, nonlinear test method, and energy analysis are separately used to prove the effectiveness of the nonstationary, nonlinear RQA measures and binary RQA index. Based on the effective implementation of a model cable-stayed bridge, the frequency-evolving RQA of the given signal indicates that the apparent dynamics are mainly controlled by its medium-frequency and low-frequency components.
Conclusively, although the novel method that combines signal decomposition with statistical RP and RQA was applied to detect the nonlinearity and non-stationarity of signals from a model cable-stayed bridge, this method has the potential to be applied to bridges in real complex operating environment. The difference is that the level of noise contained in those signals is generally high and some effective signal processing technologies such as adaptive noise reduction methods may be needed to improve the SNR.

Future Works
In future works, in addition to continuing with the topic, we will: (i) study the method for detecting and repairing the abnormal data in signals, (ii) study the fast RP and RQA algorithm, (iii) try to conduct nonlinear and non-stationary detection on very long time series from the full-scale bridge, and (iv) apply the method proposed in this paper to other structures such as high-rise buildings as well as large and long tunnels.

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