Natural frequencies and modal damping ratios identification of civil structures from ambient vibration data

Damping is a mechanism that dissipates vibration energy in dynamic systems and plays a key role in dynamic response prediction, vibration control as well as in structural health monitoring during service. In this paper a time domain and a time-scale domain approaches are used for damping estimation of engineering structures, using ambient response data only. The use of tests under ambient vibration is increasingly popular today because they allow to measure the structural response in service. In this paper we consider two engineering structures excited by ambient forces. The first structure is the 310 m tall TV tower recently constructed in the city of Nanjing in China. The second example concerns the Jinma cable-stayed bridge that connects Guangzhou and Zhaoqing in China. It is a single tower, double row cable-stayed bridge supported by 112 stay cables. Ambient vibration of each cable is carried out using accelerometers. From output data only, the modal parameter are extracted using a subspace method and the wavelet transform method.


Introduction
Modal parameters identified from the measured data can reflect the dynamic characteristics of a vibrating system and often serve as input to model updating, health monitor or damage diagnosis.To identify the modal parameters of a system, both the excitation force and the response of the system should be measured and the frequency response function between the excitation and the response need to be calculated.In laboratory cases this experimental modal analysis method is effective and practical.However, for large and heavy structures, expensive and very strong exciting equipment is needed and the testing process may damage the mechanical structure.The last ten years, more attention was paid to so-called ambient excitation structures.The structural response is measured to freely available natural forces such as wind, traffic, waves and micro earthquakes.Obviously, the excitation force cannot be measured and is not available as input measurement.The mechanical system is subject to an uncontrolled, unmeasured and non-stationary excitation.The advantage of using ambient forces is that they are more representative for the true excitation to which a structure is subjected during its lifetime.Hence, the need arises to identify modal parameters in real operational conditions, from output only measurements.In this paper the modal parameters of structures are extracted from output-data only using subspace methods [1][2][3] and the wavelet transform [4,5].Two examples using real data are presented.The first experiment concerns the TV tower in the city of Nanjing in China [6].This tower is 310 m high and the acceleration response of the structural system measured under ambient conditions is used to identify its dynamic characteristics.The second experiment concerns the Jinma cable-stayed bridge that connects Guangzhou and Zhaoqing in China [7].It is a single tower, double row cable-stayed bridge supported by 112 stay cables.Ambient vibration of each stay cable is carried out using accelerometers.The identification of modal parameters is compared using the subspace method and the wavelet transform technique.

The stochastic state-space model and the modal model
The subspace identification method assumes that the dynamic behaviour of a structure excited by ambient forces can be described by a stochastic state space model [1][2][3]: where y k is the (mx1) vector of observations, w k , v k are white noise terms representing process noise and measurement noise together with the unknown inputs; z k is the (nx1) internal state vector; A is the (nxn) state matrix describing the dynamics of the system and C is the (mxn) output matrix, translating the internal state of the system into observations.The modal parameters of a vibrating system are obtained by applying the eigenvalue decomposition of A: where Λ = diag(λ i ), i = 1, 2, . . ., n, is the diagonal matrix containing the discrete-time complex eigenvalues and Ψ contains the eigenvectors as columns.The eigenfrequencies f i and damping ratios ζ i are obtained from the eigenvalues: The mode shapes are obtained by pre-multiplying the eigenvectors with C: Our purpose is to determine the state matrix A and the output matrix C in order to obtain the modal parameters of the vibrating system.

Estimation of the state matrix A
We define the (mf x1) and (mpx1) future and past data vectors as and the (mf xmp) covariance matrix between the future and the past as: where E[•] denotes the expectation operator, H is the block Hankel matrix (a block band counter diagonal matrix) formed with the (mxm) theoretical covariance matrices: In practice we estimate sample covariance matrices from data which are used to estimate a sample block Hankel matrix Ĥ.
In order to obtain the state matrix A, two matrix decompositions of the block Hankel matrix are employed: the singular value decomposition and the factorization of Ĥ into its observability and controllability matrices.The singular value decomposition of Ĥ is Ĥ = ÛŜ VT with Û and V orthogonal matrices and Ŝ a diagonal matrix of singular values.The second factorization of the block Hankel matrix into its (mf xn) observability and (nxmp) controllability matrices, O and K is: The two factorisations of the estimated block Hankel matrix can be equated to give: implying: Two methods are used to estimate the state matrix A: the shifted block Hankel matrix method and the observability matrix method.

