Modal parameters identification of bridge by improved stochastic subspace identification method with Grubbs criterion

In the wind tunnel test of a long-span bridge model, to ensure that the dynamic characteristics of the model can satisfy the test design requirements, it is particularly important to accurately identify the modal parameters of the model. First, the stochastic subspace identification algorithm was used to analyze the modal parameters of the model in the wind tunnel test; then, Grubbs criterion was introduced to effectively eliminate outliers in the damping ratio matrix. Stochastic subspace identification algorithm with Grubbs criterion improved the accuracy of the modal parameter identification and the ability to determine system matrix order and prevented the modal omissions caused by determining the stable condition of the damping ratio in the stability diagram. Finally, Oujiang Bridge was used as an example to verify the stochastic subspace identification algorithm with Grubbs criterion and compare with the results of the finite element method. The example shows that the improved method can be effectively applied to the modal parameter identification of bridges.


Introduction
Modal parameters are key parameters for the design and establishment of aeroelastic models. Accurate modal parameters can establish models that satisfy the test requirements, and improve the efficiency of the wind-resistant design of long-span bridges. At present, the commonly used methods of modal parameter identification in engineering usually include experimental methods and working methods. 1,2 The experimental method is a traditional method of modal parameter identification. The method uses artificial excitation to make the model produce free vibration. Based on the displacement response data of the measured points, the specific modes of the bridge can be identified. Based on the hammer excitation and varied-time-base transfer function analysis, the modal parameters of the aeroelastic model of the Runyang Yangtze River Bridge were identified by Zhenshan et al. 3 Based on the initial displacement pulse excitation, the modal parameters of the aeroelastic model of the Sutong Bridge were identified by Fuyou et al. 4 The working method directly uses environmental excitations (such as vehicles, pedestrians, and wind load) acting on the structure as the input of the system. There is no need to artificially impose incentives on the structure. The modal parameter identification can be completed based on output data of the system. 5,6 Based on the random vibration frequency domain method, the modal parameters of the aeroelastic model of the third bridge in the Nanjing Yangtze River were identified by Ledong et al. 7 He Neng used the time domain method of the stochastic subspace algorithm based on environmental excitations to identify the modal parameters of the aeroelastic model of Maputo Bridge.
Neng H verified this method through vibration tests under environmental excitation. 8 Based on the fast Fourier transform (FFT) method, the modal parameters of the rigid-continuous combination beam bridge under environmental excitations were identified by Jianping and Peijuan. 9,10 In addition, based on working methods, scholars have used frequency domain decomposition (FDD) and stochastic subspace identification (SSI) to identify modal parameters. 11,12 For large scale civil structures such as long span bridges, the output-only SSI is very effective to identify and monitor these structures due to ambient vibrations. There are various SSI techniques such as covariance-driven SSI, data-driven (DATA-SSI), and a combination of other methods, such as the expectation maximization technique (EM-SSI). 13,14 The main output-only identification in the DATA-SSI techniques is the orthogonal projections performed by LQ decomposition. 15 Then, it is followed by singular value decomposition (SVD) to extract the subspace system. In DATA-SSI, there are variants that correspond to different choices to weight the matrices before factorizing projection matrix. Weng et al. 16 studied the DATA-SSI to investigate the dynamic characteristics.
In wind tunnel tests, in addition to accurately recognizing the vibration frequency of the structure, it is necessary to accurately identify the damping ratio of the structure. Due to noise during measurements, the singular value changes, resulting in changes of system matrix and ultimately the modal parameters. 17,18 Therefore, it is easy to have a small difference in frequency and mode and a large difference in damping ratio and even modal omissions. In this study, stochastic subspace identification algorithm is used to identify the modal parameters. The most important thing in stochastic subspace identification algorithm is to determine system matrix order. The system matrix order is an empirical parameter, and there is currently no unified method to determine it, which needs to be determined in advance. Three methods are commonly used to determine the system matrix order. Firstly, decompose singular value on projection matrix and get the number of singular values that are not zero. Secondly, use singular values jump method. Thirdly, use stability diagram method. Stability diagram method is used to determine the system matrix order in this paper. However, when stability diagram method is applied to stochastic subspace identification algorithm, the modal omission might be caused due to larger damping ratio dispersion. Therefore, stochastic subspace identification algorithm with Grubbs criterion is introduced to optimize the process of determining the damping ratio stability. In Grubbs criterion, the critical value to judge outliers has nothing to do with the average and standard deviation of the sample, so it is easily adjusted. The advanced algorithm can effectively avoid modal omissions and improve the accuracy of damping ratio identification and the precision of stability diagram method through numerical simulation result demonstration of Oujiang Bridge.

