Modal identification of building structures under unknown input conditions using extended Kalman filter and long-short term memory

2.1Fundamental algorithm of extended Kalman filter

The Kalman filter is a recursive filter algorithm developed by [59]. It performs estimation by updating the state vector while correcting the residual between the sensor response and predicted response. Hoshiya and Saito [28] attempted to estimate system parameters using the EKF algorithm by expressing the dynamic governing equation as a state space equation. Until recently, the EKF algorithm has been actively used in SI, and this section briefly introduces the EFK algorithm. The state space equation of a nonlinear dynamic system to which an earthquake load is applied can be expressed as shown in Eq. (1):

In Eq. (1), 𝐗⁢(t) denotes the state vector for continuous time and 𝐮¨g⁢(t) denotes the ground acceleration applied to the dynamic system. 𝐰⁢(t) denotes the process noise vector with a covariance matrix 𝐐⁢(t). The measurement equation can be expressed as follows:

In Eq. (2), 𝐘k denotes the output vector at t=k⁢Δ⁢t for discrete time and 𝐗k denotes the state vector at t=k⁢Δ⁢t. 𝐯k denotes the measurement noise vector with a measurement covariance matrix 𝐑k. The state space equation for continuous time in Eq. (1) must be expressed in discrete time for interaction with the measurement equation. When t=k⁢Δ⁢t, the state vector for discrete time can be obtained for every section (k+1)⁢Δ⁢t through integration. Therefore, the predicted state space equation can be evaluated as shown in Eq. (3):

where the superscript of the state vector 𝐗^ indicates the predicted state. In addition, the subscript k+1/k indicates an uncorrected state vector. The predicted error covariance matrix used in Kalman gain for state correction can be calculated by applying Eq. (4):

Here, 𝐅k+1/k denotes a state transient matrix, and in order to update the state transient matrix, a Jacobian matrix is calculated using Eq. (5):

The predicted measurement vector can be obtained by utilizing Eq. (6), which represents the state vector corresponding to the measurement sensor response in Eq. (3), which is the uncorrected predicted state:

In order to update the predicted state vector of Eq. (3), the Kalman gain is calculated as shown in Eq. (7). 𝐇k+1 in Eq. (7) is a Jacobian matrix that is a partial derivative of the measurement vector 𝐡⁢(𝐗k+1) with respect to the state vector 𝐗k, which is given by Eq. (10). The predicted state vector in Eq. (3) is updated as shown in Eq. (8) by correcting the residual of the sensor measurement data 𝐘k+1 and the predicted measurement vector 𝐘^k+1 with Kalman gain, and the updated error covariance matrix can be calculated by applying Eq. (9):

2.2Application of the modal extended dynamic system to the extended Kalman filter

The dynamic governing equation of the multi-degrees of freedom system for earthquake load is as shown in Eq, (11).

M is the mass matrix, C the damping matrix, and K the stiffness matrix. In addition, 𝐮¨, 𝐮˙, u, and 𝐮¨𝐠 indicated the acceleration, velocity, displacement, and ground acceleration of the system, respectively. To estimate the modal parameters, Eq. (11) is expanded into modal coordinates.

The eigenvector (mode shape) v acquired through eigen value analysis by the determinant of Eq. (12) can be represented as the modal matrix 𝚽 for all modes. After multiplying this to both side of Eq. (11) and superpositioning the mode, the modal equation for the jth mode can be presented as shown in Eq. (13).

Here, ω is the natural angular frequency and ζ the modal damping ratio. The base j denotes the jth mode. Lj is the modal participation factor, and q¨, q˙, and q represent the modal response of the acceleration, velocity, and displacement, respectively.

Applying the modal equation of Eq. (13) that was expanded to modal coordinates to the EKF algorithm requires the sensor measurement vector of Eq. (8) 𝐘k+1 for updating the state vector in EKF must be obtained as a modal response. However, acquiring the modal response without additional filtering process is difficult. Therefore, the predicted measurement vector of Eq. (6) 𝐘^k+1 must be acquirable with the system response for application in reality. Using Eq. (14), the modal response can be represented as the system response.