-The shifted block Hankel matrix method
To estimate A, it is necessary to introduce a time shifted block Hankel matrix, ← H, obtained by shifting a block column or a block row: We obtain: The state matrix A is obtained by applying the pseudo inverses of O and K, yielding: where (•) # represents the pseudo inverse of a matrix.

-The observability matrix method:
We introduce two matrices: where O ↓ (m(f − 1)xn) is the matrix obtained by deleting the last block row of O and O ↑ (m(f − 1)xn) is the matrix obtained by deleting the first block row of O.We obtain then: Finally:

Estimation of the output matrix C
Two methods are used to determine the output matrix C: the Hankel matrix method and the observability matrix method.

-The block Hankel matrix method
Consider H 1L , the first block row of the block Hankel matrix H, we have H 1L = CK.We obtain -The observability matrix method It is easy to show that C is the first block row of the observability matrix O: With the determination of the output matrix C and the eigenvector matrix Ψ, we obtain the mode shapes of the structure using Eq. ( 5).

Identification of modal parameters: practical aspect
When performing modal analysis one of the key decisions that must be upon starting to analyse the data is to decide how many modes there are in the frequency range of interest.This can be done with the stabilization diagram which is developed to track the estimates of eigenfrequencies and damping ratios, in the frequency range of interest, as a function of the model order (the size of the transition matrix A).As the model order is increased, more and more frequencies and damping ratios are estimated, but hopefully, the estimates of the physical modal parameters stabilize as the correct model order is found.For modes that are very excited in the measured data, the modal parameters stabilize at a very low model order.For modes that are poorly excited in the measured data, the modal parameters may not stabilize until a very high model order is chosen.Nevertheless, the non physical modes (in generally computational modes, or modes due to noise) do not stabilize at all during the process and can be sorted out of the modal parameter data set more easily.Note that inconsistencies (frequency shifts, leakage errors . . . ) in the measured data set obscure slightly the stabilization diagram and a tolerance in percentage is given for the stability of each of the modal parameters that are being evaluated.
In practice, we progressively increase the model order n, determine the state matrix A and form its eigenvalue decomposition to obtain the eigenfrequencies and damping ratios of the system for each n.We establish then the stabilization diagram.The output matrix C is also determined to obtain the mode shapes.In order to select the physical modes and to suppress the spurious modes, eigenfrequencies, damping ratios and mode shape vectors between two successive orders are compared.The criteria are 1% for eigenfrequencies, 10% for damping ratios and 99.5% for mode shape vectors (MAC).
Experimental results showing the effectiveness of the procedure are presented in Section 4.

