Output ‐ Only Modal Estimation Using Eigensystem Realization Algorithm with Nonstationary Data Correlation

: The conventional eigensystem realization algorithm with data correlation (ERA/DC) combines the impulse response or free response data of a structural system with the concept of correlation function to identify the modal parameter of the structural system. Previous studies have shown that the modal parameters of structural systems subjected to stationary white noise excitation can be estimated by ERA/DC from the ambient response without excitation data. This concept is extended in this paper for output ‐ only modal identification for the structural system with complex modes under ambient excitation as a nonstationary process in the form of a product model. Numerical simulations and experimental verification are used to validate the effectiveness of the proposed method for response ‐ only modal estimation, and the stabilization diagram is used with modal assurance criterion (MAC) to distinguish structural modes from fictitious modes.


Introduction
The modal estimation theory uses free response or impulse response data to identify modal parameters, or uses excitation and response data simultaneously for parametric estimation of structures. However, in many actual engineering structures, the structures are often subjected to ambient excitation, which is difficult to be available. Therefore, how to estimate the modal parameter of structural systems without excitation data for structural control, damage detection, or analysis, and model updating is an important issue in structural engineering.
The output-only modal identification can be implemented by multiple time-domain or frequency-domain methods. In 1965, Ho and Kalman [1] proposed the minimum realization algorithm, using the Markov parameter composed of the impulse response function to obtain the minimum-order state-space representation. This method could describe the system accurately but disregarded effective estimation of impulse response with noise. In 1974, Zeiger and McEwen [2] proposed the concept of Singular Value Decomposition (SVD) in conjunction with the minimum realization algorithm. In 1986, Juang and Pappa [3] developed the eigensystem realization algorithm (ERA) based on SVD and the minimum realization algorithm. The ERA is mathematically reasonable, but improper optimum dimension and size of the Hankel matrix will result in incorrect modal parameter estimation [4]. In 1988, Juang et al. [5] proposed the improved method of ERA, that is, the eigensystem realization algorithm using data correlation (ERA/DC), which reduces the effect of noise on modal estimation using the characteristics of a correlation matrix to increase the accuracy of identification. In 1993, James et al. [6] assumed the ambient excitation to be a white noise and calculated the autocorrelation and cross-correlation functions of the response signal. They found that the mathematical version of the correlation function was identical to that of free response and impulse response and the modal parameter of the structure was obtained. Chiang and Lin [7] used the correlation between structural system response signals corresponding to stationary white noise to redefine the Markov parameter and built the generalized Hankel matrix with SVD to obtain the system modal parameter. In recent years, the ERA has been extensively used in output mode estimation [8,9]; this method has been effective in identifying the modal parameter of large-scale structures, such as offshore structural systems [10] and roof overflow powerhouses [11]. In 2011, Li et al. [12] studied the Hankel matrix of the ERA and found that a larger Hankel matrix made the decomposition of noise much easier. The numerical calculations may be erroneous if the difference between values of rows and columns is too large. Caicedo [13] used a four-layer rectangular steel structure (benchmark structure) for the experiment in which the ambient vibration was measured; six modes were identified successfully by using the natural excitation technique in conjunction with eigensystem realization algorithm (NExT-ERA), including a group of similar modes. It was indicated that the number of sampling points of the correlation matrix would influence the identification result directly. Liuchong and Caiyou et al. [14] measured the rail vibration signals induced by impulse and actual train running and used the ERA for identification. The impulse and ambient vibration results were close to the computer-aided analysis result. Siringoringo and Fujino [15] used this method for output mode estimation of a suspension bridge. They identified low-frequency closely spaced modes effectively, but the high-frequency closely spaced modes could not be identified due to noise interference. Kordkheili et al. [16] used the method for mode estimation for the flat plate with similar modes and the pulley with repetition frequency mode.
The ERA/DC combines the impulse response or free response signal of a structural system with the concept of correlation function to identify the modal parameter of a structural system [7]. We can determine the modal parameter of the structural system by ERA/DC from the environmental response signal directly, without excitation signal measurement; the only applicable assumption is that the excitation is a stationary white noise. This may not match the actual environmental condition. This concept is extended in this paper to study how to perform modal parameter identification for the linear structural system with complex mode under ambient excitation with a nonstationary process. We will combine the ERA/DC with the correlation function of nonstationary ambient vibration to develop a complete theory for the nonstationary in the form of a product model. It is combined with appropriate channel-expansion technique to construct the generalized Hankel matrix containing adequate dynamic information for mode estimation, to extend the applicability of ERA/DC to nonstationary output-only structural modal estimation.