Stochastic subspace identification algorithm
Space model in the discrete state of the system Figure 1 shows a k-degree-of-freedom viscous damping system. Its vibration differential equation is defined as follows: where M 2 R k3k , C 2 R k3k , K 2 R k3k are the mass matrix, damping matrix, and stiffness matrix of the system, respectively; q(t) 2 R k is the displacement vector of the system at time instant t; F(t) 2 R k is the external excitation vector at time instant t; B 2 2 R k3m is the system input position matrix; u(t) 2 R m is the input vector of the system at time instant t; m is the number of excitation points. According to the system control theory, the state space model includes the state equation and output equation. Where the state equation is a first-order differential equation, and equation (2.1) is a second-order differential equation, which obviously does not meet the requirements. When where x(t) 2 R n is the state vector; A c 2 R n3n is the continuous state matrix, which represents the internal state of the system; B c 2 R n3m is the input matrix of continuous time, which represents the relationship between the input and the system, and n = 2k. Because the dynamic response (displacement, velocity, acceleration, etc.) of the system can be measured by the corresponding sensor, its response data y t ð Þ can be expressed as follows: where C a 2 R l3k , C v 2 R l3k , and C d 2 R l3k are the output matrices of acceleration, velocity, and displacement, respectively; y(t) 2 R l is the output vector; l is the number of monitoring points. With the introduction of x t ð Þ, equation (2.3) can be expressed as follows: where is the output matrix of continuous time, which represents the relationship between system output and system state; D c 2 R l3m is the transfer coefficient matrix, which represents the relationship between system output and system state.
Through the conjunction between equation (2.2) and equation (2.4), the state space model of continuous time of the system can be obtained as follows: However, equation (2.5) is not consistent with the actual situation and must be discretized because the data collected by the actual sensor are discrete in time instead of being continuous. Suppose that the initial condition of the system is x t ð Þ = x t 0 ð Þ when t = t 0 . Then, the general solution of x(t) can be expressed as a combination of free decay under initial conditions and forced vibration under random loads.
Researchers perform discrete sampling of continuous time series, where the sampling period is Dt. Supposing that t 0 = kDt and t = k + 1 ð ÞDt, which can be substituted into equation (2.6), we can obtain Assume that the input vector u t ð Þ of the system is a constant in a sampling period. t 0 = k + 1 ð ÞDt À t. The state vector x kDt ð Þ at time kDt is simplified to x k . Then, equation (2.7) can be rewritten as follows: Likewise, equation (2.4) can be rewritten as follows: where C = C c and D = D c . Since there must be noise interference in the input and output of actual engineering, equations (2.8) and (2.9) show that the discrete state space model of the system considering noise interference is as follows: where w 0 k and v 0 k are the input noise and output noise, respectively. Suppose that environmental excitation is an ideal white noise excitation. The excitation term in equation (2.10) is combined with the input and output noise terms, since the white noise excitation is a random excitation, cannot be measured and is similar to the input and output noise. The random state space model is obtained as follows: where A is the discrete state matrix; C is the discrete output matrix; E is the mathematical expectation factor; d pq is Kronecker Delta function; Q, S, and R are the covariance matrices of w k and v k .

Modal parameter identification
After obtaining system matrix A and output matrix C, the modal parameters of the system can be obtained by solving the eigenvalues of the system matrix as follows: where C is the eigenvector matrix of the system, whose value is a complex number, L is the diagonal matrix composed of system eigenvalue l i , and the eigenvalues are a pairwise conjugate complex number.
In the process of deriving the discrete state space equation of the system, considering the sampling limitation of the actual instrument, the continuous-time state space equation of the system must be discretized. However, the actual system is continuous, so the authors transferred the discrete system state matrix A into a continuous system state matrix A c when solving the modal parameters. The relationship between A and A c can be rewritten as follows: The relationship between the two eigenvalues can be obtained as follows: where, j is the imaginary unit. The relationship among the eigenvalue of the system, system circular frequency v and damping ratio j can be expressed as follows: According to equation (2.16), the modal frequency can be expressed as follows: the damping ratio can be expressed as follows: and the modal shape can be expressed as follows:

Improved stochastic subspace identification algorithm
In the calculation or theoretical analysis, the assumed damping form may not be completely consistent with the actual situation, and the singular value changes due to the interference of the measurement noise, which causes a change in the system state matrix and finally a change in modal parameters. The damping ratio is more sensitive to this change, so there tend to be a small difference in frequency and mode shape but a huge difference in damping ratio. Even if the tolerance value of the damping ratio is relaxed, the following situations still occur in many cases. Taking a certain bridge as an example, the SSI method is used to identify the mode through the response data of the complicated train and wind load. After the frequency stability has been determined, the first column of the damping ratio matrix is shown in Figure 2. Figure 2 shows that even if the frequencies of these stable points are considered the same frequency and most of the damping ratios fluctuate near a certain value within the allowable range of the frequency resolution, some points with damping ratios still have obvious deviations. The average damping ratio is obtained by dividing the sum of all samples by the total number of samples. If these outliers are not eliminated, the average damping ratio of the column is calculated to be 2.813%, so the column damping ratio does not satisfy the stable condition of DT \ 40%. The frequency of this order is considered a false mode and consequently eliminated, which results in modal omission.
However, these results are obviously unreasonable. This paper introduces Grubbs criterion to detect outliers to effectively eliminate outliers in the damping ratio matrix and improve the ability of the stability graph to determine the order of the system. Grubbs criterion test for outliers has two characteristics: (1) It is suitable for the case where the number of sample data points is 20-100, which has typical meaning in the stability graph method; (2) The critical value for outlier discrimination is unrelated to the average and standard deviation of the sample data, and it is easy to adjust and control.
Suppose that Y i is the sample data to be detected, where i = 1, 2, Á Á Á , N. Statistics G and critical value G a in Grubbs criterion can be expressed as follows: where Y is the average value of the sample data; s is the standard deviation of the sample data; a is the significance level; t ½a=2N, NÀ2 is the critical value of the t distribution when the degree of freedom is N À 2 and the significance level is a=2N. When G . G a , the corresponding Y m will be determined as an outlier and deleted from the sample. Since Grubbs criterion can only determine one outlier at a time, the sample data with Y m deleted as the original data must be continuously repeated until no outliers can be found.
In this paper, a = 0:05 is used in the discrimination process as the modal parameters are usually identified according to this value in practice. Therefore, the deleted value has a probability of P = 1 À a = 95% as an outlier. The scatterplot of the damping ratio obtained after removing the outliers through the Grubbs criterion is shown in Figure 3. In Figure 3, after removing two outliers with large deviations in the damping ratio, the average damping ratio can better represent the distribution of the column damping ratio.