Definitions
The wavelet transform (WT) gives time and frequency information about the analyzed data.The wavelet transform of a signal x(t) is defined as [8,9] where ψ(t) is an analyzing function called mother wavelet, a is the dilatation or scale parameter defining the analysing window stretching and b is the translation parameter localising the wavelet function in the time domain.
The WT represents the correlation between the signal x(t) and a scaled version of the function ψ(t).The idea of the WT is to decompose a signal x(t) into wavelet coefficients using the basis of wavelet functions.The decomposition is obtained locally at different time windows and frequency bands.The size of the time window is controlled by the translation parameter b while the length of the frequency band is controlled by the dilatation parameter a.Hence, one can examine the signal at different time windows and frequency bands by controlling translation and dilatation.
Since the wavelet transform is a linear representation of a signal, it follows that the WT of P signals is this property is convenient for the analysis of multi-component signals.
A number of different analyzing functions have been used in the wavelet analysis.One of the most known and widely used is the Morlet wavelet defined in the time domain as [4,5,8] where ω o is the wavelet frequency.The dilated version of the Fourier transform is In practice the value of ω o is chosen ω o 5 which meets approximately the requirements of admissibility condition [9].Note that ψ(aω) is maximum at the central frequency ω c = ω o /a and the Morlet wavelet can be viewed as a linear bandpass filter whose bandwidth is proportional to 1/a or to the central frequency.Thus, the value of the dilatation parameter a at which the wavelet filter is focused on the wavelet frequency can be determined from a = ω o /ω c .
For a given value of the dilatation parameter a the spectrum of the Morlet wavelet has a fixed bandwidth.If the analysed frequency is important the dilatation parameter becomes small and the spectrum of the Morlet wavelet function is wide.There is then a bad spectral resolution.An alternative is to modify the Morlet wavelet function by introduction of a parameter N which controls the shape of the basic wavelet: this parameter balances the time resolution and frequency resolution of the Morlet wavelet.The modified Morlet wavelet function used in this paper is: with N > 0 and whose dilated version of its Fourier transform is: The wavelet filter central frequency is ω c = ω o /a and gives then a relation between the scale parameter a and the central frequency of the modified Morlet wavelet.
The duration and bandwidth of the modified Morlet wavelet function are given by: When N tends to infinity the modified Morlet wavelet tends to e jωot which has the finest frequency resolution allowing a better resolution of closely spaced modes, but at the expense of time resolution.Indeed, increasing N will increase the frequency resolution but it decreases the time resolution.So, there always exists an optimal value of N that has the best time-frequency resolution for a certain signal localized in the time-frequency plane.This modified Morlet wavelet function offers a better compromise in terms of localization, in both time and frequency for a signal, than the traditionally Morlet wavelet function.The optimal value of N is obtained by minimizing the entropy of the wavelet transform [10].

Application of the wavelet transform to modulated signals
Consider the case of a signal x(t) modulated in amplitude: if x(t) is assumed to be asymptotic, the WT of x(t) can be obtained by means of asymptotic techniques and can be expressed as [4,5] W the point indicating a derivative.The dilatation parameter can be calculated in order to maximize ψ * (a φ(b)), that is using the Morlet wavelet (or the modified Morlet wavelet) for the dilatation a(b) = ω o / φ(b).The wavelet transform is essentially concentrated in a neighbourhood of a curve given by a(b).
Consider now the free response of a viscously damped single degree of freedom system: with B the residue magnitude, ω n the undamped natural frequency, ω d = ω n 1 − ζ 2 the damped natural frequency, ζ the viscous damping ratio.If the system is underdamped, that is if the damping ratio is smaller than 1, (in general 0 < ζ 1, so ω d ; ω n ) the signal x(t) can be considered asymptotic, and therefore the results obtained previously can be used considering: The wavelet transform of the damped sinusoid is: For a fixed value a o of the dilatation parameter the logarithm of the wavelet transform amplitude is: Thus the decay rate σ = ω n ζ of the signal can be estimated from the slope of the straight line of the logarithm of the wavelet transform modulus.The wavelet transform phase is given by: and the plot of d db Arg (W ψ [x] (a o , b)) should be constant in time and equal to the damped natural frequency ω d .The damping ratio and eigenfrequency estimation procedures, based on the wavelet transform presented above, can be extended to multi degrees of freedom systems by selecting the right value of the dilatation parameter a i corresponding to the mode of interest.Consider now the free response of a P degrees of freedom system: where B i is the residue magnitude, ζ i is the damping ratio, ω ni the undamped natural frequency and ω di the damped natural frequency associated to the i-th mode.From Eqs (23) and (34), the wavelet transform of the multi degrees of freedom system is:  The wavelet transform is a signal decomposition procedure working as a filter in the time-frequency domain: it analyzes a signal only locally at windows defined by the wavelet.Thus, multi-degrees of freedom system can be decoupled into single degree of freedom.For a fixed value of the dilatation parameter (a = a i ), which maximizes ψ * (aω di ), only the mode associated with a i gives a relevant contribution in the wavelet transform, while all the other terms are negligible.Thus the wavelet transform of each separated mode i = 1, 2, . . ., P becomes: Clearly, the wavelet transform offers a decoupling of multi degrees of freedom systems into single modes.However, Eq. ( 39) is true under the assumption of vanishing ψ * (a i ω di ) outside the interval that is, if none of the other frequencies of the system except ω i and more likely if neither ω i−1 or ω i+1 belongs to the interval The resolution of the wavelet transform is good enough to separate the i-th mode from the neighbouring modes.
Using Eq. ( 39) associated with Eqs ( 35) and ( 36), it is possible to follow the amplitude and phase variations in the time domain of each modal component and to estimate the corresponding damping ratio and eigenfrequency associated to the isolated mode.This technique requires a previous choice of the value of the dilatation parameter a i corresponding to the analyzed mode and the resolution of the wavelet transform depends on the value of this scale parameter, thus the choice of the analyzing wavelet is important.