In Eq. (14), ϕi⁢j indicates the jth mode eigenvector element from the ith degree of freedom. Additionally, n is the total number of MDOF system’s degree of freedom. Eq. (14) implies that the physical behavior of the MDOF system at the ith degree of freedom is equal to the results of the total sum of modal responses of j= 1, …, n. In Eq. (14), ϕi⁢j multiplied by q¨j separates by each degree of freedom the modal response q¨j that includes the response of all degrees of freedom. Thus, the modal response separated for each degree of freedom is represented by Eq. (15).

Based on Eq. (15), the following Eq. (16) is obtained by multiplying ϕi⁢j to both sides of Eq. (13).

In Eq. (16), Γ is the modal participation factor separated for each degree of freedom, and Eq. (16) is a modal equation separated for each degree of freedom from the modal equation of Eq. (13).

Therefore, to convert the predicted measurement vector of Eq. (6) into the form of system response for each degree of freedom, Eq. (16) was developed into an equation for the acceleration, which was reproduced into Eq. (17) as an expression of the sum of the acceleration.

The number of modes in Eq. (17) depends on the total number of the degree of freedom in the direction of the sensor’s measurement. Therefore, n can be up to as many as the degree of freedom. The predicted measurement vector 𝐘^k+1 of Eq. (6) can be determined by the following Eq. (18) through Eq. (17).

The predicted measurement vector 𝐘^k+1 of Eq. (16) is determined by the sum of modal response. Therefore, 𝐘^k+1 can be obtained by the system response instead of the modal response. This research uses the top floor response as the system response of Eq. (8). The reason for this is because in general the modal responses of the top floor responses have higher the modal participation ratio in system response with a mode of a lower degree. The following Eq. (19) is the state vector of the structure with three degrees of freedom. In this case, the mode occurs up to the third mode, for which n=3 determined for Eq. (18) and the predicted measurement vector 𝐘^k+1 is determined as ∑i=13u¨i⁢j if the acceleration response is obtained at the top floor. Additionally, to estimate the state variables of ω and ζ, the composition of the state vector that includes ω and ζ is as shown in Eq. (19).

The state space equation for Eq. (19) is as shown in Eq. (20).

2.3LSTM network for input load estimation

The EKF-LSTM approach developed in this study is a method of estimating the input load by constructing a network between the structural response and input load as an LSTM algorithm. Therefore, the input load becomes the output data in the LSTM algorithm and the structural response becomes the input data. Figure 1 shows a schematic diagram of the specific EKF-LSTM method developed in this study. In Fig. 1, the LSTM unit can be determined by the time series length or the number of time steps. In addition, in the LSTM unit, information about the time sequence of training data is transmitted through the cell state and hidden state. Regarding the learning process in Fig. 1, structural responses as long as the time series are input into the LSTM unit as input data. After passing through the learning layers of the LSTM layers, fully connected layer, and regression layer, the predicted output and output of the input load, which constitute the output data, are compared with a loss function, and the weight of the LSTM unit is updated using the optimizer. As learning by batch progresses, learning state information of time series data between batches is delivered by stateful. Therefore, the time series of the training data is learned dependently regardless of the batch. Further details on the LSTM algorithm can be found in the literature [60, 61]. Furthermore, the loss function used in this research as the mean square error is determined by Eq. (21).

In Eq. (21), S is the sequence length of the train data and R the dimension of the output data.

In addition, the adaptive moment estimation (Adam) optimizer applied for optimizing the weight value was proposed by [62]. The weight value update of the train network is determined by the following Eq. (22).

θ is the weight, bias parameter, η the learning rate, m^ the estimate of the first moment (the mean) of the gradients, and v^ the second moment (the uncentered variance) of the gradients. ϵ as the default value is proposed as 10-8[62, 63].

Figure 1.

Schematic of LSTM for the estimation process of input loads.

2.4EKF-LSTM