Research Method
The equation of motion of a Degree Of Freedom (DOF) linear vibration system with nonstationary excitation can be expressed as (1) where , , and are the mass, damping, and stiffness matrices of the vibration system, respectively, is a nonstationary white noise with the product model, , is the stationary displacement response when the system is excited by a stationary white noise , is a deterministic envelope function. If is a time-varying function with the secular change, meaning 0 and 0, then and , where and are the response signals of the system corresponding to stationary white noise , the system response can be expressed as a nonstationary model coincident with envelope function of nonstationary excitation, and Equation (1) can be rewritten as: (2) In practical application, the dynamic response data of structures are in discrete form, so the continuous state equation is transformed to discrete state space: where , , and are state, output, and input vectors, respectively. , , and are the system, input, and output matrices, respectively.
is the sampling point sequence. As is a nonstationary white noise with zero mean, the Hankel matrix is formed of the system output vector . The nonstationary correlation function matrix between response data, in different time-delayed times, τ is defined as: According to Equation (3) and Equation (4), can be expressed as where is the covariance matrix. In , is the stationary white part of nonstationary white noise in the form of a product model, is the displacement response when the system is excited by a stationary white noise. As is uncorrelated with , 0. If the excitation can be expressed as the product model of the nonstationary white noise, and the nonstationary response of the system can be expressed as a nonstationary process coincident with the envelope function of excitation, the nonstationary correlation function matrix can be treated as an approximately stationary correlation function matrix. In the analysis process, if the stochastic process can be reasonably assumed to be the whole process, the statistical value of correlation can be derived from a single (delay is long enough) sample function, so for , the calculation can be approximated as where and correspond to the and of traditional ERA/DC. To reduce the effect of noise and to obtain a relatively accurate mode identification result, the data correlation matrix is used to build the generalized Hankel matrix : where is defined as the generalized observation matrix, is the generalized control matrix, is the channel number of the measurement response signal. and are the extended channel number of Hankel matrix and generalized Hankel matrix , respectively. To make the number of identified mode not less than the number of mode to be identified, the sum of and a shall be larger than or equal to the mode number to be identified; and determine how many correlation matrices of different time differences are to be analyzed in the generalized correlation Hankel matrix. After the new data correlation matrix is defined, the theoretical structure of ERA/DC is used to obtain the state matrix: In Equation (10), the system matrix A contains the characteristics of the system, so eigenvalue decomposition is performed for A, eigenvalues , , … , and the corresponding eigenvectors ( , , … , ) can be obtained, and the characteristic matrix can be expressed as The modal coordinate transformation relation is (13) Equation (13) is substituted in Equation (6) to obtain (14) where is the modal output matrix. As the discrete data are used for calculation in the process, the obtained eigenvalue shall be transformed from discrete system eigenvalue into the eigenvalue of the continuous system, and the transformation relation is expressed as follows: where is the real part of eigenvalue, is the imaginary part of eigenvalue, Δ is the sampling time; the natural frequency and damping ratio can be derived from eigenvalue : According to a prior reference [7], the ERA/DC has been able to redefine the Markov parameter according to the data correlation between the ambient responses of structures subjected to stationary white noise. The generalized Hankel matrix is built and combined with SVD and then the system modal parameter can be obtained, whereas the traditional ERA/DC is only applicable to the system impulse response. This concept is extended in this study and it is mathematically proven that if the ambient excitation on the system can be properly expressed as the nonstationary process in the form of a product model with a slow-varying envelope function, the data correlation analysis between nonstationary excitation and the corresponding response can be approximated as stationary correlation calculation. Afterward, with appropriate extended channel technology, a generalized Hankel matrix containing adequate dynamic information is built for mode estimation to extend the applicability of ERA/DC to nonstationary output-only structural modal estimation.