The TV tower of Nanjing in China
Figure 1(a) shows the main structure of the TV tower in the city of Nanjing in China [6].This tower is 310 m high and the acceleration response of the structural system measured under ambient conditions is used to identify its dynamic characteristics.The accelerometers are installed on the tower at two sets of different locations, as shown in Fig. 1(a), to measure the ambient vibrations of the system.The sensors at the first set of locations are concentrated on the upper part of the structure since this part is more flexible resulting in more vibration than the lower part, while those at the second set of locations are distributed along the height of the whole structure.The accelerometers are placed as close as possible to the centre of the cross section of the tower in order to minimise the response component due to torsional vibration.Acceleration records are obtained simultaneously in one direction each time, with a sampling time interval of 0.03125 second and a total recording time of approximately 600 seconds.
Figure 1(b) shows the stabilization diagram of the tower using the observability matrix method.From this plot we extract eigenfrequencies and damping ratios of the tower.Only accelerometer number two on the top of the tower has been used, this aceelerometer being very excited.To eliminate spurious poles we use Eq. ( 21).The identified modal parameters for the five first modes are presented in the Table 1.
The WT estimation technique operates on the free response of the system.A method converting random responses to free decay responses is the random decrement technique [12].Its basic concept is that the acceleration response can be decomposed into free vibration component and forced vibration component.The free vibration component can be obtained by a special averaging procedure of measurements which remove the random part, leaving its deterministic part.We apply the random decrement technique to the accelerometer number 2 to obtain its free response (see Fig. 2).We apply then the WT method.Figure 3 shows the wavelet transform amplitude using accelerometer number 2; the five first modes are visible from this plot and in particular from Fig. 3(b).The modal parameters are obtained using Eqs (35) and (36) and from the plots shown in Fig. 4.
Table 1 presents also the modal parameters obtained using this accelerometer by application of the wavelet transform.In generally the results are similar to those obtained by the subspace method.However the third mode presents a damping ratio slightly different.The difference comes from spurious poles which appear with the subspace identification method and are not eliminated automatically.

