Online Modal Identification of Concrete Dams Using the Subspace Tracking-Based Method

To investigate the time-varying dynamic characteristics of concrete dams under the excitation of large earthquakes for online structural health monitoring and damage evaluation, an online modal identification procedure based on strong-motion records is proposed. The online modal identification of concrete dams is expressed as a subspace tracking problem, and a newly developed recursive stochastic subspace identification (RSSI) method based on the generalized yet another subspace tracker (GYAST) algorithm, which exploits both the accuracy of the subspace identification and fast computational capability, is used to extract the time-varying modal parameters of concrete dams during earthquakes. With the simulated vibration response records, a numerical example is used to verify the accuracy, robustness, and efficiency of the proposed GYAST-based, time-varying modal identification method. Then, the realistic strong-motion records of the Pacoima arch dam are analysed using the proposed modal identification procedure, and the time-varying characteristics of the concrete arch dam during three different earthquakes are analysed.


Introduction
e modal parameters of concrete dams, including the natural frequencies, damping ratios, and modal shape vectors, extracted from earthquake response records can be used for structural antiearthquake capacity analysis, health monitoring, and postearthquake structural damage evaluation. ese parameters have obvious physical meanings and represent the practical dynamic characteristics of structures. Some researchers have made efforts to address the modal identification problem of concrete dams based on earthquake response records [1][2][3][4][5][6]. However, the modal identification methods adopted in the past studies are restricted to linear systems and stationary processes, and using such methods would require the assumption that the damreservoir-foundation system is linear and time invariant for the duration of a vibration measurement. e damreservoir-foundation system under the excitation of large earthquakes may present some time-varying characteristics, resulting from the time-varying character of the material parameters, boundary conditions, dam-reservoirfoundation interactions, and structural damage [7]. Using the offline modal identification procedure to extract the constant modal parameters cannot track the time-varying dynamic properties of a structure during a large earthquake, which may be important information for online structural health monitoring. erefore, the online modal identification method of concrete dams is an important research topic. e structural online modal identification problem is a research topic with many challenges that still attract significant attention for its significance in aerospace engineering [8], mechanical engineering [9], and civil engineering [10], among other fields. For time-varying structures, the assumption of the classic modal identification methods that the structure is a time-invariant system is no longer applicable. us, Liu [11] proposed the concept of "pseudomodal parameters," which use the time-freezing concept. ree types of methods have been developed for time-varying modal identification, i.e., the time series model method and the recursive stochastic subspace identification (RSSI) method in the time domain and the time-frequency decomposition method in the time-frequency domain [12,13]. For the time series model method, Gong [14] used adaptive forgetting through multiple models (AFMM) to track the time-varying modal parameters of tall buildings during a large earthquake and evaluated the structural health according to the modal identification results. Poulimenos and Fassois [15] performed an important literature review of different time-varying modal identification methods using different nonstationary time series models and point out their pros and cons. Klepka and Uhl [16] used the timefrequency domain decomposition technique to analyse the structural vibration response and identified the modal parameters. Goethals et al. [17] and Nezam Sarmadi and Venkatasubramanian [18] used the projection approximation subspace tracking (PAST) algorithm-based RSSI (PAST-RSSI), and Loh et al. [19,20] used the extended instrumental variable version of the PAST algorithm-based RSSI (EIV-PAST-RSSI) to carry out the online modal identification of a power system, an in-flight flutter, a tall building, and a television tower. Among these methods, the RSSI method has attracted the most attention. e advantage of the RSSI method is mainly reflected in its strong numerical stability and strong robustness to noise interference. erefore, it is suitable for the online modal identification of structures under random vibration excitation such as earthquakes. Although the online modal identification problem has been studied in many other fields, very few research works have been reported that use the online modal identification for concrete dams [21]. e dam-reservoirfoundation system is characterized by its large size, large number of degrees of freedom, and closely spaced modes, and the vibration measurement of concrete dams is usually performed with a low signal-to-noise-ratio (SNR) as a result of external disturbances. is brings great difficulties to the online modal identification of concrete dams.
In this work, the RSSI method is used to extract the modal parameters of concrete dams based on the strongmotion records of the structure. e online version eigenvalue decomposition (EVD) is solved as a subspace tracking problem using the generalized yet another subspace tracker (GYAST) algorithm [22]. GYAST is a newly proposed advanced algorithm which makes several robust modifications of the original subspace tracking algorithm. e YAST algorithm relies on an interesting idea of optimally extracting the updated subspace weighting matrix in each step. In the GYAST algorithm, computation reduction of optimal subspace extraction is achieved by an approximation. It caused the GYAST has the advantages of a low computational complexity, high computational efficiency, and strong robustness. Combining GYAST with RSSI is expected to improve the accuracy and efficiency of the online modal identification. e online time-varying modal identification procedure based on the GYAST-RSSI is verified using a numerical example to evaluate the algorithm identification accuracy, robustness, and computation efficiency. Finally, the time-varying modal identification of the Pacoima arch dam is carried out using the proposed time-varying modal identification procedure and the strong-motion records of three different earthquakes to track the time-varying characteristics of the arch dam during earthquakes.