Results and Discussion
In a dynamic test for the structure with external excitation, the modal parameter can be estimated according to the measured data of excitation and response of the structural system. However, it is sometimes difficult to perform a dynamic test for a large structure in practice. It is unlikely to measure the complete modal information of the actual structure. Therefore, to confirm the effectiveness of the developed algorithm theory, numerical simulation and tests are usually used to validate the feasibility of the proposed method.

Six-DOF Chain Model of a Cantilever Beam
A six-DOF linear chain model system of a cantilever beam is shown in Figure 1. The mass matrix , stiffness matrix , and damping matrix of the system are expressed as follows: The damping matrix of the system is proportional damping (i.e., the linear combination of mass matrix and stiffness matrix ). It should be mentioned that the proportional damping matrices are suitable for modeling the behavior of most structural systems, in which the damping mechanism is distributed rather uniformly throughout the structures. The six DOF are excited by nonstationary white noise using the Newmark method and the initial state is stationary. To obtain the displacement response of the time domain, the time-domain sampling interval is 0.01 s, the number of sampling points is 131,072.
The nonstationary white noise discussed in this paper will be imported into the nonstationary process in the form of a product model: is a stationary white noise. The nonstationary white noise used in this paper, as shown in Figure 2, is composed of a stationary white noise with power spectrum density 0.02 rad s ⁄ ⁄ and envelope function Γ 4 . .
, as shown in Figures 3 and 4.    In the course of system identification and to avoid omission of modal identification, a higher order is used for identification. There will be a fictitious mode in the process. Therefore, we can employ the values of modal assurance criterion (MAC) to check the agreement between the identified and the theoretical mode shapes, and use the stabilization diagram to filter the fictitious mode and to decide the order for calculation. Figure 5 is the stabilization diagram of the six-DOF chain model system. It is roughly observed that the frequencies 5, 13, 20, 27, 31, and 33 rad/s can be identified steadily when the order number is 6. The six modes of corresponding frequencies are identified as system modes, and the identification order number α is 6. The displacement response of each DOF is shown in Figure 6. Table 1 shows the result of modal parameter identification using ERA/DC. As seen, the frequency errors are less than 1%, and the maximum error of damping ratio is 5.92%. Figure 7 shows the identified mode shapes, which are approximately coincident with theoretical mode shapes, and the MAC values are larger than 0.99, meaning the method is effective on identification.

Six-DOF Chain Model with Heavy Damping
To further validate the reliability of this method, this paper uses the Group 2 numerical model which is also a six-DOF chain model, as shown in Figure 8 Figure 9 is the stabilization diagram of the system, roughly showing five frequencies of about 8, 19, 27, 31, and 43 rad/s. It is observed that the identification state of the mode of 31 rad/s is relatively unstable before order 10, so the identification order α is 12. The modal parameter identification can be performed effectively in this order, and the results are shown in Table 2 and Figure 10. Table 2 shows the result of modal parameter identification using ERA/DC for a six-DOF chain model with heavy damping. The sixth mode has a higher error in identification of natural frequency, which is 14.6%, the others are lower than 2%, and the damping ratio errors are less than 30%. Figure 10 shows the corresponding exact and identified mode shapes of this structure; the fourth measurement DOF in the fifth mode shows phase contrary sign. The computed MAC value of the fifth mode is lower than 0.70, and the rest are approximately coincident with exact mode shapes.