Dynamic characteristic calculation
We use Oujiang Bridge as an example. Oujiang Bridge is a concrete cable-stayed bridge with two towers and two cable planes. The span layout is (52 + 90 + 300 + 90 + 52) m, and the bridge length is 584 m. The main beam adopts a single-box double-chamber equalheight prestressed concrete box beam with a full cross-  section of 13 m and a height of 4.0 m at the center. The cables of the bridge adopt a double-cable plane system and are arranged in a fan shape. The main tower adopts a diamond shape. The tower above the deck adopts an inverted Y shape, and the tower column below the deck shrinks into a diamond shape. The total height of the pylon above the bottom of the tower is 112.5 m and 118 m, and the tower height above the beam top is 75 m. The main content of this paper is the modal parameters of the main girder. Figure 4 shows the main mode and frequency of Oujiang Bridge. The sampling frequency of each measuring point is 100 Hz, and the sampling time of each signal is 50 s. Figure 6 shows the acceleration response curves of the middle of the mid-span (measuring point 4) under two conditions. Due to space limitations, the acceleration response curves of the remaining measuring points are no longer drawn one by one.

Modal parameter identification
The above improved stochastic subspace identification is programmed in MATLAB, and the vertical bending, lateral bending, and torsion modal parameters of the Oujiang Bridge are identified. When the vertical and lateral bending modal parameters are identified, the acceleration response signal of the bridge under the action of fluctuating wind loads is used after the train has crossed the bridge. When the torsional modal parameter is identified, the torsional mode is not obvious under the excitation of only fluctuating wind, so the analysis signal uses the acceleration response signal of the bridge under the combined action of a fluctuating wind load and a vehicle load. In the calculation, the number of rows of Hankel matrix is selected as 2I = 200, and the number of columns is J = 4000. The minimum and maximum orders of the system matrix are N min = 2 and N max = 100. The frequency resolution is 0.1 in the stability graph criterion, and the minimum number of frequency stable points is 30. The damping ratio tolerance is 0.4, and the similarity factor of the  maximum mode is 0.05. The significance level of Grubbs criterion is a = 0:05 when the outliers of the damping ratio are eliminated.
The obtained vertical bending, lateral bending and torsion modal parameters are shown in Table 1. The stability diagram is shown in Figure 7. The main mode shape curves are shown in Figures 8 to 10 after three sample curves have been fitted.
The identified frequencies and mode shapes are consistent with the finite element results from the results of the stochastic subspace modal parameter identification by Grubbs criterion. The first-order vertical bending  appears as false modes that originally belonged to the longitudinal drift of the main beam, but there are no false modes or mode omissions in the lateral bending and torsion. The stability graph method has the advantages of system order determination and removal of false modes. In addition, because the stochastic subspace identification is completely data-driven, it must only determine a few parameters to obtain a better, more reliable result through convenient calculations.

Conclusion
During the test, it is very important to identify the modal parameters of each order in a certain frequency range, where the actual vibration response of the structure can be predicted. The specific conclusions are as follows: 1. The stability diagram method can effectively determine the order of the system, and proper adjustment of the parameters in Grubbs criterion can eliminate false modes and avoid modal omissions. 2. The modal parameters analyzed by the stochastic subspace algorithm with the introduction of Grubbs criterion can effectively improve the accuracy of modal parameter identification of the system. 3. In the improved stochastic subspace identification algorithm, the average damping ratio in the first column of the ratio matrix of bridge damping changes from 2.813% to 1.486%, which improves the accuracy of the stability diagram method and can be better applied to identify the modal parameters of bridges subject to environmental excitations.
Although the improved algorithm has achieved some excellent results in wind tunnel test, damping ratio identification is still a difficult point in modal   parameter identification. Damping ratio identification algorithm should be further improved in subsequent research.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.