Linear Time-Varying Modal Identification
Using the GYAST-Based RSSI

Time-Varying State Space Expression of the Dam-
Reservoir-Foundation System. In the design stage of a concrete dam, the structure is generally assumed to be linear, and the observed vibration response of the structure is stationary. However, for some reasons, the damreservoir-foundation system may be time-varying. As shown in Figure 1, the dynamic constitutive model of the concrete and the base rock may be nonlinear. Additionally, the joints on the dam body may close or open with variations in the amplitude of the earthquake excitation, which will change the stiffness of the connections between different dam monoliths. e effects of the structurefoundation interaction may also introduce some timevarying factors into the dynamic features of concrete dams [23]. Moreover, a damage structure normally exhibits nonlinear dynamic behaviours and a time-dependent stiffness, and the damping variations over time lead to the time-varying modal parameters of the system. In this study, the dam-reservoir-foundation system is assumed to be a linear time-varying system. e "frozen-time" assumption that consists of modelling a linear time-varying system as a piecewise linear time-invariant system is adopted. e parameters will be assumed to vary in a stepwise way; i.e., they will change abruptly at instants separated by small time intervals and will remain constant during these time intervals. If the time intervals are sufficiently small, this model can be considered an approximation of the continuously varying case [24]. e expression of this system will be discussed in this section.
It should be noted that the material characteristics, the boundary conditions, the damping, etc. may be affected by some environmental variables, such as the water level, temperature, and humidity. e impact of these variables will not be taken into account in this study since the observation time of an earthquake response record is not very long. e equation of motion for an n degrees-of-freedom time-varying linear system S with time-varying structural parameters can be expressed as e response can be found by solving the di erential equations using initial conditions, but in general, there is no closed-form solution set to these equations. us, an approximation is introduced to solve the vibration problem. In other words, the instantaneous modal parameters are obtained at each time instant under the slowvarying assumption of the structural system. According to this assumption, the continuous "frozen" time-varying structure is made up of the mass, damping, and sti ness of every instantaneous moment τ ∈ [0, t end ], which can be expressed as follows: τ t k , k 1, 2, . . . , N, where t end is the end time of the time series of the vibration response signal, S ′ (τ) represents the time-invariant structure of t τ, and S ′ is a set of the time-freezing descriptions of the linear time-varying structure. Analytically, during the time interval [t k , t k+1 ), the parameter matrices of equation ( e state space expression for the time-varying system at discrete time intervals is as follows [26]: where A k and G k are the discrete state space matrix and observation matrix of the time interval [t k , t k+1 ), respectively. z k is a vector of state variables, and the observation vector y k is a vector of l channels and can be any or some of the three types of structural vibration responses, i.e., the displacement, velocity, and acceleration. u k is the stochastic excitation, and v k is the observational noise. e expression of the structural vibration response for a time-varying system during a time instant t ∈ [t k , t k+1 ) is discussed rst. Since the system during the interval [t k , t k+1 ) is assumed to be time invariant, a basic expression for the linear vibration system is applicable during the time interval.
To identify the complex modal shape vector, the doublesized mixing model is adopted [27]. It is well known that the shift of the phase angle between the displacement and velocity and the velocity and acceleration is 90°. e acceleration response is generally measured for the structural strong-motion response observation system. Since the physical parameters vary slowly with time, component y 0 [t] is an asymptotic signal. To track the time-varying characteristics of a structure with n degrees-of-freedom, the Hilbert transformation is performed to obtain the pseudoresponse of the structure with a 90°transferred phase angle y 90 [t].
en, in the discrete time domain, the augmented response vector y(k) ∈ R 2l×1 can be obtained as follows: Equation (4) is equivalent to the representation of the original signal in the analytical form y ′ [t] as follows: are the real and imaginary parts of the analytical signal, respectively, and j −1 √ is the imaginary unit; H · { } is the mathematical operator of the Hilbert transformation.
Based on the derivation by McNeill [28] and Masuda et al. [29], the free acceleration vibration of a structure can be expressed as Force-deflection relations for joint [25] F n Opening Closing Earthquake wave (b) tangential movement [25].
where k represents the kth time interval, i.e., and is composed of the modal acceleration responses of all the activated modes. e terms € q R [k, t] and € q I [k, t] are the real and imaginary parts of the modal acceleration response, respectively. Re · { } is the real part of a number.
If the acceleration response is measured and the output selection matrix is C a , the analytical expression of the acceleration response signal is Based on the coordinate transformation shown in equation (6), the augmented response vector y(t) shown in equation (4) is expressed as Assume the observational noise components are independently identically distributed with a variance σ 2 n . en, the covariance matrix of y[t] can be expressed as where the superscript H , denotes the conjugate transpose of a matrix and R € q € q (k, 0) is the covariance matrix of the modal acceleration responses. For a structure system with weak damping, After the discrete sampling of the physical acceleration of the l channels at N identical time intervals, a time-delayed vector Y(k) ∈ R 2pl×1 , which is composed of y(k) and its time-delayed data, can be written as follows: in which the number of time delays should satisfy 2pl > n.
en, the Hankel matrix H 0 (k) ∈ R 2pl×2pl can be defined as follows: Based on the derivation of the free acceleration response shown in equation (9), the Hankel matrix H 0 (k) can be decomposed into the following form: where . . , λ n0 are the eigenvalues of the continuous system matrix.
For the free vibration response of a structure system without damping, the modal acceleration responses of the different modes are independent of each other. However, the damping effect and the disturbance of the noise result in the matrix R € q € q (k) being approximately a diagonal matrix. Since the damping of concrete structures is usually weak, the error caused by the approximation is acceptable.
At the time instant t k , new observational data y(k) are available, and the Hankel matrix H 0 (k) should be updated online. One way to update the Hankel matrix is called the exponential forgetting method. Using this matrix, the Hankel matrix H 0 (k) at the time t � t k is expressed as in which β is the forgetting factor, with a suggested value between 0.995 and 0.999 [30].

Modal Subspace Tracking Using GYAST.
For each time step, the EVD of H 0 (k) (where H 0 (k) is a square matrix) should be performed at each time instant as follows: (15) If Q 1 (k) contains the first n dominant eigenvectors of H 0 (k), it can be taken as the dominant signal subspace, which represents the contributions of different modes. Δ 1 (k) is a diagonal matrix that is composed of the first n dominant eigenvalues of H 0 (k). Q 2 (k) and Δ 2 (k) correspond to the remaining eigenvectors and eigenvalues. However, the total computation time will then be very large and is not suitable for online applications. us, the online version of the EVD algorithm should be adopted. Since only the left part of the EVD is needed to retrieve an estimate of the observability matrix, an economical (so-called thin) EVD, which computes only that left part, can be used. Using the subspace tracking procedure, the GYAST, in each time step, a sufficient approximation of the dominant subspace, namely, can be obtained. It has been proven that, in a stationary situation, the desired steady-state weights Π(k) � Q 1 (k) H can be obtained by maximizing the mean square value of the output of the combiners, given by in which the operator trace (·) stands for the trace of a matrix.
To ensure the orthogonality of Π(k), a maximization optimization problem with constraints is expressed as follows: subject to : (18) To reduce the computational cost of updating Π(k) online, the ideal subspace projection-(SP-) type algorithm limits the search space to the range space of Π(k − 1) plus one additional search direction, i.e., Y(k). en, the n-dimensional range space of Π(k) is found as a subspace of the (n + 1)-dimensional space spanned by V(k) ∈ R 2pl×(n+1) : If Π(k) is an orthonormal basis of the (n + 1)-dimensional range space V(k), then Π(k) can be expressed as en, replacing equation (20) with the optimization problem shown in equation (18), the following new optimization problem is obtained: . For the GYAST algorithm, an approximation is introduced that calculates the modal acceleration responses en, the orthonormal basis Π(k) of the augmented (n + 1)-dimensional range space span(V(k)) is given by where Based on equation (14) and equations (22) and (23) For the subspace tracking algorithm, GYAST, to track the signal subspace of H 0 (k), four basic steps must be performed, which are (i) calculating the orthonormal basis Π(k) of the range space V(t), (ii) constructing the matrix R € q € q (k), (iii) conducting the online updating of Π(k), and (iv) conducting the online updating of R € q € q (k). Using equations (14) and (22)-(24) and the following equations (25a)-(25k), these calculations can be realized online in a recursive way, and Π(k) can be updated online for each time instant [23,31]: Shock and Vibration e ′ (k)

Flowchart of the Online Modal Identi cation Method
Using the GYAST-RSSI. e proposed online modal identi cation procedure based on the GYAST-RSSI is shown in Figure 2. Some basic steps of this procedure are summarized as follows: (1) Initialize some parameters of the procedure, such as the number of time delays p, the system order n, the forgetting factor β, and the initial data length N 1 . e system order n can be determined using the SVD spectrum method. e number of time delays should satisfy 2pl > n. e forgetting factor β is set within the range of [0.995-0.999]. rough some numerical experiments, the choice of N 1 will have an e ect on the calculation result of the initial Hankel matrix, but considering that GYAST is an iterative algorithm, the result of H(N 1 ) only a ects the results of the previous steps. Here, for convenience, set N 1 100 and k N 1 and use the initial strong-motion records of l channels with data length N 1 to calculate the initial Hankel matrix H 0 (k) as H 0 (k) (1/ p) (2) GYAST is a recursive subspace tracking algorithm, and it should rst be initialized. Perform EVD on the Hankel matrix H 0 (k) and initialize the orthonormal matrix Π(t) as follows: σ n (k) 2 smallest eigenvalue of H 0 (k).
(3) Set k k + 1 and update Π(k) using the GYAST algorithm to calculate the orthonormal basis Π(k), H(k), R € q € q (k), etc., using equation (14) and equations (22)- (24). (4) Since Π(k) ∼ Ο i,k , the two matrices have the same eigenvalues and eigenvectors, and so we set Ο i,k Π(k). From the above recursive procedures, the orthonormal matrix Π(k) can be updated. en, at the kth time interval, the discrete state space matrix A k and the observation matrix G k are realized using the following expressions: where Ο i,k O i,k (1 : 2l(p − 1), :) and Ο i,k O i,k (2l + 1 : 2lp, :). (5) After we obtain the modal identi cation results at time instant t k , some meaningless mathematical poles should rst be discarded. If the identi ed natural frequencies f i,k and damping ratios ξ i,k at time instant t k satisfy the following criteria, then mode i is deemed to be a real structural mode: Initialize some parameters, such as number of delays p, system order n, forgetting factor β, initial data with length N 1 Set k = N 1 and form the initial Hankel matrix H 0 (N 1 ) using the initial data with length N 1 Perform EVD on the Hankel matrix H 0 (N 1 ) and initialize Π(k) Obtain the estimation of matrices A and G Calculate natural frequencies, damping ratios, and modal shape vecotors k < N?

Yes
Update H 0 (k), R qq (k) R qq (k), Π (k), and 6 Shock and Vibration in which f min and f max are the minimum and maximum natural frequencies of the active modes of the structure and ξ min and ξ max are the minimum and maximum damping ratios. (6) Continue with step 3∼step 5 of the above calculation procedure until reaching the end of the strongmotion records.

Numerical Example.
Using the mass-spring-damper system with 4 degrees-of-freedom shown in Figure 3 Case 2: continuous linear change of k 1 : e initial state of the system is zero, and the support excitation a g (t) is simulated using a white noise signal. e structural vibration response is obtained by adding some white noise to the calculated responses using a Runge-Kutta algorithm (using MATLAB code ODE45). e sampling frequency of the simulated structural vibration response is 1000 Hz. Figures 4 and 5 show the support excitation and the acceleration response of m 1 for two simulation cases of k 1 .
Setting β � 0.995, p � 50, and N 1 � 100 and using the simulated structural vibration responses (SNR � 40 dB) of case 1 and case 2, the time-varying modal identification results for the frequencies and damping ratios are shown in Figures 6 and 7. e identification accuracy of the frequencies is sufficient, but the identification accuracy of the damping ratios needs to be further improved. e forgetting factor is an important parameter in the GYAST-RSSI method proposed in this paper. To verify its influence on the recognition result, different values of β are selected for a parameter sensitivity analysis. Under simulation case 1, the analysis results are shown in Figure 8. It can be seen from Figure 8 that when β is in the interval of [0.995 0.999], the result does not change much. In this example, the result of β � 0.995 is better.
Based on the simulated vibration response (SNR � 40 dB), the identification results for frequencies f 1 and f 2 of mode 1 and mode 2 using the EIV-PAST-RSSI and GYAST-RSSI methods are plotted in Figures 9  and 10 to make a comparison. Figures 9 and 10 show that the difference between the identification results of the EIV-PAST-RSSI and GYAST-RSSI methods is small, but the identification results of GYAST-RSSI are more stable.
To verify the robustness of the proposed time-varying modal identification method, simulated vibration response signals with different SNRs are used to perform the timevarying modal identification. e modal identification results of the first two frequency orders f 1 and f 2 for simulation cases 1 and 2 are shown in Figures 11 and 12, respectively. Figures 11 and 12 show that, with the increase in the noise level (decrease in the SNR), the identification accuracy of the frequencies is still very high. e modal assurance criterion (MAC) indexes between the identified modal shape vectors and the theoretical modal shape vectors are calculated and shown in Figure 13. In addition to some sudden changes in the MAC value at 4 s, 8 s, and 12 s, which is caused by the instability of the algorithm, the identification results of the modal shape vectors indicate the sufficient accuracy of the GYAST-RSSI method of tracking the time-varying structural modal shape vectors.
To evaluate the identification accuracy of the proposed time-varying modal identification method, the average relative error (ARE) of the frequencies is calculated and is shown in Table 1. e ARE of the certain order frequency is defined as (31) in which f t m and f c m are the identification results and the theoretical value of the certain order frequency at time instant m, respectively.
From Table 1, it can be seen that the identification accuracy is good for the two time-varying modal identification methods since the maximum ARE is 2.11%. With the decrease in the SNR, the identification accuracy for the frequencies has become low. Moreover, the GYAST-RSSI has better identification than the EIV-PAST-RSSI method in most cases.

Time-Varying Modal Identification of the Pacoima Arch Dam.
e Pacoima dam is a 113 m high arch dam located near Los Angeles, California. Since the Pacoima arch dam is near the famous San Fernando earthquake-prone area, it has experienced several earthquakes of different magnitudes. e dam experienced two severe earthquakes: the San Fernando earthquake (1971) and the Northridge earthquake (1994). Due to the damage caused by the San Fernando earthquake and the opening of the joints near the thrust block at the left abutment, rehabilitation operations were Shock and Vibration 7 support excitation   performed afterward, and an array of 17 accelerometers were installed on the dam in 1977 [32]. e arrangement of the 10 acceleration sensors with 17 measurement channels is shown in Figure 14. Among the 17 measurement channels, channels 1∼8 are located on the dam body and channels 9∼17 are located on the base rock. e measurement directions of channels 1, 2, 5, 6,7,8,9,12, and 15 are in the radial direction; channels 4, 11, 14, and 17 are in the tangential direction; and channels 3, 10, 13, and 16 are in the vertical direction. After the Northridge earthquake, the strong-motion observation system was updated, and the strong-motion responses of the three di erent earthquakes were recorded. e information of the three earthquakes is shown in Table 2.
e strong-motion records of the Pacoima arch dam for channels 1, 2, 5, 6, 7, and 8 in the radial direction during the Newhall earthquake are shown in Figure 15. e sampling frequency of the strong-motion monitoring data is 200 Hz.
us, the cuto frequency of the data is 50 Hz. Time-invariant modal identi cation results of the Pacoima dam using di erent methods are shown in Table 3.
ese identi cation results are used to make a comparison with the following time-varying modal identication results.   Channels 9∼17 are located on the base rock. e strongmotion records of channel 9 in the radial direction during the three major California earthquakes and the corresponding power spectrums are shown in Figure 15. As seen from Figure 15, for the San Fernando earthquake and the Chino Hills earthquake, the dominant frequencies of the earthquake excitations are less than 3.32 Hz. It can be seen from Table 3 that the dominant frequencies of the earthquake excitation are     less than the rst-order modal frequency of the dam. However, for the Newhall earthquake, the range of the dominant frequencies of earthquake excitation is much wider than those of the other two earthquakes. erefore, for the Newhall earthquake, the nal modal identi cation results of the dam should not include these dominant frequencies of excitation shown in Figure 16(f).
By setting β 0.995, p 100, and N 1 200, the structural frequencies and damping ratios are identi ed using the proposed time-varying modal identi cation method based on GYAST-RSSI, and the identi cation results are shown in Figures 17-19. e vibration response of channel 1 is plotted to show the variation of the modal parameters with the acceleration amplitude. From Figures 17-19, the timevarying characteristics of the arch dam during earthquakes are clearly observed, especially when the amplitude of the acceleration response is relatively large. According to Figures 17-19, the frequency after the earthquake tends to be  stable, indicating that none of the three earthquakes may cause any damage to the dam.
For the three earthquakes, the frequencies of the dam show a decreasing trend after the earthquake, which may be caused by the solid-water interaction and the reduction of the structural sti ness. As was shown in some previous works, the solid-water interaction can be simulated by the simple adding-mass method. e additional mass that is added to the structural system will reduce the frequencies of the damreservoir-foundation system. However, the decrease in the frequencies can also be caused by the reduction of the structural sti ness during earthquakes, which needs to be     Figure 20. Figure 20 shows the time-varying characteristics of the arch dam's modal shape vectors during the earthquake.

Conclusions
e online modal identi cation method of concrete dams proposed in this work has been veri ed using a numerical example and an engineering example. e GYAST is a good   19.92Hz (f ) Figure 16: e strong-motion records of channel 9 and the corresponding power spectrums during three major California earthquakes: (a) the strong-motion record of channel 9 during the San Fernando earthquake; (b) the power spectrum of the strong-motion record of channel 9 during the San Fernando earthquake; (c) the strong-motion record of channel 9 during the Chino hill earthquake; (d) the power spectrum of the strong-motion record of channel 9 during the Chino hill earthquake; (e) the strong-motion record of channel 9 during the Newhall earthquake; (f ) the power spectrum of the strong-motion record of channel 9 during the Newhall earthquake. tool to track the signal subspaces of different modes to reduce the disturbance of the observation noise and stochastic excitation since it shows good identification accuracy, robustness, and computation efficiency. e identification accuracy of the frequencies and modal shape vectors is sufficient, while for the damping ratios, the performance of the proposed online modal identification method still needs to be improved. Some possible ways to obtain a better identification result of the damping include improving the robustness of the subspace tracking algorithms, signal denoising, and adopting input-output modal identification methods. e time-varying modal parameters identified using the strong-motion records of concrete dams can be used as important features to study the real dynamic characteristics of structures during large earthquakes, which is our research direction in the future.
Data Availability e strong earthquake records data used to support the findings of this study may be released upon application to the Center for Engineering Strong Motion Data at https:// strongmotioncenter.org. e MATLAB code for the GYAST algorithm online to support the findings of this study is available from the Professor Mostafa Arjomandi-Lari [22]. e numerical simulation and other data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.