Identifying time-varying frequency and amplitude online in real-life structural vibrations is an essential topic of data processing in structural health monitoring. This paper proposes a novel method for this task. We assume that structural vibration signals are stationary in a short time, thus a spectral analysis method called amplitude and phase estimation is conducted to obtain the amplitude spectrum at corresponding time window, and a post- processing technique is proposed to extract the modal frequency and amplitude from the spectrum automatically. The extracted frequency and amplitude could be regarded as the average value of instantaneous frequency and the average value of instantaneous amplitude during the window, respectively. Due to the instability of measured structural vibrations and the uncertainty of spectral shapes under ambient excitation, Kalman filtering is introduced by taking the signal that reconstructed from the identified frequencies and amplitudes as the prediction to enhance the reliability and quality (signal-to-noise ratio) of the next measured signals. Numerical study is performed to inspect the performance of the proposed method. It is also employed to analyze the vibration signals of actual structures, i.e., a cable of a cable-stayed bridge, a hanger of an arch bridge and the main girder of a suspension bridge. The results show its potential to track frequency and amplitude in structural vibrations Preprint submitted to Engineering Structures November 27, 2021 under environmental measurements. The method is supposed to provide fundamental sup- port for further information obtaining and high-level decision making for structural health monitoring systems. Keywords: online frequency tracking, frequency-squeezing postprocessing, Kalman filtering, structural health monitoring 1. Introduction In recent years, many large-scale civil structures such as bridges, buildings and dams have been equipped with structural health monitoring (SHM) systems, and massive monitoring data have been collected for condition evaluation and damage diagnosis [1, 2]. Frequency as a natural property of structures is one of the most fundamental and important indexes that should be concerned in processing monitoring data. As actual structures in service, especially for large-span bridges bearing varying live loads, exhibit nonstationary and non- linear dynamic characteristics, many researches have been focused on the identification of time-varying frequency. The structural vibration response is generally modeled as the superposition of multiple amplitude modulated (AM) and frequency modulated (FM) components: x(t) =KX k=1xk(t) =KX k=1ak(t) cos ( ϕk(t)), (1) where ak(t) and ϕk(t) represent the instantaneous amplitude (IA) and instantaneous phase (IP), respectively. The phase can be written as ϕk(t) =Rt 0ωk(t)dt+θk, and then the instantaneous frequency (IF) is generally defined as the first derivative of the phase function [3]: ωk(t) =d dtϕk(t). (2) As for structure response, in Eqs. (1) and (2), the amplitude is time-varying mainly due to the randomness of ambient excitation and the energy dissipation; and the frequency is time-varying mainly due to the uncertainty (e.g., constraint uncertainty caused by tem- 2 perature, mass uncertainty caused by passing vehicles) and the degeneration (e.g., stiffness decrease caused by section damage, material corrosion, etc.) of structure itself. Uncertainty results in the fluctuations of frequency around a referring value (i.e., natural frequency) while degeneration results in the permanent natural frequency change that cannot recovery automatically. These fluctuations and evolution laws of frequency and amplitude can be used for fatigue analysis, damage detection [4, 5], etc. Also, they can serve as the extracted underlying features for data-driven structural diagnosis techniques [6]. There are several commonly studied ways for IF identification. One class of method is based on mode decomposition. First, the multi-component signal is decomposed into mono-components called intrinsic mode functions (IMFs) in relative literature, which can be realized by empirical mode decomposition (EMD) [7], analytical mode decomposition (AMD) [8], empirical wavelet transform (EWT) [9], variational mode decomposition (VMD) [10], etc. For more discussions and comparisons on those or other mode decomposition methods, we recommend the review works [11, 12, 13, 14] to interested readers. Second, Hilbert transform (HT) is generally employed to obtain the analytical form zk(t) of IMF xk(t) [15]: zk(t) =xk(t) +ivk(t), v k(t) =H(xk(t)) :=1 πZ∞ −∞xk(τ) t−τdτ, (3) where Hdenotes the HT. The IA is obtained as the envelop of the analytical signal, ak(t) = p x2 k(t) +v2 k(t), and the phase is calculated as ϕk(t) = arctan( vk(t)/xk(t)). According to Eq. (2), it’s able to obtain the frequency as ωk(t) =d dt arctanvk(t) xk(t) =xk(t) ˙vk(t)−˙xk(t)vk(t) x2 k(t) +v2 k(t). (4) Using Eq. (4) to calculate IF avoids phase unwrapping instead of using Eq. (2) directly after obtaining the phase. For this kind of technical route, A. Baccigalupi and A. Liccardo [16], P. Ni et al. [17], Y. Xin et al. [18] respectively employed EMD, VMD and AMD to decompose the measured multi-component responses into limited number of IMFs, and the IF was identified by using HT according to Eqs.(2) and (4). Due to the end effect introduced 3 by HT resulting Gibbs phenomenon in IF function, P. F. Pai and R. Zhong et al. proposed the sliding-window tracking method [19] and the conjugate-pair decomposition method [20] for IF identification after the acquisition of IMFs by EMD. In addition, A. Cicone et al. [21] developed a new definition for IF by mapping the IMF xk(t) toθk(t) in polar coordinate to alleviate the end effect and gained better ability to characterize the local property of a signal. Besides, if the IMF is to be a regular harmonic, the energy scheme, Teager-Kaiser algorithm (TKA) [22], is also usually used, which is applicable to online computing but sensitive to noise. The second class of method for IF identification is based on the time-frequency analysis (TFA) and consists of two steps as well. First, by TFA approaches, the original 1D signal is transformed into a 2D time-frequency (TF) plane. Then the embedded IF trajectories ωk(t) (i.e., frequency value at each instant) in the TF plane need to be extracted for the purpose of further quantitative analysis such as natural frequency based cable force calculating [23, 24]. For the first step, the signal processing community provides a wealth of TFA tools for disclosing IF trajectories in the TF plane. Widely used include the linear TFA methods such as short time Fourier transform (STFT) and wavelet transform (WT), and quadratic TFA methods such as Winger-Vile distribution (WVD) and its variants, Cohen class dis- tributions [25, 26]. Due to the principle of Heisenberg uncertainty, the TF representation (TFR) or spectrogram is unable to achieve high resolution in both the time and frequency directions. Many efforts have been conducted to improve the readability of TFR to present more compact and clearer IF trajectories, and this is helpful to the frequency extraction in the next step. Thus, parametric TFA methods such as chirplet transform (CT) [27] and polynomial CT (PCT) [28], which can be regarded as preprocessing techniques, were pro- posed by specifying the parametric form of TF atoms to match the IF’s variation law of a nonstationary signal. Besides, several postprocessing techniques developed recently can greatly concentrate the energy dispersed around the IF traces in the TF plane. Reassign- ment method (RM) is one of the widely adopted, which relocates the coefficients of STFT along both frequency- and time-axis by replacing them with centroids of the local energy 4 distribution [29]. And one drawback of RM is that it doesn’t support signal reconstruction. Synchrosqueezing transform (SST) proposed by Dubechies et al. [30] can be seen as a spe- cial case of RM by squeezing TF coefficients only along frequency-axis and it’s allowable for signal reconstruction. Differ from RM and SST squeezing spectrogram, synchroextracting transform (SET) proposed by G. Yu [31] only extracts the coefficients of STFT most related to IF trajectories and removes the others, furtherly enhancing the concentration of TFR. Let’s denote the TFR results by those TFA methods as TFR method (t, ω),method ∈ {STFT ,WT,SST,SET...}, (5a) i.e., TFR STFT(t, ω) =Z∞ −∞x(τ)g∗(τ−t)e−jωτdτ. (5b) For the second step, the post-task of extracting IF from the TFR is generally resolved based on the fact that coefficients at the locations of IF trajectories are greater than that at other positions around [32, 33], i.e., ˆωk(t) = arg max ωTFR method (t, ω), (6) so it’s also called ridge detection [34, 35]. To keep the identified IF curves continuous and smooth, many ridge detection algorithms such as dynamic path optimization [36, 37], image processing technique [38], Markov chain Monte Carlo algorithm [34], maximum gradient method [39], ant colony optimization [40], etc., have been presented recently. Of course, these two classes of methods mentioned above are not completely independent. On the one hand, TFA methods that support signal reconstruction can also be used for mode decomposition. On the other hand, the decomposed IMFs by mode decomposition methods can also be transferred into TF plane for further IF extraction[41]. These reviewed methods are summarized in Fig. 1. Existing studies provide powerful and plenty of ways for IF identification. However, nowadays, in the context of SHM, we are not only required to identify time-varying parame- 5 Figure 1: Review (incomplete) of existing methods of IF identification. ters for indicating structural conditions, but also realize this task in an online and real-time manner in order to support timelier anomaly detection and decision-making[42]. This has caused the attentions of some researchers and related approaches for online modal identifica- tion and damage detection has been proposed [43, 44, 45, 46, 47, 48]. However, most existing methods as shown in Fig. 1 are realized offline with batch processing. To be adaptive to online data analysis, some modifications to them are necessary. Nevertheless, when handling with operational data, especially for the data generated by large-scale civil structures, there are some problems that need to be faced with for above methods, either realized offline or supposedly online. In online monitoring environment, there exist very weak excitation pe- riods which means a low signal-to-noise-ratio (SNR). Obviously, since the noise is dispersed among the decomposed components except the residual term, the identified IF according to Eqs. (2) and (4) will be extremely affected. A low SNR also leads the IF trajectory in TF plane submerged in the jungle of noise, increasing the difficulty of ridge detection. And because of the randomness of excitation, the modal components participating in vibration and the number of them are varying in different periods. This becomes a further obstacle for some mode decomposition methods and ridge detection algorithms that need to specify the number of signal components in advance. Additionally, existing methods are most veri- fied on numerical and experimental signals, lacking validity proof directly applied to actual long-time continuous signals. 6 This paper presents a novel approach for online frequency and amplitude identification of time-varying structural vibrations. The approach consists of two aspects, one is the au- tomatic natural frequency and amplitude extracting from the spectrum obtained by the method of amplitude and phase estimation (APES) [49] that possesses higher resolution than FT-based methods; another is embedding this extracting algorithm for a single an- alytical step into a Kalman filtering (KF) framework to obtain more robust and reliable continuous tracking results. Both numerical and real-world signals prove the effectiveness of this proposed approach. The rest of this paper is organized as follows. Section 2 defines the problem of concern and Section 3 gives a full description of the proposed solving methodology. In Section 4, synthetic signals are used to illustrate the performance of the proposed method. In Section 5, three long vibration signals that each lasts for 24 hours belong to three actual bridge structures are analyzed. The last Section 6 is the conclusion of the whole work. 2. Problem definition The signal recorded by the sensor of an SHM system is the discretization of realistic structural vibration, suppose the sampling frequency is Fs, then the signal value at instant tis denoted as xt=x(t) :=x(t·1/Fs). This paper focuses on the online identification of frequency and amplitude and the data processing scheme is shown in Fig. 2. Unlike offline data processing: first, only past data are available; Second, assuming the start time is t= 0, once a sufficiently long data series (denote asx0:= [x0, x1,···, xN]Tin Fig. 2(a)), we called computing frame, has been generated, the analysis of this computing frame has to start immediately; Third, the computing frame in the previous step will be updated by the new samples, called updating block, forming a new computing frame and following the analysis of it. Here, the length of computing frame is fixed asN, and denote the length of updating block as S∈[1, N). Thus, as illustrated in Fig. 2(a), the online data processing is in the manner of updating the Nlengthed computing frame with Snewest samples on and on. Each computing frame xk:= [xkS, xkS+1,···, xkS+N−1] is 7 immediately analyzed to obtain the frequency and amplitude. In our paper, this is realized by using APES [49] to estimate the amplitude spectrum first and then using a postprocessing technique to extract the natural frequency and amplitude, which is introduced in Section 3.1. Figure 2: Online data processing scheme. (a) Fixed length computing frame and block updating model. (b) Illustration of analytical time-lag. As shown in Fig. 2(b), the analytical result is outputted at the end of the computing frame, i.e., instants at [ N−1, S+N−1,···, kS+N−1,···]. The outputting frequency isFs/SHz. Since APES assumes the signal is stationary, the amplitude spectrum obtained by APES characterizes the average properties of the signal during the frame length, which means, the identified parameters at the end of the computing frame are actually around the center of it, i.e., instants at [ N/2, S+N/2,···, kS+N/2,···]. Therefore, there exists an analysis delay about half the length of the computing frame. As for TFA methods which perform offline, the IF trajectory will also be delayed for half the window interval unless zeros are padded on both ends [50]. The delay in online analysis can be attributed to the cold start. More precisely, we define the analysis delay as the difference between the time of obtaining the parameter and the time that it actually locates, and take the time cost for data transmission and computation into accounts, then the total analysis delay τanalysis can be expressed as τanalysis =τmethod +τcomputation +τtransmission ≈τmethod , (7) 8 where τmethod represents the lag introduced by the method itself, e.g., in this paper, τmethod≈ N/2 as mentioned above, τtransmission andτcomputation represent the lags caused by data transmission and computation. Among them, τmethod ≫τcomputation ≫τtransmission , the analysis delay τanalysis is dominated by τmethod . In order to enhance the real-time performance of identification, a method that a small computing frame length Nis available should be adopted. This is also one of the points that we employ APES for spectra estimate in this paper. For vibration signals of civil structures, the analysis delay is several seconds of proposed method (longer for lower frequency struc- tures). What’s more, to ensure the continuity of data processing and avoid accumulation, it has to be satisfied that τcomputation +τtransmission < S/F s. (8) So,Sis a trade-off between the output frequency Fs/Sand the threshold of computation costS/F s. 3. Methodology First, for every computing frame, APES is employed to estimate its spectrum and fre- quency squeezing postprocessing (FSP) developed in our previous work [51] is used to con- centrate the spectral lines distribution. To keep the paper as self-contained as possible, we will briefly introduce the methods of APES and FSP in section 3.1. Then, a new method is proposed to automatically extract frequency and amplitude from the spectra obtained by APES and FSP (section 3.2). Last, KF is embedded in the online analysis for better parameters tracking (section 3.3). 3.1. Recaps of APES and FSP 3.1.1. APES method APES was proposed by L. Li and P. Stoica in 1996 [49], which is a nonparametric adaptive filter-bank method with higher resolution than FT based methods [52]. Besides, it is able 9 to adopt shorter computing frame for exact frequency and amplitude estimation [51]. Consider a discrete one-dimensional (1D) time series {x(n)}N−1 n=0, for a frequency ω(in Hz) of interest, it can be decomposed into a sinusoidal term with amplitude α(ω) and a residual term εω(n) that contains all interference caused by sinusoidal components at other frequencies as well as the random noise: x(n) =α(ω)ej2πωn+εω(n), n = 0,1,···, N−1. (9) Arrange the Nconsecutive samples in an M×LHankel matrix: Y= y0y1···yL−1 = x(0) x(1)···x(L−1) x(1) x(2)··· x(L) ............ x(M−1)x(M)···x(N−1) , (10) where M≤N/2 is a user-defined parameter, and L=N+ 1−M. In the matched-filterbank interpretation of APES, the core is to design a data dependent finite-impulse-response (FIR) filter h(ω)∈CM×1that passes the sinusoid with frequency ω undistorted and suppresses sinusoids with any other frequency as much as possible [53]: min h(ω),α(ω)L−1X l=0|hH(ω)yl−α(ω)ej2∗πωl|2s.t.fH M(ω)h(ω) = 1 , (11) where fM(ω) = 1ej2πω···ej2πω(M−1)T , (12) and (·)Hand (·)Tdenote conjugate transpose and transpose of a vector or matrix, respec- tively. The closed solution to Eq. (5) is h(ω) =Q(ω)−1fM(ω) fH M(ω)Q(ω)−1fM(ω), (13) 10 α(ω) =hH(ω)g(ω) =fH M(ω)Q(ω)−1g(ω) fH M(ω)Q(ω)−1fM(ω), (14) where g(ω) =1 LL−1X l=0yle−j2πωl=1 LYf∗ L(ω), (15) Q(ω) =1 LYYH−g(ω)gH(ω) =R−g(ω)gH(ω), (16) and (·)∗denotes conjugate of a vector or matrix. For more details about the derivation of APES, please refer to literature [52, 53]. And we recommend the matrix implementation of it when coding as illustrated in [51] for time- saving, which is also listed in Algorithm 1 here. Algorithm 1 APES. 1:Input : Data series {x(n)}N−1 n=0, calculating frequency ω= [ω1, ω2,···, ωK]T, filter order Mand sampling rate Fs 2:Calculate parameter L=N+ 1−M, matrices FM= [fM(ω1),fM(ω2),···,fM(ωK)] andFL= [fL(ω1),fL(ω1),···,fL(ωK)] by Eq. (12) 3:Construct Hankel matrix Yby Eq. (10) and calculate R=L−1(YYH) 4:Calculate vectors: a= diag( FH MR−1FM),b=L−1diag(FH MR−1YF∗ L) and c= L−2diag(FT LYHR−1YF∗ L) 5:Calculate amplitude α=b÷(a−a×c+b×b∗), here ’ ÷’ and ’ ×’ denote dot products (element-wise) 6:Output : Frequency ωand amplitude α 3.1.2. FSP method FSP as a postprocessing technique for concentrating the spectral lines consists of three steps: (1) amplitude pretreatment, (2) frequency squeezing and (3) amplitude restoration and zero setting. First, pretreat the APES spectrum α= [α1, α2,···, αK]Tby normalization and expo- nentiation: ¯αi= (αi/max(α))m, i= 1,2,···, K, m > 1. (17) Second, squeeze uniformly distributed spectral lines by considering consecutive 2 P+ 1 spectral lines (¯ αi−P,···,¯αi,···,¯αi+P) centered on frequency ωiand replacing the value of 11 ωiwith the centroid of the 2 P+ 1 spectral lines in frequency direction: ωn+1 i=min(i+P,K)X k=max(1 ,i−P)ωn k¯αk/min(i+P,K)X k=max(1 ,i−P)¯αk, i= 1,2,···, K, P ≥1. (18) By iteratively conducting Eq. (18), the original spectral lines will be gradually concentrated around the positions of local highest peaks. Here, the convergence criteria is set as: the frequency difference between adjacent steps should be less than a threshold δs, i.e.,∥ωn+1− ωn∥2≤Kδs. A maximum step is also set to enforce the stop of iterations. Third, restore the amplitude that was changed in step 1 by assigning the original am- plitude vector αto the squeezed frequency vector ˆωin subscript order. Meanwhile, set the amplitudes at both ends of the frequency vector ˆωand that corresponding to frequencies at clusters’ edges as zeros to obtain more clear result presentation. That is, ˆαi=  0, i∈Ω ={1, K} ∪ { (j, j+ 1)|ˆωj+1−ˆωj> δω}, αi, i∈ {1,2,···, K} \Ω,(19) where δωis a threshold used to judge the edges of frequency clusters. The pseudo-code of FSP is also listed in Algorithm 2 here. Algorithm 2 FSP. 1:Input : Frequency ω=ω1= [ω1, ω2,···, ωK]Tand amplitude α= [α1, α2,···, αk]Tby Algorithm 1, order m, step size P, thresholds δsandδω, max iterations I,n= 0,ω0=0 2:Amplitude pretreatment by Eq. (17) 3:Frequency squeezing: 4:while ∥ωn+1−ωn∥2> Kδ sandn < I do 5:n=n+ 1 6: Moving frequency by Eq. (18) 7:end while 8:Amplitude restoration and zero setting by Eq. (19) 9:Output : processed frequency ˆω=ωn+1and amplitude ˆα 3.2. Automatic frequency and amplitude extraction from the spectrum For the purpose of quantitatively analyzing the changing laws of frequency and ampli- tude, we have to extract them from the estimated spectrum. This is usually realized by 12 peak analysis with a proper selection of parameters such as minimum peak separation, min- imum peak height, etc. However, when handling with realistic data online, those pre-fixed parameters are not suitable for all spectrum shapes, which means not every picked peak corresponds to a natural frequency all the time. Thus, we proposed a novel approach for automatically extracting frequency and amplitude based on the original and the processed spectra. First, it’s found that the spectral lines of APES spectrum (Fig. 3(b)) are redistributed into several groups (Fig. 3(c)) after the processing of FSP, and each group corresponds to a potential location of natural frequency. So, we separate the frequency domain into Psub- intervals. In each sub-interval, the spectral lines are clustered within the maximum distance less than a threshold δΩ, as shown in the following equation: Figure 3: FSP and the extraction of frequency and amplitude. (a) Data series, which is measured from the vibration of a cable of a cable-stayed bridge. (b) APES spectrum of (a), the calculating frequency ω= [0,0.01,0.02,···,25]Hz. (c) FSP result of (b): original spectrum (purple line); spectral lines by FSP (black line); spectral lines by FSP without zero setting of amplitude (blue dashed line). (d) Extracted natural spectral lines by Algorithm 3. The points under the spectral lines in (b) ∼(d) show the distribution of frequencies. 13   Ω ={1,2,···, K}= Ω 1∪ ··· ∪ Ωi∪ ··· ∪ ΩP, Ωi={j| |ˆωj−ˆωi0| ≤δΩ, j∈Ω},where i0=  1, i = 1, max(Ω i−1) + 1, i≥2.(20) Note that the frequency ˆωstill keeps monotonically increasing after FSP, i.e., ˆ ω1≤ˆω2≤ ··· ≤ ˆωK, and the number of sub-intervals P(i.e., the number of extracted frequency and amplitude) is automatically determined. Second, we locate the maximum amplitude of each group Ω i: Ωmax=n mi|mi= arg max m(ˆαm), m∈Ωi, i= 1,2,···, Po , (21) and index miof the maximum amplitude will be regarded as the position of natural frequency if it satisfies 1) ˆ αmi̸= 0 and 2) the distance between original frequency and processed frequency |ˆωmi−ωmi|is less than a threshold δd, i.e., Ωs={i|ˆαi̸= 0 & |ˆωi−ωi|< δd, i∈Ωmax}. (22) The index set of Ω sin Eq.(22) is still just the potential locations where natural frequencies are, because there is no guarantee that the spectral lines will not cluster at the lobes caused by noise, especially in the case of low SNR. Finally, extracting frequencies and amplitudes from the original frequency vector ωand amplitude vector αaccording to the indexes set Ω s:   ωs={ωi|i∈Ωs}, αs={αi|i∈Ωs}.(23) Here, the subscript ’ s’ in Ω s,ωsandαsmeans ’sparse’, because most non-natural frequen- 14 cies and amplitudes have been removed. And we will denote the natural frequencies and amplitudes as ωnandαnin the following, obviously, ’ n’ means ’natural’. If the spectrum enjoys a well shape, each frequency extracted in Eq. (23) will correspond to a natural frequency, as shown in Fig. 3(d). In addition, the frequency deviation in the processed spectrum ( ˆω,ˆα) by FSP as shown in 3(c) is eliminated by Eq. (23), which relocates the biased peaks to the true position in original APES spectrum ( ω,α). Actually, after the processing of FSP, the original spectral lines distributed with an equal spacing of ∆ ω=ωi+1−ωiwill be extremely compact (e.g., at the precision of 0.001Hz or smaller) at local peaks where probably the natural frequencies are. By using this feature, the following manner of cluster separation is suggested which is more efficient than Eq. (20): ˇωi=⌊ˆωi×10n⌋/10n, i= 1,2,···, K, (24a) ´ω= unique( ˇω) ={ˇω1,ˇω2,···,ˇωP}, P≪K, (24b) Ωi={j|ˇωj= ´ωi, j∈Ω}, i= 1,2,···, P. (24c) Eq. (24a) means the frequency ˆ ωiis truncated to the nth decimal place, and ⌊·⌋represents round down. In Eq. (24b), operation unique( ·) means to delete the same elements in ˇω, leaving only the Pdifferent elements. Eq. (24c) divides the original frequency interval into Psub-intervals. The pseudo-code of frequency and amplitude extraction is given in Algorithm 3. Algorithm 3 Automatic frequency and amplitude extraction. 1:Input : Original frequency ωand amplitude αby Algorithm 1, processed frequency ˆω amplitude ˆαby Algorithm 2, threshold δΩor order n, threshold δd 2:Separate the frequency domain into Psub-intervals by Eq. (20) or Eq. (24) 3:Find index of the maximum amplitude of each sub-intervals by Eq. (21) and determine if the index corresponds to a natural frequency by Eq. (22) 4:Extract frequency and amplitude through the index by Eq. (23) 5:Output : Extracted frequency ωsand amplitude αs 15 3.3.Online frequency and amplitude tracking embedded with KF By now, using Algorithms 1 ∼3, it’s able to automatically identify frequency and ampli- tude in a single step. However, when dealing with actual continuous and long-term signals of structural vibrations, it’s hard to guarantee that the vibrations in every computing step are well excited by the environment. The existing signal frames with very low SNR will cause poor identification results mixed with non-natural frequency. In order to improve the naive online extracting results, i.e., each computing frame is processed separately, we apply KF to the tracking framework. As shown in Fig. 2, denote the measured kth (k= 0,1,···) computing frame at time tk:=kS+N−1 asxk:= [xkS, xkS+1,···, xkS+N−1], and denote the procedure of frequency and amplitude extraction from xkusing Algorithms 1 ∼3 as FAE( ·): (ωs,k,αs,k) = FAE( xk), (25) where, ωs,kandαs,krepresent the estimated frequency and amplitude of the kth comput- ing frame, respectively, and the name of operation ’FAE’ means ’frequency and amplitude extraction’. Assume that the extracted frequency and amplitude between adjacent computing frames (Ssamples spaced) do not change very much: ωs,k+1∼=ωs,k,αs,k+1∼=αs,k. (26) Therefore, according to Eq. (2), the vibration signal of ( k+ 1)th computing frame can be written as xk+1,n= Re Pk+1X i=1(αs,k+1,i)ej2π(ωs,k+1,i)n! ∼=Re PkX i=1(αs,k,i)ej2π(ωs,k,i)n! , n= (k+ 1)S,(k+ 1)S+ 1,···,(k+ 1)S+N−1,(27) where ωs,k,iandαs,k,irepresent the ith element of the estimated frequency vector ωs,kand 16 amplitude vector αs,kat the kth computing frame, respectively; Pk≥1 is the number of estimated parameters ωs,kandαs,k; Re(·) returns the real part of a complex number. And we denote the construction operation expressed by Eq. (27) as xk+1= Construct( ωs,k,αs,k). (28) Thus, taking the collected computing frame xkas state vector, the prediction and mea- surement equations can be written as   xk+1=F(xk) +wk, zk =Hxk+vk,(29) where, F:xk→xk+1is realized by Eqs. (25) and (28), i.e., F(xk) = Construct(FAE( xk)); measurement matrix His selected as unit matrix IN×Nto make a full observation; ukis the zero-mean process white noise with covariance matrix E[ ukuT k] =Qk;vkis the zero- mean measurement white noise with covariance matrix E[ vkvT k] =Rk; and E[ ·] stands for expectation. AsFis a nonlinear function, generally, unscented Kalman filter (UKF) [54] or extended Kalman filter (EKF) [55] are resorted to for state estimation. For UKF, sigma points χare generated first[56], and then perform unscented transform on these sigma points χthrough the nonlinear function Fyielding a set of transformed points Y:=F(χ) to compute the mean and covariance of the new estimate. In this procedure, function F composed of the operations of FAE( ·) and Construct( ·) has to be conducted many times, which is extremely time-consuming and is not suitable for online estimate. EKF is more efficient compared to UKF, which handles nonlinearity by locally linearizing the system at the current estimate. To linearize F(xk), we have to calculate its Jacobian matrix at xk, i.e.,∂F(x) ∂x|xk. However, this is unsolved observing the form of function F. Considering that the state vector xk∈RN×1in Eq. (29) represents Nconsecutive observations of a 17 vibration response (e.g. acceleration) at a measurement point of a structure. In this sense, xkis a series of realizations of a single random variable. It’s reasonable to assume that the uncertainty Pof state xkkeeps the same after the mapping of F, which equivalently means F:=∂F(x) ∂x|xk∼=IN×N, (30) and the covariance matrices of process white noise and measurement white noise are reduced to a variance: Qk→Qk:= E[ukTuk],Rk→Rk:= E[vkTvk]. Therefore, a standard liner Kalman filter is formalized as follows: (ωs,k,αs,k) = FAE( ˆxk|k), (31a) ˆxk+1|k=F(ˆxk|k) = Construct( ωs,k,αs,k), (31b) ˜yk+1=zk+1−ˆxk+1|k, (31c) ek+1=∥ˆxk+1|k∥2 2/(∥ˆxk+1|k∥2 2+∥zk+1∥2 2), (31d) dk+1= 1−ˆxT k+1|kzk+1/q (ˆxT k+1|kˆxk+1|k)(zT k+1zk+1), (31e) Qk+1= (ρ(1−ek+1) + (1 −ρ)dk+1)σ˜yk+1, (31f) Rk+1= (ρek+1+ (1−ρ)(1−dk+1))σ˜yk+1, (31g) Pk+1|k=Pk|k+Qk+1, (31h) Kk+1=Pk+1|k/(Pk+1|k+Rk+1), (31i) ˆxk+1|k+1=ˆxk+1|k+Kk+1˜yk+1, (31j) Pk+1|k+1= (1−Kk+1)Pk+1|k. (31k) The process function Fis separated into two parts of Eqs. (31a) and (31b) to explicitly express the extraction of frequencies and amplitudes. Eq. (31c) computes the residual between prediction ˆxk+1|kand measurement zk+1. Eqs. (31d) ∼(31g) aim to computing the variance Qk+1of process noise and the covariance Rk+1of measurement noise adaptively based on the following considerations: 18 (1) If the measurement zk+1enjoys a high energy (defined as ∥zk+1∥2 2), which means the structure is well excitated and the measured signal has a high SNR, thus we should be inclined to trust the measurement to capture changes in the signal’s spectral structure. Otherwise, we should be inclined to the prediction ˆxk+1|kto keep the continuity of identification results. This defines the energy metric eas Eq. (31d). (2) If the prediction ˆxk+1|kis similar to the measurement zk+1, we are inclined to the prediction, because it is reconstructed by extracted frequencies and amplitudes, mean- ing the noise term has been greatly filtered. Otherwise, we should be inclined to the measurement zk+1to ensure the sensitivity to change perception. This defines the distance metric das Eq. (31e). Therefore, Qk+1(process noise’s variance) and Rk+1(measurement noise’s variance) are designed relating to above energy metric eand distance metric das shown in Eq. (31f) and Eq. (31g), where 0 ≤ρ≤1 is a weight between the energy metric and distance metric and σˆyis the variance of residual term ˆy. Eqs. (31h) and (31i) are to obtain the Kalman gain according to the process uncertainty and the measurement uncertainty. Eq. (31j) gives the state estimate and Eq. (31k) gives the estimate of process uncertainty according to Riccati equation [57]. The form of Eq.(31) can be analogous to the standard univariate Kalman filtering. As shown in Eq. (31), we take the structural vibration response at a certain measured point in time-domain as the state variable, and the actually concerned natural frequency and amplitude in frequency-domain are exported from the state estimate. The frequency ωs,kand corresponding amplitude αs,kidentified by Eq. (31a) or Eq. (25) are arranged in order of frequency from the smallest to the largest, i.e., ωs,k,1< ω s,k,2<···< ω s,k,P k. Due to the randomness of environment, there is no guarantee that all modes are excited efficiently always, and we couldn’t make sure that the identified frequency ωs,kand amplitude αs,kat different steps are in alignment. Alignment here has two meanings: (1) the number Pkof identified frequency ωs,kin each step keeps the same: Pk=Pk+1; (2) at the same indexed 19 position i, the element ωs,k,iof vector ωs,kin each step is of the same order. If frequencies in different steps are not aligned all the time, it is hard to conduct further quantitative analysis but for result visualization. Fig. 4 demonstrates the frequency alignment schema designed in this paper. First, the reference frequency ωrand the allowable variation range δahave to be appointed in advance. For each reference position i, select the candidates from extracted frequency ωs,k+1, which are in the supposed interval ( ωr,i−δa,i, ωr,i+δa,i): Ωi={j| |ωs,k+1,j−ωr,i|< δa,i, j∈ {1,2,···, Pk+1}}, i= 1,2,···, Pr. (32) In the above equation, Pk+1means the number of extracted frequency ωs,k+1at the k+ 1th computing frame, Prmeans the number of reference frequency ωr, Ωiis the index set of candidates correspond to ωr,i. Figure 4: Illustration of frequency alignment in Eqs. 32 and 33. ’f & a’ in the figure means ’frequency and amplitude’. Second, determine true natural frequencies and amplitudes from the candidates: if there is no candidate for ωr,i, i.e., Ω i=∅, the frequency will be copied from the last step, and the amplitude will be assigned as half of that in the last step; If the candidates are unique for ωr,i, they will be directly set as the natural frequency and amplitude; If there are multiple candidates for ωr,i, first select the candidate with maximum amplitude, Ω i,max = arg max j{αs,k+1,j, j∈Ωi}, and the candidate that is nearest to the frequency ωn,k,iin the last step, Ω i,nearest = arg min j{|ωs,k+1,j−ωn,k,i|, j∈Ωi}, then if the maximum amplitude 20 αs,k+1,Ωi,max is greater than twice the amplitude αn,k,iin the last step, the frequency and amplitude will be set according to the maximum candidate Ω i,max, otherwise, they are set according to the nearest candidate Ω i,nearst . Those are expressed in the following equation: ωn,k+1,i=  ωn,k,i Ωi=∅, ωs,k+1,Ωi card(Ω i) = 1 ,  ωs,k+1,Ωi,max αs,k+1,Ωi,max≥2αn,k,i ωs,k+1,Ωi,nearest otherwisecard(Ω i)≥2,i= 1,2,···, Pr, (33a) αn,k+1,i=  0.5αn,k,i Ωi=∅, αs,k+1,Ωi card(Ω i) = 1 ,  αs,k+1,Ωi,max αs,k+1,Ωi,max≥2αn,k,i αs,k+1,Ωi,nearest otherwisecard(Ω i)≥2,i= 1,2,···, Pr, (33b) where, card( ·) returns the number of elements in the set. For the situation where the time-varying frequency fluctuates at a certain level, the reference value of frequency can be determined by the statistical average based on the analysis results of historical data, and it can be updated according to the latest analysis results to capture the slow change. For the situation where the frequency shows obvious trend variation, we can set the extracted results in the last step as the reference value of frequency, i.e., ωr←ωn,k. Algorithm 4 gives the pseudo-code of frequency and amplitude tracking embedded with Kalman filtering. 21 Algorithm 4 Frequency and amplitude tracking using APES postprocessing and KL. 1:Input : Initial frequency ωn,0and amplitude αn,0, updating block length S, weight ρ, reference frequency ωrand allowable frequency variation range δa 2:fork= 0,1,2,···do 3: Predict the state ˆxk+1|kby Eq. (31b) 4: Obtain the measurement zk+1= [xkS, xkS+1,···, xN+kS−1] 5: Calculate the residual signal ˜yk+1by Eq. (31c) 6: Calculate the variance of process noise Qk+1and the variance of measurement noise Rk+1by Eqs. (31d) ∼(31g) 7: Calculate the gain Kk+1by Eqs. (31h) ∼(31i) 8: Update the state ˆxk+1|k+1by Eq. (31j) and the process uncertainty Pk+1|k+1by Eq. (31k) 9: Extracted the frequency ωs,k+1and amplitude αs,k+1by Eq. (31a), i.e., Algorithms 1∼3 10: Frequency alignment by Eqs. (32) and (33) to obtain the natural frequency ωn,k+1 and amplitude αn,k+1. 11:end for 12:Output : Time-varying natural frequency [ ωn,0,ωn,1,···] and amplitude [ αn,0,αn,1,···] 3.4.Section summary 3.4.1. Parameter settings in the algorithms The proposed approach involves Algorithms 1 ∼4. A specific description of the parameter settings in those Algorithms is provided in Appendix A. 3.4.2. Usages The proposed method consists of Algorithms 1 ∼4 could be employed for different pur- poses in different ways, in addition to the online frequency and amplitude tracking, which is the most focus of this paper, as illustrated in Fig. 5. 1. Using Algorithms 1 and 2 (APES+FSP) for TFR obtaining. This combination is equivalent to SST or SET in terms of result exhibition, as shown in Fig. 5(a). And this is exactly the work we done in [51]. Because the spectral lines of a TFR slice are concentrated at natural frequencies by frequency squeezing, so we called the result squeezed TFR. 2. Using Algorithms 1 ∼3 (defined as FAE in Eq. (25)) for sparse TFR obtaining. In this combination, the potential spectral lines corresponding to natural frequencies are further extracted on the base of squeezed TFR. That is, the number of spectral lines 22 Alg.1 + Alg.2 +Alg.3 (FAE)(APES+FSP)+Alg.4 (FAE+KL)timepast future Frequency TimeTime Frequency Time Frequency Time (a) (b) (c)squeezed TFR sparse TFR natural TFRFrequencyFrequency FrequencyFigure 5: Usages of the proposed Algorithms. (a) Algorithms 1 and 2 for squeezed TFR. (b) Algorithms 1∼3 for sparse TFR. (c) Algorithms 1 ∼4 for natural frequency and amplitude tracking. of each slice is reduced from KtoPk, and Pk≪K, so we called the result here sparse TFR, as shown in Fig. 5(b). 3. Using Algorithms 1 ∼4 (FAE+KL) for online frequency and amplitude tracking. In which, the Kalman filtering enhances the reliability of sparse TFR and the adopted frequency alignment removes all non-natural spectral lines, as shown in Fig. 5(c). We illustrate the usages of 1 and 2 here is that sometimes the whole TFR plotting gives us more comprehensive and convincing results than the only surviving frequency trajectories. However, the latter supports further quantitative analysis. 3.4.3. Characteristics The main characteristics of the proposed method are described below. 1. The method is applicable to online frequency and amplitude tracking with a small anal- ysis lag and high resolution due to the superiority of APES. And it’s data-dependent 23 in two aspects: one is the data-dependent ’basis’ of APES spectra analysis and the other is the adaptive noise variance setting in Kalman filtering. 2. The framework is easily altered and to derive other specific approaches, such as replac- ing APES with other spectral analysis methods, adopting other frequency alignment strategies, establishing the state equation in Kalman filtering with other forms. 3. The signal reconstruction using the identified frequency and amplitude can be regarded as the filtering of measured signal. 4. The frequency alignment requests that the different modes are well separated, this is true for most civil structures. 4. Numerical verification Because the true values of frequency and amplitude of realistic signals are unknown, before diving to them, we investigate a numerical example to explicitly illustrate the per- formance of the proposed method on accuracy and noise robustness. The considered synthetic signal with a sampling rate of 50Hz consists of three components and an additive Gaussian white noise term: x(t) =x1(t) +x2(t) +x3(t) +n(t),0≤t≤40, x1(t) = 3 cos(2 πt+ 2 sin(0 .25πt)), x 2(t) = 4 e−0.05tcos(5 πt−0.005πt2), x3(t) =  3.0 cos(8 πt),0≤t <10,20≤t <30, 1.5 cos(8 πt),10≤t <20,30≤t <40,n(t)∼ N(0, σ2 n),(34) where σndenotes the standard deviation of the white nose. In this section, a low noise level (σn= 0.5) case in which SNR=21.52dB and a high noise level ( σn= 2.0) case in which SNR=-2.57dB are analyzed. And the SNR is calculated as SNR = 10 log10∥signal∥2 2 ∥noise∥2 2. (35) 24 Fig. 6 plots the components, noise term, synthetic signal and IF laws defined in Eq. (34). The IA functions can be seen from the envelops of components in Figs. 6(a) ∼(c). 010203040 (a)Time/s-4-2024Component, x1 010203040 (b)Time/s-4-2024Component, x2 010203040 (c)Time/s-4-2024Component, x3 010203040 (d)Time/s-10010Noise, n 010203040 (e)Time/s-20020Synthetic signal, x 010203040 (f)Time/s024Frequency/Hzf1f2f3 Figure 6: Numerical signal in Eq. (34). (a) ∼(c) Components 1 ∼3. (d) Noise term (showed σn= 2.0). (e) Synthetic signal. (d) IF law. We compare the proposed method with SET to validate its effectiveness. The frequency and amplitude in SET are obtained by multi-ridge detection algorithm as reported in litera- ture [58]. Meanwhile, we compare the results of proposed method when KF is activated and deactivated to show the performance of KF. Figs.7 and 8 plot the relative error of analysis results under the two noise levels. The relative error is defined as Relative error =|Estimated value −True value | True value×100% . (36) Fig. 7 shows the analysis results corresponding to the case of low Gaussian white noise. In this situation, SET enjoys more accurate frequency and amplitude identifications in most times. For the 1st component, its frequency changes according to the sinusoidal law with a period of 8s. It changes by 1/2 circle in a computing frame (N=200). In this case where the frequency varies rapidly, the assumption of slow change is not satisfied, neither in our method nor in SET, so the errors are relatively high. For the 2nd component, because the SNR decreases with the attenuation of amplitude, the error increases gradually. For the 3rd component, when the amplitude changes suddenly, the amplitude error increases significantly. 25 0 10203040 (a) Time/s01020Error of amplitude (%)1st 0 10203040 (b) Time/s050100Error of amplitude (%)2nd 0 10203040 (c) Time/s0204060Error of amplitude (%)3rdSET FAE without KF FAE with KF 0 10203040 (d) Time/s0510Error of frequency (%)1st 0 10203040 (e) Time/s012Error of frequency (%)2nd 0 10203040 (f) Time/s00.51Error of frequency (%)3rdFigure 7: Analysis results of signal in Eq. (34) with a low noise level ( σn=0.5, SNR=21.52dB). (a) ∼(c) Error of 1st ∼3rd amplitude.(d) ∼(f) Error of 1st ∼3rd frequency. 0 10203040 (a) Time/s050100Error of amplitude (%)1st 0 10203040 (b) Time/s050100150200Error of amplitude (%)2nd 0 10203040 (c) Time/s050100Error of amplitude (%)3rd SET FAE without KF FAE with KF 0 10203040 (d) Time/s05101520Error of frequency (%)1st 0 10203040 (e) Time/s050100Error of frequency (%)2nd 0 10203040 (f) Time/s05101520Error of frequency (%)3rd Figure 8: Analysis results of signal in Eq. (34) with a high noise level ( σn=2.0, SNR=-2.57dB). (a) ∼(c) Error of 1st ∼3rd amplitude. (d) ∼(f) Error of 1st ∼3rd frequency. However, at the high noise level, as shown in Fig. 8, our method exhibits better tracking accuracy and stability. For example, the 2nd and 3rd frequencies obtained by SET are much worse than the proposed method. This can be explained intuitively by the corresponding TF spectrogram in Fig. 9. The TFR obtained by SET are blurrier resulting in the difficulty in ridge detection. While the trajectories in the sparse TFR obtained by FAE (see in Fig.5(b)) are clearer, which is benefit to frequency extraction. What’s more, it’s found that the adoption of KF increase the robustness of the results, e.g., the frequency and amplitude in 26 Figs. 8(a) and (d). Figure 9: (a) The TFR obtained by SET. (b) The sparse TFR obtained by FAE proposed. Although the method of SET can achieve online analysis as the manner we illustrated in Fig. 1. When we dive into more details, the differents between it and the proposed method come out. In the latter, frequency and amplitude are extracted from the spectra by APES and FSP. The signal for spectral analysis every time is Nlength and there is only one frequency and one amplitude to output for each order. While SET belongs to TF analysis methods, the frequency and amplitude are extracted from the spectrogram. The signal for TF analysis every time is N+kS(Nis the window’s size, Sis the moving spacing of window) length and there are ( k+1) frequencies and ( k+1) amplitudes to output for each order. This could be seen as online mini-batch processing of streaming data. Obviously, the method based on spectrum is more real-time with smaller analysis-lag. The values of algorithms’ parameters in the above numerical cases are summarized in Table B.1. 5. Realistic application In this section, we will apply the proposed method to track the natural frequency and amplitude from the vibration signals of three actual structures recorded by their SHM sys- tems. They are the ambient vibrations of a cable of a cable-stayed bridge (Section 5.1), a hanger of a half-through arch bridge (Section 5.2) and the girder of a suspension bridge (Section 5.3), as shown in Fig. 10. In each case, the analyzed signal lasts for 24 hours. The parameter settings in the algorithms are also summarized in Table B.1 and will not be described at the time of discussing these cases in detail. 27 Figure 10: Actual structures that generate the vibrations. (a) A cable of a cable-stayed bridge, span layout: (105+97)m. (b) A hanger of a half-through arch bridge, span layout: (100+450+100)m. (c) The girder of a suspension bridge, span layout: (578+1650)m. 5.1. Cable vibrations of a cable-stayed bridge The vibration signals of a cable of a cable-stayed bridge (Fig. 10(a)) that last from 0:00 at midnight to the next day (Fig. 11(a)) is analyzed. The cable is 102.4357m long with a linear mass of 40.846kg/m, and the cable force measured once is 2496kN. The first 6 orders of frequencies and amplitudes of interests are tracked by the proposed method. Fig. 11(b) plots the time-varying frequencies, and the results during 2.5 ∼2.52h and 8.5∼8.52h are zoomed in Figs. 11(c) and (d). The signals corresponding to these two intervals represent the situations of relatively small and large vibrations, respectively. For the time between 2.5 ∼2.52h at night, there is supposed to be less traffic on the bridge deck than the time between 8.5 ∼8.52h at day, so the signal during the former interval has lower vibration energy. As we can see from Figs. 11(c) and (d), the proposed method is able to track the frequency continuously no matter under the low vibration energy or high. Actually, there still exist instants corresponding to some computing frames, where the shape of amplitude spectrum obtained by APES is not regular enough to exactly extract the natural frequencies and amplitudes of interest. Here a regular spectral shape means that the lobes at all the natural frequencies are higher than the lobes caused by noise. And the non-regular spectral shape is not only going to appear in the case of weak environmental excitation, but also in the case of high environmental excitation because not all modes are excited well. The non-regular spectral shape results in the raw extracted parameters not- aligned and then the frequency-alignment in Eq. (33) is activated to make the identified frequencies aligned and continuous. In the operation of frequency-alignment, if there are no candidates within the allowable variation range of reference value, the frequency in current 28 step is inherited from the last step while the amplitude is taken as half of value in the last step. Thus, in Figs. 11(c) and (d), some straight segments with amplitude close to zero are observable. These frequencies with near zero amplitudes indicate the corresponding modes are not well excited by the ambient. And we lost the tracking of the changing frequency during these periods since they are copied from the last step. Once the mode is excited by the ambient again, the method will recover the tracking of it. Figure 11: (a) Acceleration (Acc.) time history of a cable of a cable-stayed bridge in 24 hours. (b) Tracked first six natural frequencies. (c) Enlargement of (b) during 2.5 ∼2.52h . (d) Enlargement of (b) during 8.5∼8.52h. Furthermore, Figs. 12 and 13 redraw the identified results in Figs. 11(c) and (d) in different panels separately according to the order, which provide a clearer view revealing the frequency fluctuations with respect to magnitudes of the identified amplitude. It’s easy to find that, if the amplitude enjoys a continuous and moderately large magnitude (Figs. 12(b) and 13(b)), the corresponding frequency enjoys a smaller variance. Otherwise, if the amplitude is small and changes quickly (Figs. 12(a), (f) and 13(f)), the corresponding frequency is more discrete. This is reasonable since our method assumes that the signal during the analytical time window of Nlength is stationary. So, the relatively large variation of frequency may be more likely to be caused by the fast or abruptly changing amplitude, not the variation of structural characteristic itself. This also suggests that the frequency 29 results with large amplitudes (i.e., dominant modes) are more reliable. Figure 12: Frequency and amplitude of each order corresponding to the result in Fig. 11(c). (a) ∼(f) Results from 1st to 6th order. Figure 13: Frequency and amplitude of each order corresponding to the result in Fig. 11(d). (a) ∼(f) Results from 1st to 6th order. 5.2. Hanger vibrations of a half-through arch bridge In the last subsection, we have discussed the results of a cable vibration in detail. Here we give the analysis results of the 24-hour vibrations of a hanger of a half-through arch bridge (Fig. 10(b)) in brief since the major conclusions are similar. The hanger is 41.2269m long with a linear mass of 14.0kg/m, and the cable force measured once is 666kN. Fig. 14(a) is the time history of the hanger vibrations and Fig. 14(b) shows the obtained frequencies of the first 4 modes. For a clearer view of the frequency changes, an enlargement of Fig. 14(b) in the time interval between 2.5h to 2.52h is given in Fig. 14(c). What’s more, as what Figs. 12 and 13 do, Fig. 14(d) plots the frequency as well as the amplitude of each order separately, note that the limitations of the frequency and amplitude axes of each panel have been set to be consistent for comparison purpose. 30 Figure 14: (a) Acceleration time history of a hanger of a half-through arch bridge in 24 hours. (b) Tracked first four natural frequencies. (c) Enlargement of (b) during 2.5 ∼2.52h. (d) Frequency and amplitude of each order corresponding to the result in (c). As we can see, for this about 41m hanger, which is more well excited than the cable in the last subsection (compare Fig. 14(a) with 11(a)), showing more convincing frequency identification results (compare Fig. 14(c) with 11(c)). And Fig. 14(d) provides more proof that the frequency variation may also be caused by the rapidly varying amplitude not only by the real variation of itself. And the former is attributed to the assumption that the amplitude is slowly varying in an analytical window used in spectral analysis. And this assumption is common in many existing methods no matter whether to mode decomposition or TFA. 5.3. Girder vibrations of a suspension bridge Lastly, we analyze the vibration signals generated by the girder of a suspension bridge, and the signals are recorded by an accelerometer installed at 1/4 of the main span, as shown in Fig. 10(c). Compared with the frequency and amplitude tracking of the cable or hanger, whose vibrational modes are well separated by almost equal distances, the same objective for the long-span suspension bridge becomes more challenging because its dominated modes are concentrated in the low frequency region and it’s possible to happen mode crossover. Similarly, the acceleration time history of the main girder is depicted in Fig. 15(a). And Fig. 15(b) is a scatter graph of the tracked results, showing the frequencies at different 31 instants and the brighter color means larger amplitude corresponded. As we can see, the suspension bridge exhibits 8 close vertical modes under 0.5Hz, and the frequency trajectories fluctuate around the early measured natural frequencies, which can be obtained from the Table 1 of the paper by H. Li et, al. [59]. Also, we list the 8 measured natural frequencies here. They are 0.0953, 0.1328, 0.1825, 0.2301, 0.2757, 0.3237, 0.4351 and 0.4918Hz. Observe the interval marked in Fig. 15, it is found that the bridge experienced vortex-induced vibrations (VIVs). The occurrence of VIV can be observed from both the acceleration time history in Fig. 15(a) and the extracted frequencies and amplitudes in Fig. 15(b). And the latter gives us more information about the energy distribution among different modes and which mode is participating in the VIV. Figure 15: (a) Acceleration time history of the girder of a suspension bridge in 24 hours. (b) Tracked results in the low frequency zone. We enlarge the results during the VIV as shown in Fig. 16(a) and plot the time-varying amplitude in Fig. 16(b). It’s interesting to find that there exists an interchange for the modes participating in the VIV. At the beginning, the vibration energy is concentrated on the 6th mode, then the 6th mode’s energy is progressively reduced while the 5th mode’s energy is steadily raised until the 5th mode becomes the dominant one. Similarly, at the end of the disappearance of VIV on the 5th mode, the vibration energy is concentrated 32 on the 6th mode again for a short while. This should be the product of the interactions between structure and wind. Since VIV is generated when the frequency of vortex shedding caused by wind flow passing through the bridge section is close to a certain natural frequency of the structure under specific wind field conditions, this mode conversion may be due to the change of vortex shedding frequency caused by the change of wind field conditions. Additionally, from Fig. 16, we can determine when the VIV happens, which mode is most exciting and how long the VIV lasts. These discoveries are going to offer a new fundamental understanding of VIV of large-scale suspension bridge for the monitoring, perception and vibration control. Of course, this is beyond the scope of this paper and will be one of the focuses of our future work. Figure 16: (a) Enlargement of the results in Fig. 16(b) during 18 ∼23h where VIVs occurred. (b) The tracked amplitudes corresponding to the first 8 mode during 18 ∼23h. Furthermore, we pull our eyes back to the tracking results of the time-varying frequency (Figs. 15(b) and 15(a)). The proposed method is able to obtain the time series of frequency and amplitude of each mode. However, for this main girder with close modes in the low frequency zone, under environmental excitation, the mode mixing or coinciding are inevitable in the identified frequencies. Thus, the frequency trajectories in Figs. 15(b) and 15(a) at some intervals appear turbulent. On the other hand, to what extent this ’turbulence’ is a 33 true reflection of structural characteristics is worth further investigation. This phenomenon is slighter in the cases of cable and hanger discussed above because their modes are more separated. 6. Conclusion For the purpose of tracking natural frequency and amplitude of civil structures from their monitoring data, this paper develops a method based on APES spectrum postprocessing and KF. Its effectiveness is validated through the analysis of numerical examples and the comparison with SET, exhibiting more accurate and robust results in the case of low SNR. What’s more, we employ the proposed method to analyze the vibration signals generated by actual engineering structures. While proving its ability to extract frequency and amplitude online in environmental vibrations, some main conclusions are reached as follows: 1. Under the random ambient excitation, the magnitude and changing rate of ampli- tude of different modes are quite different. For modes with large and slowly changed amplitude, the variance of identified frequency is small (Fig. 13(b)), and for modes with small and fast changed amplitude (not well excited), the variance of identified frequency is large (Fig. 13(f)). This is because the severe amplitude modulation is going to submerge the information of structural oscillations and bring errors in the identified frequency positions. On the other hand, those frequencies with relatively large and slow-changing amplitudes are more reliable. In other words, the identified amplitude can be used as a measure of the reliability of identified frequency. 2. Under environmental excitation, it is uncertain which modes are well excited and the SNR will be quite low sometimes resulting the characteristic information of vibrations immerged in noises. This brings challenges to the online tracking of frequency and amplitude. The proposed method utilizes KF and frequency alignment scheme to keep the robustness of the tracking results and it shows positive performance in long-term analysis. However, for the modes that are very closely spaced, the tracking results may be mixed in the case of large frequency variation is caused by the varying ambient. 34 3. The frequency alignment scheme provided in this paper is adaptive to separated modes, which is true for most engineering structures. When handling with the situation of mode crossover, it has to be newly designed. Our future work will be focused on the more applications of the proposed method, such as the damage diagnosis and abnormal event detection, e.g., the VIV monitoring of suspension bridge discussed in the Section 5.3. Appendix A. Discussion on parameter setting in Algorithms 1 ∼4 Appendix A.1. Parameters in Algorithm 1: APES a. Computing frame length Nshould be selected according to the longest period of signal components, and Nis recommended to be 5 to 10 times the longest period. b. Calculating frequency ω= [ω1, ω2,···, ωK]Tis assigned according to frequency in- terval of interest, and 0 ≤ωi≤Fs/2. Generally, the elements of frequency vector ωare uniformly distributed, and the spacing ∆ ω:=ωk+1−ωkis set according the requirement of precision, such as 0.01Hz, 0.001Hz. c. Filter order Mshould be selected as a trade-off between the resolution and the statis- tical stability, and it is suggested N/4≤M≤N/2 [52, 60]. Appendix A.2. Parameters in Algorithm 2: FSP a. Order m(Eq. (17)) is used to take the mth power of the normalized amplitude, which is beneficial to promote convergence and effectiveness of frequency squeezing, and is suggested as 30 ∼60. b. Step size P(Eq. (18)) is a key parameter for frequency squeezing. A large value ofPis helpful for spectrum concentration. However, a too large Pwill lead to the missing of natural frequency with relatively small amplitude. Generally, the frequency span composed of the 2 P+ 1 spectral lines shouldn’t include more than 2 natural frequencies. 35 c. Threshold δω(Eq. (19)) is used as a criterion to judge the edges of frequency clusters. Based on the observation that frequency points will be extremely gathered at the natural frequencies after several steps of frequency squeezing (Eq. 18), a small δω should be selected and we recommend it as δω= 0.01∆ω, where ∆ ω:=ωi+1−ωi. d. Threshold δsis used to stop iteration of frequency squeezing, a small value such as 1e-4∼1e-8 could be set. e. Maximum iterations Iis for compulsively stopping the iteration in case the frequency squeezing fails to converge in finite steps. In this paper, it is set as I= 50. Appendix A.3. Parameters in Algorithm 3: Automatic frequency and amplitude extraction a. Threshold δΩ(Eq. (20)) is used to separate the processed frequency ˆωintoPsub- intervals in which the maximum distance between individuals is less than δΩ. And we recommend it as Fs/N, which is the precision of fast Fourier Transform (FFT). b. Order n(Eq. (24a)) is used to truncate decimal places, and it’s suggested as n= max(1 ,⌊log10(1/∆ω)⌋). c.δd(Eq. (22)) determines a threshold to judge if the frequency deviation caused by FSP is acceptable, and it’s also suggested to be Fs/N. Appendix A.4. Parameters in Algorithm 4: Frequency and amplitude tracking using APES postprocessing and KF a. Initial frequency ωn,0and amplitude αn,0can be obtained from historical data, and ωn,0can also be assigned as the reference frequency ωr. b. Updating block length Scontrols the output frequency of results, i.e., a smaller S means a higher time resolution. Meanwhile, the computations in a single step should be completed within Ssamples, so Scan’t be too small (see in Eq. 8). And we recommended it as 1/10 to 1/5 times the computing frame length N. c. Weight ρ(Eqs. (31f), (31g) balances the energy metric and distance metric for updat- ing covariances of process noise and measurement noise. If we emphasize the energy 36 metric (more faith in prediction), set a large ρ; if we emphasize the distance metric (more faith in measurement), set a small ρ. And it’s 0.5 used in this paper. d.δais the allowable range of frequency variation with respect to reference value ωr. It could be a vector with the same length as the reference frequency vector ωr, that is, each order of reference frequency ωr,ihas its own allowable variation range; or it could be a scalar, that is, all order reference frequencies have the same allowable variation range. Appendix B. Specific parameter setting in the numerical and realistic cases Table B.1: Parameter setting of the proposed method in analyzing numerical and realistic cases. Alg. Para. Sec. 4 Sec. 5.1 Sec. 5.2 Sec. 5.3 1N 200 150 100 4000 ω [0,0.01,···,5] [0 ,0.01,···,8] [0 ,0.01,···,12] [0 ,0.001,···,0.5] M 33 37 37 666 Fs 50 50 50 50 2m 50 50 50 50 P 25 60 100 20 δω 1e-4 1e-4 1e-4 1e-5 δs 1e-8 1e-6 1e-6 1e-6 I 50 50 50 50 3n 2 2 2 3 δΩ 0.25 0.3333 0.5 0.0125 δd 0.25 0.3333 0.5 0.0125 4S 10 15 10 100 ωr ωn,k 1.25[1,2,3,4,5,6] 2.50[1,2,3,4][.095,.133,.183,.230, .275,.324,.435,.492] ωn,0 [1.1,2.49,4.0] 1.25[1,2,3,4,5,6] 2.50[1,2,3,4][.095,.133,.183,.230, .275,.324,.435,.492] αn,0 [3.0,3.71,3.0] [1,1,1,1,1,1] [1,1,1,1] [1,1,1,1,1,1,1,1] ρ 0.5 0.5 0.5 0.5 δa 0.1 0.3125 0.625 0.0238 Acknowledgments This work is supported by the Opening Foundation of National Key Laboratory of Bridge Structural Health and Safety; the National Nature Science Foundation of China (Grant No. 5187849); the National Key R&D Program of China (2017YFF0205605). 37