Seven-DOF Car Model
To demonstrate this method applicable to complex structure systems, a car model with two groups of closely spaced modes is considered, as shown in Figure 11. The system is a seven-DOF system, including car body bounce, pitch, and roll. Each of the four wheels represents a DOF and has the dynamic behavior of bounce; related parameters are shown in Table 3. The mass matrix , stiffness matrix , and damping matrix of the system are expressed as follows.
In this system, as the external force on the actual car body is close to stationary white noise, the DOF of the car body is excited by stationary white. The DOF of four wheels is excited by nonstationary white noise, the time domain sampling interval is 0.01 s, and the number of sampling points is 131,072.   Table 4 and Figure 13. Table 4 shows the result of modal parameter identification for the seven-DOF car model by ERA/DC. The frequency errors are less than 5%, and the damping ratio error is less than 30%. Figure 13 shows the mode shape identification result. The first mode is pitch mode, the second mode is car body bounce mode, the third mode is roll mode, the fourth and fifth modes are front-wheel bounce mode, the sixth and seventh modes are rear-wheel bounce mode, and the identification result is approximately coincident with theory. Distance from center of gravity to fore axle 1.5415 Distance from center of gravity to rear axle

Twenty-DOF Chain Model with Many Groups of Closely Spaced Modes
To further examine the effectiveness of the present method for a structural system with more DOFs, we conduct output-only modal estimation using a 20-DOF chain model rigidly fixed at both ends containing the mechanical properties of each 1kg mass and all 600 N/m spring constants. The damping matrix is assumed to be proportional to a combination of the mass and the stiffness matrices as given by 0.05 0.001 / . This structure is subjected to a vibration practically recorded at Sun-Moon Lake on September 21, 1999 when the Chi-Chi Earthquake with a moment magnitude of 7.6 occurred in central Taiwan. The sampling interval and period of this seismic record are Δt = 0.005 sec and T =59.995 sec , respectively. A sample of the seismic record, which serves as the excitation input acting on the sixth mass of the model, is shown in Figure 14. Figure 15 is the stabilization diagram of the system, roughly showing the first 12 natural frequencies of the identified structural modes. The system considered has many groups of closely spaced modes, as listed in Table 5. The results of identification are also summarized in Table 5, which shows that the errors in natural frequencies are less than 5% and the errors in damping are larger. From Table 5, 12 out of the 20 vibration modes are identified more accurately with MAC ≥ 0.9 as the threshold. Because the structural response generally has lower sensitivity to damping ratios and mode shapes than to the natural frequencies, the errors of identified damping ratios and mode shapes are somewhat larger. Since the contribution of the high-order vibration modes to the structural response is somewhat less than that of the low-order vibration modes, the high-order vibration modes are not identified as accurately as the low-order vibration modes in general.

Experimental Validation of Free Beam
To validate the effectiveness of the method proposed in this paper, an actual beam structure of free boundary is used for the experiment as shown in Figure 16. The Brüel & Kjaer RT Pro Photon 7.0 data acquisition system, PCB 208C01 force sensor, and PCB 352B10 piezoelectric accelerometer are used for measuring response signal. The simulated nonstationary excitation is imported into the Modal Shop K20070E01 vibration exciter through Teledyne LeCroy T3AFG40 signal generator and the vibration exciter excites the free beam structure. Finally, the excitation signal and response signal are obtained by force sensor and accelerometer, as shown in Figures 17 and 18. Currently, a roving accelerometer testing based on output-only modal estimation is implemented to extract the practical modal properties of the beam structure; nine measurement positions on the realistic aluminum alloy beam were marked and the shaker excitation impacts acted as the middle (fifth) location of the beam.
The length, width, and height of the beam structure used for this paper are 300 mm, 20 mm, and 15 mm, respectively, the mass is 245.7 g, and the material is 6061-T6 aluminum alloy. According to the theory of engineering vibration, the natural frequencies of the first three groups of modes are about 865.42 Hz, 2385.61 Hz, and 4676.75 Hz, as listed in Table 6. The data are used as a reference for experimental modal analysis (EMA). The natural frequencies of modes and mode shapes match the analytic solution, meaning the experimental results are reliable. For the EMA procedure, the Brüel & Kjaer RT Pro Photon 7.0 data acquisition system is employed to respectively extract the excitation data induced from a 086C01 PCB impulse hammer impacting on the structure as well as the response data obtained from a 352B10 PCB piezoelectric accelerometer, and the frequency response function of each measurement degree of freedom can be obtained through roving accelerometer testing to perform modal estimation.
This paper uses the modal parameter obtained by EMA as a reference frame and uses the ambient response to compare the identification results of ERA/DC, as shown in Table 6 and Figure 19. It is observed that the frequency errors are less than 4%, and the maximum error of damping ratio is 48.52%. Figure 19 shows the identified mode shapes which are approximately coincident with theoretical mode shapes, and the modal assurance criteria (MAC) are larger than 0.77. This means that the method is effective on identification in practical application.