The general flowchart for EKF-LSTM is as shown in Fig. 2. To compose the LSTM network that can estimate the input load, the structural response and the input load, which is the ground acceleration, are acquired from the SI target structure. The network using the obtained structural response and input load as training data is trained through LSTM as shown in Fig. 1. Training proceeds only for the initially set max epoch. The initial state vector is established to use the EKF method for SI. The input load that was estimated through the LSTM network trained during the calculation of the predicted state space equation is used. This study is composed of a state vector that includes the modal responses, frequency, and damping ratio for estimating the modal parameters as shown in Eq. (19). In addition, as mentioned in Section 2.2, the sensor measurement data 𝐘k+1 for application to Eq. (8) uses only the top floor responses for proceeding with the EKF.

Figure 2.

Flowchart for EKF-LSTM.

Figure 3.

Training data obtained in numerical verification model: (a) System parameters of 3-DOF model; (b) Structural response used for the input; (c) Ground acceleration used for the output.

3.13-DOF numerical model

To verify the proposed method, a dynamic system model with three 3-DOF was constructed and structural responses were obtained through Newmark numerical analysis. The specific system parameters for the model included a mass of 5 kg for each of the first, second, and third floors, with each floor having a stiffness of 1000 N/m. When the model was subjected to eigen value decomposition (EVD), the natural angular frequency was approximately 6.29 rad/s, 17.63 rad/s, and 25.48 rad/s for the first, second, and third modes, respectively, and the damping ratio for each mode was assumed to be 1%. The applied load was a white-noise input load. Furthermore, the sampling frequency is 100 Hz. Figure 3 shows the training data used as input and output data. To divide the training data shown in Fig. 3 into the input and output datasets, 200 pieces of data were generated for each dataset and used for LSTM learning. The learning rate was 0.001, minibatch size was 20, maximum epoch was 5,000, and maximum number of iterations was 50,000. The gradual decrease of the loss function in Fig. 4 indicates that the LSTM network was trained normally.

Figure 4.

Loss function of long short-term memory (LSTM).

Figure 5.

Estimated ground acceleration from LSTM: (a) Comparison between the estimated and the exact ground accelerations, and (b) Scatter plot of predicted and exact results.

Figure 5a compares the estimated and exact ground acceleration results. As shown on the time axis of Fig. 3b, the training was conducted with a time series length of 100 s, but the LSTM algorithm showed results over 10 s. As such, the LSTM algorithm enabled estimation regardless of time. To compare the similarity between the exact ground acceleration and input load estimated by LSTM, the RMS results were compared. The RMS results can be used as a good index for representing the average for the signal that intersects the 0-axis. The root-mean-square (RMS) result for the exact ground acceleration was 9.8739e-05 m/s2, and the estimated result was 1.0336e-04 m/s2. It can be observed that the relative error of the two for RMS is estimated to be approximately 4.6%. Figure 5b shows a scatter plot between the results estimated by LSTM and the observed (measured) results. The more the predicted results in the scatter plot match the observed results, the more linear the plot. In addition, if an estimated result in Fig. 5a shows a shifted bias error off the axis of 0, the scatter plot result will appear as a parallel translation result. In Fig 5b, the results estimated by using the LSTM algorithm show a similar gradient to the observed results, with no bias error. Furthermore, the R-squared value, representing the scale of regression between the predicted and observed results, was 0.9201, indicating a similarity of approximately 92%. Therefore, the input load estimation using the LSTM network proposed in this study proceeded normally, suggesting that EKF-based SI was possible under unknown input using LSTM.

Figure 6.

Natural frequency estimation and relative error plots. (a) Natural frequency and (b) Modal damping ratio.