The Jinma cable-stayed bridge
The subspace method and the wavelet transform method were applied to the analysis of stay cables of the Jinma cable-stayed bridge (Fig. 5a, b), that connects Guangzhou and Zhaoqing in Guangdong Province, China.It is a  single tower, double row cable-stayed bridge, supported by 28*4 = 112 stay cables.The stay cables were excited by ambient vibrations essentially due to wind.Inputs could evidently not be measured, so only acceleration data are available.A full description of the test can be found in [7].
The subspace identification method is used to obtain modal parameters of a cable (cable number 10) and they are presented in Table 2.The instability is assessed by the Scruton number [11] defined as S c,i = ζi.ρρaD 2 , where ζ i is the damping ratio for each mode, ρ is the mass of the cable per meter, ρ a is the density of the air and D is the cable diameter.High values of the Scruton number tend to suppress the oscillation and bring up the start of instability at high wind speeds.Considering ρ = 66.94 kg/m, ρ a = 1.2 kg/m 3 and D = 0.203 m, the Scruton number for each mode is presented in Table 2.
The identification procedure of modal parameters using the WT is carried out as previously.We apply the random decrement technique to the cable number 10 to obtain its free response (see Fig. 6).We apply then the WT method.
The wavelet transform amplitude is shown in Fig. 7(a) and the dilatation parameter a i for each eigen-mode is obtained from Fig. 7(b).The estimated modal parameters for the eight first modes of the cable number 10 are shown in Table 3, with the Scruton number for each mode.
The fundamental frequency of the cable number 10 has been obtained using the subspace method and the WT method (Fig. 8(a)).We obtain approximately the same value: f 0 = 0.697 Hz.This value is obtained from the slope of the straight line.The cable tension can be estimated by the approximated expression T = 4ρL 2 f 2 0 , where L = 180,174 m, ρ = 66.94 kg/m.We obtain T = 4227,593 kN.This cable tension can be considered as a reference and used as an indicator for monitoring of stay cables in this cable-stayed bridge.shows fundamental frequencies of all stay cables (on upstream and downstream side).These frequencies vary between 0.533 Hz for the longest cable and 2.703 Hz for the shortest cable.The maximum and minimum cable forces for the Jinma bridge are then: T max = 5052,073 kN (cable number 57), T min = 2490,653 kN (cable number 84).These cable forces can be considered as reference tensions and used as indicators in the field of health monitoring process.

Conclusion
It is shown in the paper that the subspace method and the wavelet transform method can be effectively employed in operational modal analysis.The results demonstrate that the automatic estraction of damped natural frequencies and modal damping ratios can be successfully performed from ambient vibration data and eventually used as a non destructive health monitoring technique.The eigenfrequencies and damping ratios of a tower have been extracted and can be used as reference for studying the safety comportment of this tower.A cable-stayed bridge has been studied  and stay cables eigenfrequencies and damping ratios could be used to assess the health of cables of this bridge.It is shown how operational modal analysis applied to the dynamic data of stay cables provide useful information to determine cable force and the current condition of stay cables accurately.A Scruton number has been computed and used as an indicator to prevent rain-wind induced vibration.A monitoring system which records the vibrations of stay cables and a software that extracts automatically the force cables and the Scruton number is under investigation.
Large structures tend to present large motions, and therefore nonlinearities.The use of the wavelet transform for the identification of nonlinearities on damping and stiffness is under investigation.However, the concept established in the paper, where the hypothesis of linear vibrations is used, can be extended to nonlinear vibrations [13].Indeed, if we consider the response of a nonlinear oscillator and apply the wavelet transform, it is shown [13] that the amplitude

Figure 8 (
Figure 8(b)  shows fundamental frequencies of all stay cables (on upstream and downstream side).These frequencies vary between 0.533 Hz for the longest cable and 2.703 Hz for the shortest cable.The maximum and minimum cable forces for the Jinma bridge are then: T max = 5052,073 kN (cable number 57), T min = 2490,653 kN (cable number 84).These cable forces can be considered as reference tensions and used as indicators in the field of health monitoring process.

Fig. 7 .
Fig. 7. (a) Wavelet transform amplitude; (b) Determination of dilatation parameters a i .of the wavelet transform is related to nonlinear damping coefficients and the phase of the wavelet transform is related to to nonlinear stiffness parameters.From the amplitude and the phase of the wavelet transform we can identify nonlinear damping and nonlinear stiffness in vibrating systems.

Fig. 8 .
Fig. 8. (a) Fundamental cable frequency; (b) Fundamental cable frequency of all stay cables for the Jinma Bridge.

Table 1
Natural frequencies and damping ratios using the subspace method and the WT method

Table 2
Natural frequencies, damping ratios and Scruton number using the subspace method

Table 3
Natural frequencies, damping ratios and Scruton number using the WT method