Discussion
The system order α is selected, which refers to the stabilization diagram in this study. The number of frequencies in the signal is evaluated by spectrum and whether the order is sufficient for identification is judged according to the stabilization diagram, but only the frequency can be used as a reference. The accurate damping ratio or mode shape cannot be identified by the order. Reference [6] indicates that parameter α (number of columns of Hankel matrix) and β (number of rows of Hankel matrix) shall be intercompared, and the difference shall not be too large. Therefore, the parameter β must be adjusted appropriately when selecting α.
The ERA/DC is likely to be influenced by noise on the identification of high-order vibration mode, and closely spaced modes are unlikely to be identified [12]. The closely spaced modes of the model selected in this study exist in relatively high-order vibration modes. The numerical simulation process is free of additional noise factors in this study, and the damping ratio of high-order closely spaced vibration modes still has a relatively high identification error, so it is not completely ascribed to noise.

Conclusion
The concept of ERA/DC is extended to nonstationary ambient response in this study. If the ambient excitation on the system can be appropriately expressed as the nonstationary process in the form of a product model with a slowly-varying envelope function and stationary white noise, the system response can be expressed as a nonstationary process coincident with the envelope function of external force. The data correlation between nonstationary excitation and corresponding response can be theoretically shown approximate to stationary data correlation. Afterward, the appropriate extended channel technology is used to build the generalized Hankel matrix containing adequate dynamic information for mode estimation, to extend the applicability of ERA/DC to nonstationary output-only structural modal estimation. The findings of this paper can be concluded as follows.
(1) This paper deduces that the correlation matrix of excitation and response signals can be made into the form identical with stationary white noise under nonstationary one in the form of a product model. With this characteristic, the ambient response can be directly used for ERA/DC by using channel-expansion technique. It is not necessary to transform ambient response into free-response through additional data processing before modal identification, preventing unnecessary distortion in vibration data to lead errors in parametric estimation of structural system..
(2) The ambient excitation considered in this paper can be treated as a nonstationary process, which is the chosen nonstationary white noise in the form of a product model of a simple mathematical form and convenient for processing. According to numerical simulation, the structural mode can be identified effectively in most cases. Also, if the ambient excitation meets the characteristic of periodic vibration, the method proposed in this paper can be used for identifying modal parameter.
(3) The ERA/DC is performed directly using ambient response data only for modal parameter identification through the method proposed in this paper, the natural frequency identification result is good. In terms of the identified mode shape, as the response data remains complete, the MAC value approaches 1 meaning the agreement with the actual mode shape is good. In terms of damping ratio, the error changes significantly, adequate response points and time-delayed sampling is required to augment system matrix order and to increase the damping ratio identification accuracy.
(4) As the nonstationary response only is directly used for modal estiamtion, the identification result may contain structural modes, excitation modes, and fictitious modes. The fictitious mode is eliminated by using the concept of stabilization diagram, and the agreement between the identified mode shapes and actual structural mode shapes is checked by using modal assurance criterion.