As it is necessary to set the initial values of the state variables of the initial state vector to execute EKF-based SI, the initial values for modal responses up to the third mode were set to zero, the natural angular frequencies and the modal damping ratios were set to ω1= 10 rad/s, ω2= 25 rad/s, ω3= 30 rad/s and ζ1=ζ2=ζ3= 0.5, respectively. The initial error covariance matrix P0 was set to diag{1, 1, 1010, 1010, 1010, 1, 1, 1010, 1010, 1010, 1, 1, 1010, 1010, 1010}, process noise covariance matrix Q to 10-4I15, and measurement noise covariance matrix R to 10-4I1. The input load was estimated using LSTM for 100 s, equal to the time series length of the training data, and was applied to EKF for SI. Figure 6 presents the SI results obtained by using the EKF-LSTM method. Figure 6a shows only the convergence results up to 30 s, as the obtained natural angular frequency is a constant value after 30 s. The RMS natural angular frequency after 20 s estimated by the EKF was identified to be 6.29 rad/s, 17.64 rad/s, and 25.60 rad/s for the first, second, and third modes, respectively. The EVD results were approximately 6.29 rad/s, 17.63 rad/s, and 25.48 rad/s, with an identification error of less than 1%. The modal damping ratio was identified as 0.0101, 0.0099, and 0.0098 for the first, second, and third modes, respectively, after 20 s. The damping ratio was identified with a maximum relative error of 2%. Considering the error at the fourth decimal place in the damping result, even the damping result, which showed a relatively large error, also provided similar estimation results, proving the effectiveness of the EKF-LSTM method developed in this study. Regarding the convergence speed of the state vector, the convergence speed to the correct value of the natural angular frequency and modal damping ratio can be determined by the P0 matrix. If the P0 matrix values are greater than the determined values in this research, the convergence speed could be improved, but the state vector stability will decrease. Moreover, the Q and R matrices are initial parameters determined by the process noise and measurement noise, and they are one of the critical parameters along with the P0 matrix that determine the state vector convergence. Accordingly, the P0, Q, and R matrices were determined through the customary trial and error method. The individual optimization method of the initially set EKF parameters can be found in the following reference [41].

Figure 7.

System response estimation: (a) 1st modal response of displacement, (b) 2nd modal response of displacement, (c) 3rd modal response of displacement.

Table 1 presents the identification results using the EKF-LSTM and conventional EKF methods. To examine the impact of noise on the EKF-LSTM method, results for the case of the signal noise ratio (SNR) or the input and structural response being 20 dB was included. The LSTM network for estimating the input load when SNR = 20 dB was additionally composed, and the training data having SNR = 20 dB and no noise was applied. The learning data and max epoch for training was identical to the value set above. Additionally, the error listed on the table was relative error, which was indicative for the results of identification without noise. The results of EKF-LSTM and conventional EKF methods were obtained within the relative error of 1%. No significant difference was presented between the two methods, other than the EKF method using accurate values obtained from numerical analysis of the input load but the EKF-LSTM method using input load estimated by the LSTM network not generating noticeable error from the estimated results. Thereby, the accuracy of the input load estimation by the LSTM network was demonstrated, which accentuates the advantages of the use of the EKF-LSTM method for unknown inputs. Additionally, the application of the LSTM method eliminates computational time issues that may occur during input load estimation. In the case of SNR = 20 dB noise being included, the error rate of convergence results tended to increase in all primary to tertiary modes. In particular, the error of the third mode increased significantly, which was attributed to the noise and high uncertainty of the third mode response because the top floor response was used. Further, the error of the EKF-LSTM and conventional EKF methods increased. In this study, noise was not included in the previously trained LSTM network. Therefore, the LSTM network that was trained without noise did not generate good results for responses including noise. However, the LSTM network that uses responses with and without noise as training data was able to estimate the input load up to regression scale of R square 93% both when SNR = 20 dB and noise was absent. Further studies on the impact of noise are required in detail for the application of this method in real buildings because noise problems can occur due to environmental impact. Specific analysis regarding this issue is discussed in subsequent research on real building application.

Furthermore, the EKF-LSTM method can estimate the modal response because it is composed of state variables, as indicated in Eq. (19). Figure 7 presents the modal displacement responses of the first through third modes estimated on the third floor. Figure 7a shows the first modal displacement response, which is already similar to the system displacement response with only the first mode. This indicates that the contribution of the first modal response to the system response was large. Because the sum of the first through third modal displacement responses represents the system displacement response, the RMS result of adding all the modal responses obtained from the EKF-LSTM method was estimated to be 3.3149e-06 m. The RMS result of the system displacement response was 3.2834e-06 m, and the relative error with the estimated result was identified as being less than 0.96%. This result was smaller than the 4.6% error of the input load estimated by using the LSTM algorithm. Because the RMS error of the input load was approximately 4.6%, the estimated response could also have an error of approximately 4%, but the correction effect obtained through the Kalman gain of the EKF algorithm seemed to have led to a more accurate estimation. Thus, the EKF-LSTM method proposed in this study was also effective in system response estimation.

Figure 8.

Steel-frame structure model design: (a) Front view, (b) Side view, (c) Experiment model.

4.1Experiment verification using steel-frame structure

In this section, the structural response was obtained by performing excitation experiments on the steel frame structure model, and the effectiveness of the EKF-LSTM method was verified using this response. The design information about the steel frame structure used in the study is provided in Fig. 8. The accelerometer used was a PCB333B50 sensor from PCB, which had a sensitivity of 1000 mV/g and a full-scale range of ± 5 g. The data logger was an IOtech 640u from NI Instruments. The shaker was an APS400 from APS Dynamics capable of long strokes and low-frequency excitations. To verify the EKF-LSTM method, measurement data were acquired at a sampling frequency of 256 Hz, and the measurements were performed for approximately 300 s.

The responses for approximately 235 s were used as the training data by selecting only the responses with white noise from the measurement data for the initial 300 s. Figure 9 shows the acquired training data.

The obtained response was subjected to a fast Fourier transform (FFT) to obtain the natural frequency and mode vector of the steel frame structure model. When the FFT was performed, the natural frequencies were obtained as 4.48 Hz, 14.08 Hz, and 20.09 Hz for the first, second, and third modes, corresponding to 28.148 rad/s, 88.467 rad/s, and 126.229 rad/s, respectively. Because the mode vector result corresponding to the natural frequency showed the first through third eigen shapes in Fig. 10b, the measurements seemed to have been performed accurately.

Figure 9.

Training data: (a) Acceleration for the input and (b) ground acceleration for the output.

Figure 10.

Modal parameters of steel frame model.

The modal assurance criterion (MAC) was calculated based on the acquired mode vector to apply the FDD method, and the modal damping ratio was identified using a spectral bell function of MAC 0.9 or higher. The modal damping ratios obtained through the FDD method were 0.0113, 0.0055, and 0.0033 for the first, second, and third modes, respectively. To use the training data in Fig. 9 for learning, the 60,000 data points for the total time series were divided into groups of 400 data points, yielding 150 training data sets. Among them, 30 data sets (20%) were used as test data. As for the hyperparameter settings for learning, the learning rate was set to 0.0001, batch size to 60, and maximum epoch to 200,000. Figure 11 shows the loss function of the trained LSTM network. As the learning progressed 400,000 times, the loss function curve gradually converged, suggesting that the learning had progressed normally.

There were 30 test datasets in total, and considering that each dataset had 400 data points and the sampling frequency was 256 Hz, each dataset corresponded to approximately 1.56 s. Therefore, the test data had a time series length of approximately 46 s, and Fig. 12a shows the ground acceleration predicted by the LSTM network for 46 s. In the R-squared result representing the regression scale of the scatter plot in Fig. 12b, approximately 93% of the curve fitting is shown.

Figure 11.

Loss function plot as a function of the number of iterations.

Figure 12.

Estimated ground acceleration based on LSTM: (a) Acceleration plot as a function of time, and (b) Scatter plot.

SI was performed using the predicted ground acceleration in Fig. 12 for the EKF. To execute the EKF, the initial modal responses up to the third mode were set to zero, the natural angular frequencies and the modal damping ratios were set to ω1= 40 rad/s, ω2= 100 rad/s, ω3= 150 rad/s, and ζ1=ζ2=ζ3= 0.5. The initial error covariance matrix P was set to diag{1, 1, 106, 106, 106, 1, 1, 106, 106, 106, 1, 1, 106, 106, 106}, process noise covariance matrix Q to 10-4I4, and measurement noise covariance matrix R to 10-1I1.

Figure 13.

Modal parameters of steel frame structure (a) Estimated natural frequencies, and (b) Modal damping ratio.

Figure 14.

Estimated responses: (a) system and estimated acceleration responses, and corresponding responses for the (b) 1st, (c) 2nd, (d) 3rd modal accelerations

Figure 13 shows the modal parameters estimated by the EKF-LSTM method. As shown in Fig. 13a, the natural angular frequency was calculated to be 28.148 rad/s, 88.467 rad/s, and 126.229 rad/s for the first, second, and third modes, respectively. The predicted first natural angular frequency was 28.36 rad/s when the RMS was averaged after 10 s, showing a relative error of 0.75% from the results obtained by FFT. In addition, the predicted second natural angular frequency was 88.71 rad/s, indicating a relative error of approximately 0.2%. The predicted third natural angular frequency showed a slower convergence rate than the first and second modes, which may have been due to the low contribution of the third mode to the top floor response, as described in Section 2.2. The predicted third natural angular frequency was 126.882 rad/s from the RMS average at approximately 30 s, showing a relative error of approximately 0.5%. Figure 13b depicts the modal damping ratios, and the RMS results were identified as being 0.0105 and 0.0051 for the first and second modes, respectively, after approximately 10 s. These results indicated a difference of approximately 7% between the results identified by the FDD method and the maximum relative error. Considering the uncertainty of the modal damping ratio, the EKF-LSTM method seemed to show effective performance in the identification of the modal damping ratio. The third mode converged more slowly than the other modes, showing an unstable convergence curve, similar to the natural angular frequency result. This instability may have been due to the fact that the third mode had a low contribution to the system response. The RMS result from approximately 25 s after convergence was identified to be 0.0028.

Figure 14 shows the modal acceleration response estimated by the EKF. In the modal responses in Fig. 14b–d, the mode contribution to the system response changes with increasing mode order. In Fig. 14a, the response by adding all the first to third modal responses appears in a shape similar to the system acceleration response. The relative error between the average RMS result in Fig. 14a, which is the sum of the estimated modal responses, and the system response was calculated to be less than 1%. Table 2 summarizes the modal parameter estimation results by comparing the EKF-LSTM and FDD methods. To explain further, the modal response should show the shape of a sinusoidal wave, where the response scale becomes smaller for a higher mode. This tendency was not observed in Fig. 14b–d as the data were shuffled when composing the training data.

Figure 15.

Computational time for each method.

Figure 15 compares the computational time according to data length of the EKF-LSTM and FDD methods. In the application of EKF-LSTM, data was accumulated and the change in computational time had a parabolic form. The FDD method showed a linear pattern in the change in computational time, and this method could be considered to be more efficient. However, the results of the modal damping ratio of the FDD method shows a continuous decrease which approximated to about 0.01 from 200 s and above. Therefore, data for 200 s or above must be used for the application of the FDD method. In the application of the EKF-LSTM method, results of about 0.01 is obtained relatively constantly from about 10 s data. Thus, the computation of EKF-LSTM method was concluded to be more efficient. This method based on the EKF algorithm is able to estimate data real-time for each discrete time delta t. Meanwhile, the FDD method uses all measurement data at once for modal analysis, and identification is impossible until the analysis is complete.

To summarize the results obtained for the steel frame structure model using the EKF-LSTM method, the estimated natural angular frequency showed an error of less than 1%; the modal damping ratio, which has a relatively large error, shows an error within about 15%. Furthermore, the sum of the modal response showed an estimation result of less than 1% error compared with the displacement response. Among the results of the modal parameters, the modal damping ratio showed a rather large error of 15%, but it was still judged to be reliable considering the uncertainty of the modal damping ratio. Therefore, the effectiveness of the EKF-LSTM method proposed in this study is demonstrated in actual structural applications.