Automated operational modal analysis (OMA) is attractive and has been extensively used to replace traditional OMA, which requires much empirical observation and engineers ’judgment. Still, the uncer- tainties on modal parameters and spurious modes are challenging to estimate under the field condi- tions. For addressing this challenge, this research proposed an automated modal identification approach. The proposed approach consists of two steps: (1) modal analysis using covariance-driven stochastic subspace algorithm; (2) automated interpretation of the stabilization diagram. An additionaluncertainty criterion is employed to initially remove as many spurious modes as possible. A novel threshold calculation for clustering is proposed with incorporating uncertainty of modal parameters and the weighting factor. An improved self-adaptive clustering with new distance calculation is usedto group physical modes, followed by the final step of robust outlier detection to select outlying modes. The proposed automated approach requires minimum human intervention. Two field tests of the footbridge and a post-tensioned concrete bridge are used to verify the proposed approach. A modal tracking was used for continuously measured data for demonstrating the applicability of the approach. Results show the proposed approach has fairly good performance and be suitable for auto-mated OMA and long-term health monitoring.ARTICLE HISTORY Received 7 September 2020 Revised 16 January 2021 Accepted 2 February 2021 KEYWORDS Automated interpretation;clustering; long-term health monitoring; minimumhuman intervention;Operational modal analysis; threshold calculation; uncertainty criterion 1. Introduction Over the last few decades, long-term structural health moni- toring (SHM) using vibration responses has been applied tocivil infrastructures such as bridges and buildings (Chen & Ni,2018). It is widely recognized that continuous monitor- ing is useful in tracking changes in dynamic modal parame-ters: frequencies, damping ratios, and mode shapes.Inparticular, operational modal analysis (OMA) is paid greatattention because it can be implemented efficiently, econom- ically, and safely, avoiding interruption to the normal oper- ation of observed structure and requiring no artificialloadings (Brownjohn, Magalhaes, Caetano, & Cunha, 2010). Therefore, only ambient excitation (i.e., wind and traffic) isused to identify dynamic modal parameters by measuring vibration responses. Moreover, dynamic modal parameters are beneficial for: (1) finite element model updating (Zeng& Kim, 2020); (2) vibration-based damage detection (Honfi, Bj€ornsson, Ivanov, & Leander, 2020); (3) real-time alarm systems for bridges and buildings (Chen et al., 2011). It can be incorporated with maintenance systems to observe variations of modal parameters over time and trackabnormality to prevent disastrous failure at an early stage.However, it requires a huge amount of recorded data anddata analysis within a short period of time. Therefore, it is labor-intensive, even impractical, to process massive meas- ured data and identify modal parameters by manual inter-vention and engineers ’experience. Much user intervention on a large amount of vibration data can be an obstacle in areal application. For overcoming this, it is essential to have an automated evaluation of structural conditions in almost real-time. In recent years, there are growing researches focusing on the development of automated OMA techniques. Generally,OMA techniques can be categories into two groups: directmethods and indirect methods. Direct methods measurevibration responses of observed structures to identify modal parameters, requiring networks of sensors installed on the target structures. Many system identification algorithmshave been developed to be compatible with direct methods,such as stochastic subspace identification (SSI) (Cabboi,Magalh ~aes, Gentile, & Cunha, 2017), eigensystem realization algorithm (Yang, Yi, Qu, Li, & Liu, 2019), poly-least squares complex frequency domain method (Poly-MAX)(Magalh ~aes, Reynders, Cunha, & De Roeck, 2009), fast Bayesian FFT method (Au, 2011b ), etc. On the other hand, indirect methods are referred to as the primary use of mobile sensors installed inside in thepassing-by vehicles to perform modal identification. No sen-sors are needed to be deployed on structures as the indirectmethods. Recently, much progress has been made in termsof system identification using networks of mobile sensors.For example, mass-modeled vehicle scanning method forbridges (Yang, Lin, & Yau, 2004; Yang & Sun, 2020), hybrid sensors based method (Marulanda, Caicedo, & Thomson, 2017), extended structural identification using expectation optimization (STRIDEX) (Matarazzo & Pakzad, 2018), matrix completion method (Sadeghi Eshkevari, Pakzad, CONTACT Young Hoon Kim yhkim005@louisville.edu /C2232021 Informa UK Limited, trading as Taylor & Francis GroupSTRUCTURE AND INFRASTRUCTURE ENGINEERING https://doi.org/10.1080/15732479.2021.1919151 Tak/C19a/C20c, & Matarazzo, 2020), crowdsensing based method (Matarazzo et al., 2020). A detailed literature review on dir- ect and indirect methods can be found in Malekjafarian,Mcgetrick, and Obrien ( 2015). The present work belongs to the class of direct methods. Vibration response is measured from sensors installed oninvestigated structures. Among various system identificationalgorithms in the direct methods, SSI has been widelyapplied to diverse structures to perform modal analysis. It offers accurate identification results and simple implementa- tion, which are important attributes accounting for its popu-larity. In addition, due to SSI ’s explicit mathematical nature, SSI tends to be more suited for automated modal identifica-tion. However, the major challenge in SSI is spurious modesappear in outputs. Commonly, spurious modes consist of pure mathematical (i.e., non-physical) and noise modes (Reynders, Pintelon, &De Roeck, 2008). The most common strategy to deal with this challenge is to construct a stabilization diagram, a plotof model order vs. frequency for an extensive range ofmodel order. In the stabilization diagram, physical modes are referred to as those poles that cross most of the model orders consistently. Therefore, physical modes should begraphically recognized and homogeneously distributed alongvertical alignments in the stabilization diagram (Cabboi etal.,2017). On the contrary, spurious modes appear in the stabilization diagram in a scattered way. Spurious modes areeliminated in a manual analysis depending on empirical dis-covery and engineers ’judgment, which is subjective, time- consuming, and leads to incorrect modal identification. For addressing this issue, a variety of methods are pro- posed in the literature to automatically interpret stabilization diagrams and remove spurious modes. In general, the pro- cess can be divided into three steps: 1. Step 1: Apply the modal validation criteria to eliminate as many spurious modes as possible in the stabiliza-tion diagram 2. Step 2: Group modes with similar characteristics, i.e., frequencies, damping ratios, and mode shapes by clus-tering strategies 3. Step 3: Detect outliers in each cluster to improve the accuracy of modal parameters and select representativeof each cluster Several methods aiming at minimizing human involve- ment in the interpretation of the stabilization diagram havebeen developed. For example, in step one, many modal val-idation criteria are proposed to detect spurious modes inthe stabilization diagram. These criteria include hard crite-ria, which yield a binary answer, and soft criteria, whichyield a certain range of values. Reynders, Houbrechts, andDe Roeck ( 2012) thoroughly reviewed and summarized hard and soft criteria. However, conventional modal validationcriteria limitedly remove a certain number of spuriousmodes; many spurious modes, which still remain in the sta-bilization diagram, affect parameter estimates ’accuracy andimposes a computational burden to the following step (clus- tering process). In step 2, various clustering strategies are widely employed to group modes with similar characteristics. Hierarchical clustering has been extensively applied by many researchers and is considered as the most natural approach (Reynders et al., 2012). Hierarchical clustering has a significant advantage of allowing a good selection of phys- ical clusters. However, the main drawbacks include a user- defined tree cutoff distance and human intervention withdemanding computational cost. Furthermore, the hierarch- ical clustering is sensitive to outliers. Another popular clus- tering strategy is partitioning methods, often referred to as K-means clustering (Neu, Janser, Khatibi, & Orifici, 2017). K-means clustering has the benefit of being fast processing. However, the number of clusters has to be predefined, and it is sensitive to cluster seeds (initial centroid). By merging the benefits of hierarchical clustering and K-means cluster- ing to overcome some of their limitations, self-adaptive clus- tering is recently proposed (Cabboi et al., 2017; Fan, Li, & Hao, 2019). The self-adaptive clustering has outstanding fea- tures: (1) simple implementation; (2) fast computation; (3) no need for the number of clusters; (4) clustering thresholdis iteratively trained during the clustering process. While it still starts with a user-defined clustering thresh- old, which requires some level of human intervention. Somemethods are proposed to automatically calculate clustering threshold based on statistical properties, i.e., mean and standard derivation or median, of the distance between two closed poles in the stabilization diagram (Magalh ~aes et al., 2009; Reynders et al., 2012; Yang et al., 2019). However, these methods do not consider uncertainty on modal parameters and inaccuracy of mode shapes. In practice, modal parameters ’ uncertainty is inevitable due to modeling error and measurement noise; it can be a more reasonable approach to consider uncertainty when calculating the clus- tering threshold. Also, measurement on mode shapes is less accurate than that on frequencies. Thus, a weighting factor can reduce the inaccuracy of mode shape difference onthreshold calculation (Boroschek & Bilbao, 2019). In step 3, some outliers are undesirably involved in iden- tified physical clusters; this phenomenon is pronounced in adamping ratio with a scattered nature. Most outlier detec- tion techniques need to define limit bounds, such as box- plot outlier detection (Yang et al., 2019). A bound-free out- lier detection method is needed to improve the accuracy of parameter estimations. In summary, challenges to the cur-rent automated interpretation of the stabilization diagram are listed as follows: 1. Conventional modal validation criteria are inefficient resulting in high computational cost in the cluster- ing process. 2. The clustering threshold and distance calculation in the clustering process does not consider the uncertainty of parameters and the weighting factor. 3. Uncertainties on identified modal parameters and phys- ical clusters are unavailable.2 J. ZENG AND Y. HOON KIM 4. Outlier detection requires to define limit bounds. This article attempts to address the aforementioned chal- lenges. A two-stage framework for automated operational modal identification is proposed. Figure 1 shows the flow- chart of the proposed framework: (1) modal analysis using covariance-driven reference-based SSI (SSI-cov/ref); (2) two- stage automated interpretation of stabilization diagram. Inthe first place, SSI-cov/ref is adopted to perform modal ana- lysis and construct a stabilization diagram. Subsequently, a two-stage automated analysis for the stabilization diagram iscarried out. At the pre-processing stage, besides applying conven- tional modal validation criteria, such as damping ratio checkand modal complexity check, to eliminate spurious modes, a new supplementary criterion: uncertainty criterion, which is also applied for further removal of spurious modes. At theclustering stage, a novel threshold calculation, which incor- porates the uncertainty of modal parameters and weighting factor, is proposed. An improved self-adaptive clusteringwith new distance calculation is then employed to group modes with similar characteristics and identify physical clus- ters. Finally, robust outlier detection is implemented toexclude outliers. The average of each cluster ’s elements is chosen as representative frequency, damping ratio, and mode shape. Two field tests were used to evaluate the proposed approach. The first field example is the Dowling Hall Footbridge located at Turfs University, a two-span steelframe bridge with 144 ft (44 m) long and 12 ft (2.7 m) width having a reinforced concrete deck. The second field test example is the Z24 bridge located in Switzerland, a post-ten-sioned concrete bridge with a main span of 100 ft (30 m)and two sides span of 46 ft (14 m), considered as a bench- mark in the research community. This article is organized as follows: In Section 2, the background of SSI-cov/ref is briefly introduced. In Section 3, a two-stage approach for proposed automated modal identification is presented. In Section 4, the capability of theproposed approach is validated by two field tests along withthe modal tracking results. Finally, the conclusion is pre- sented in Section 5. 2. Background of SSI SSI has been extensively spread over the field of OMA dur-ing the past few decades accounting for its quick implemen-tation and high accuracy. In this article, covariance-drivenreference-based SSI (SSI-cov/ref) is employed to reduce thedimensions of the output matrix and computational cost.The detailed theoretical fundamentals of SSI-cov/refl arefully described in the literature (Peeters & De Roeck, 1999). Briefly, SSI-cov/ref is developed based on assuming a linear and stationary Ndegree of freedoms (DOFs) system with a dynamic motion characterized by the discrete-time state-space equation: x kþ1¼Axkþxk yk¼Cxkþvk(1) where subscript kdenotes time step; A2Rn/C2ndenotes sys- tem state matrix with ( n¼2N);C2Rl/C2nis an output matrix, lis defined as the number of measured signals; xk2 Rn/C21andyk2Rl/C21are discrete-time state vector and meas- ured response vector, respectively; xk2Rn/C21is a process white noise vector. vk2Rl/C21is the measurement of the white noise vector. Peeters and De Roeck ( 1999) reported that modal parameters can be identified from system matri-cesAandC: Two main SSI preparation parameters significantly affect the accuracy of identification results: (1) model order; (2)time lag, i:Unfortunately, the value of model order and i, which yield the best identification results are never known(Fan et al., 2019; Ubertini, Gentile, & Materazzi, 2013). In practice, it is necessary to over-specify model order to cover Figure 1. A framework of the proposed automated approach. Figure 2. A flowchart of the proposed two-stage approach of automated modal identification.STRUCTURE AND INFRASTRUCTURE ENGINEERING 3 weakly-excited modes, but spurious modes increase with model order increases. These spurious modes must besingled out in the subsequent procedure. On the other hand, the value of idetermines the size of the response covariance function. The smaller imay fail to identify the fundamental mode, but the larger value of iyields more spurious modes and increases computational time. The value of imay be chosen at least estimated value as follows (Fan et al., 2019): i/C21T i=t (2) where Tidenotes fundamental period (s); tdenotes sam- pling interval (s). 3. A Two-stage automated modal identification In this section, a two-stage framework for automatedmodal identification is proposed. The flowchart of the entire automated process in detail is presented in Figure 2 . Section 3.1 presented the pre-processing stage for auto- mated identification, including conventional modal valid-ation criteria ( Section 3.1.1 ) and a new additional uncertainty criterion ( Section 3.1.2 ). Subsequently, the clus- tering stage is introduced in Section 3.2 .An e w l yp r o p o s e d threshold calculation for clustering is provided in Section 3.2.1. An improved self-adaptive clustering is then employed to determine physical clusters in Section 3.2.2 . Finally, robust outlier detection is performed to improve the accuracy of modal parameter estimates in Section 3.2.3 . The pseudocode of the proposed automated framework isalso provided in Appendix A . 3.1. The pre-processing stage A full stabilization diagram looks very complex and contains the spurious modes. At the pre-processing stage, differentvalidation criteria are utilized to initially screen spurious modes, clarifying the stabilization diagram and consequently accelerating automated interpretation of the stabiliza-tion diagram. 3.1.1. Modal validation criteria Firstly, for civil engineering structures, a negative or high damping ratio hardly appears in practice. Therefore, adamping ratio with less than 0 and higher than 10% is dis-carded (Cabboi et al., 2017; Fan et al., 2019). Additionally, two popular modal validation criteria are used to measure the complexity of mode shape vectors, namely, modal phase collinearity (MPC) and mean phase deviation (MPD). Thesetwo indicators have been utilized by various researchers to distinguish physical modes from spurious modes (Neu et al.,2017; Reynders et al., 2012). The real (Re) and imaginary (Im) part of mode shapes display a linear correlation, which can be assessed by the MPC indicator. The value of MPCfor the tth mode shape, / t, is expressed as (Reynders et al.,2012): The kth component of ~/tis given: ~/t,k¼/t,k/C0RL k¼1/t,k L, Lis the number of components in /t: a¼kImð~/tÞk2/C0kReð~/tÞk2 2Re ~/T t/C16/C17 Imð~/T tÞ(4) c¼arctan ðajjþsignaðÞﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1þa2p Þ (5) MPC values are dimensionless; they lie within the range of 0 and 1. MPC value closer to 1 indicates that modeshape, / t, is more collinear and ‘monophase, ’which is usu- ally regarded as a physical mode. With regard to MPD, it represents the phase degree of each identified mode shape vector. The value of MPD/90lies between 0 and 1. A smaller quantity of MPD impliesthat mode shape vector is more likely to be physical. Adetailed discussion can be found in Reynders et al. ( 2012). For the tth mode shape, / t, the mean phase (MP) is defined as: MP/tðÞ¼arghminkIm/tðÞ/C0tanhRe/tðÞk2 1þtanh/C18/C19 (6) where his phase angle in degree, Eq. (6) can be solved by the least square as: MP/tðÞ¼arctan/C0V12 V22/C18/C19 , USVT¼Re/tðÞIm/tðÞ ½/C138 (7) where USVTis singular value decomposition, V12and V22 denotes elements (1,2) and (2,2) of Vmatrix, respectively. Then, MPD can be determined as: MPD /tðÞ¼RL k¼1xkarccos ðRe/t,k/C0/C1V22/C0Im/t,k/C0/C1V12Þ RN k¼1xk(8) where xkis a weighting factor that equals to the kth com- ponent of the tth mode shape, /t: The selection of threshold values of MPC and MPD depends on measurement conditions and dynamic vibrationproperties. For a structure with clear linear behavior andhigh signal-to-noise ratio, the threshold of MPC and MPDcan be conservatively chosen as 0.7 and 0.3, which impliesthat modes whose MPC are less than 0.7 and MPD exceed 0.3 are regarded as spurious modes. Conversely, the thresh- o l do fM P Ca n dM P Da r ec h o s e na s0 . 3a n d0 . 7i nt h eMPC /tðÞ¼kReð~/tÞk2þ1 aRe ~/T t/C16/C17 Im ~/T t/C16/C17 2a2þ1 ðÞ sin2c/C01/C0/C1 kReð~/tÞk2þkImð~/tÞk2(3)4 J. ZENG AND Y. HOON KIM case of structures with complex behavior (Cabboi et al., 2017 ; Fan et al., 2019 ). Two representative field tests with complex measurement conditions are used to valid-ate the methods. Thus, the values of 0.3 and 0.7 areselected as a threshold for MPC and MPD, respectively,in this study. 3.1.2. Uncertainty criterion Although conventional modal validation criteria remove cer-tain spurious modes, many spurious modes still remain inthe stabilization diagram, slowing down the following pro-cess (herein clustering process). More effective validationcriteria should be adopted to delete as many spurious modesas possible. This study employed supplementary uncertaintycriteria at the pre-processing stage to further eliminatespurious modes. Uncertainty on modal parameters by SSI mainly arise from five sources: (1) finite number of data sample; (2)unmeasured excitation and measurement noise modeled aswhite noise; (3) assumption of linear and stationary behav-ior; (4) imperfect filter of data; (5) incorrect choice of modelorder (Reynders et al., 2008). Reynders et al. ( 2008) initially developed the uncertainty computation based on the propa-gation of first-order perturbation from measured data tomodal parameters. Also, some validation and application aresummarized in Reynders et al. ( 2016). Later, D €ohler, Lam, and Mevel ( 2013) significantly improved the computational efficiency of uncertainty atmultiple model orders, which has been applied in variousstructures (D €ohler et al., 2013). Uncertainty quantification can provide information to measure the accuracy of identi-fied modal parameters. It is the fact that the uncertainty ofphysical modes is smaller than those of spurious modes.Based on this information, coefficient of variation (COV)(standard derivation/mean) with respect to frequency maybe used to distinguish physical modes from spurious modes(D€ohler & Mevel, 2013). Some research has introduced uncertainty features in the stabilization diagram, but uncertainty criterion is not usedor does not contribute to further automated modal proced-ure. General procedures of uncertainty computation aresummarized as follows: /C15Input parameters: the number of block rows in Hankel matrix, q;the amount of data blocks, n b;the range of model order, ðnmin,nmaxÞ; /C15Compute Hankel matrix, H, system state matrix, Aand output matrix, CinEq. (1) , as well as observability matrix, O, based on SSI-cov/ref, then compute transform matrix, T /C15Compute covariance and sensitivity of subspace matrix from SSI-cov/ref, given by ^RHcovandJO,H, respectively /C15Compute sensitivity and covariance of system state matrix, A, and output matrix, Cfrom SSI-cov/ref, given byJA,O,JC,OandRA,C, respectively /C15For each mode iat successive modal order, compute sen- sitivity matrix: Jfi,A,Jni,Aand J/i,A:Finally, compute covariance of modal parameters, frequency, fi, dampingratio, ni;mode shape, /i:covfi ni/C20/C21 ,fj nj/C20/C21/C18/C19 and covReð/iÞ lmð/iÞ/C20/C21 ,Reð/jÞ lmð/jÞ/C20/C21/C18/C19 Comprehensive derivation for uncertainty estimation can be found in D €ohler and Mevel ( 2013). When the COV in the frequency is chosen as a threshold (herein, 2%), modeswith the COV of frequency larger than the threshold will be discarded. 3.2. Clustering stage The clustering stage is sequentially performed to assemble modes based on similarities in modal parameters in this sec- tion. A novel method is proposed to calculate the clusteringthreshold; an improved self-adaptive clustering is then used to identify physical clusters. Finally, robust outlier detection is implemented, and each representative of modal parame-ters is determined. 3.2.1. Automated computation of clustering threshold Typically, two kinds of thresholds are usually adopted for clustering: (1) static threshold; (2) automatically computed threshold. A static threshold relies on the engineers ’judg- ment. Also, during long-term health monitoring, a well-defined static threshold may be suitable for some initial datasets; however, there is no guarantee that the static threshold will be keeping appropriate for all datasets. This ismore challenging in the case of handling massive datasets. In this study, a novel method is proposed to calculate the clustering threshold based on possible physical modes at thepre-processing stage. First, the mutual distance between the two modes is defined as: distance ¼ Fi/C0Fj/C12/C12/C12/C12 maxðFi,FjÞþx1/C0ðUT iUjÞ2 ðUT iUiÞðUT jUjÞ ! (9) where Fi¼fiþ2rfi,Fj¼fjþ2rfj,fiand fjareith and jth identified frequency at a pre-processing stage, respectively; rfiandrfjare corresponding standard derivation, respect- ively;Ui¼/iþ2r/i,Uj¼/jþ2r/j,/iand/jareith and jth mode shapes at a pre-processing stage, respectively; r/i andr/jare corresponding standard derivation, respectively. xis a weighting factor of mode shape differ- ence, x¼r/iþr/j ðÞ 2: Equation (9) does not consider the damping ratio differ- ence because it is difficult to accurately measure the damp- ing ratio in practice. In addition, there is a high probabilityof two different modes with a similar damping ratio. A weighting factor, x,i nEq. (9) represents different participa- tion for frequency difference and mode shape difference.Generally, mode shape is measured with limited sensors, yielding missing components of mode shape; frequency is usually measured with an accurate level. Therefore, the useofxcan reduce the effect of measurement inaccuracy ofSTRUCTURE AND INFRASTRUCTURE ENGINEERING 5 mode shapes on distance calculation (Boroschek & Bilbao, 2019). An uncertainty quantification using standard deriv- ation is used to form a weighting matrix for Finite ElementModel Updating (Yang & Lam, 2018). Similarly, this work adopted the average of the standard derivation of mode shapes to define x:Furthermore, as uncertainty on modal parameters is inevitable in practice, it is more reasonable to incorporate uncertainty in distance calculation. Here, twostandard derivations are considered in Eq. (9) . At the next model order, the mutual distance between one mode and all other modes is computed by Eq. (9) , then the minimum distance is determined. Assuming nmodes have been identified at the pre-processing stage, each mode has its minimum mutual distance with forming a minimum distance vector, V¼ðd 1 min,d2 min,d3 min/C1/C1/C1dn minÞ,(ndenotes the number of modes, dmindenotes the minimum distance between one mode and all other modes). Finally, the sum of mean and two standard derivations of Vare used to com- pute the clustering threshold, d(Reynders et al., 2012): d¼lþ2r (10) Generally, modal features are usually assumed to follow Gaussian normal distribution, such as frequency, damping ratio and mode shape (MAC value) (Au, 2011b ). In this study, Eq. (9) defines modal distance which is the sum of frequency difference and mode shape difference between two modes. Therefore, modal distance turns out to be Gaussian normal distribution. Two standard derivations inEq. (10) guarantee the distance between two modes should be captured within a 95% confidence interval with the assumption of Gaussian distribution. It is worth mentioning that conventional threshold calcu- lation for clustering does not consider the uncertainty of modal parameters and weighting factor. A newly proposed threshold calculation in Eq. (9) incorporates both uncer- tainty and weighted factor, which is considered a contribu- tion of this study to the advancement of automated OMA. 3.2.2. Mode clustering Mode clustering starts with a calculated threshold in Section 3.2.1. to group individual physical modes with similar modal characteristics. This study adopts self-adaptive clustering (Cabboi et al., 2017) to accomplish automated process. But different from original Cabboi et al. ’s work, a weighted dis- tance with an uncertainty of modal parameters is proposed.The ith weighed distance at model order, n, is defined as: d n,i¼~Fz/C0Fn,i/C12/C12/C12/C12 ~Fz ! þc1/C0MAC ð~Uz,Un,iÞ/C0/C1 (11) where ~Fz¼~fzþ2r~fz,~Uz¼~/zþ2r~/z,Fn,i¼fn,iþ2rfn,i, Un,i¼/n,iþ2r/n,i:~fzand ~/zare mean frequency and mean mode shape at the zthcluster, respectively; r~fzand r~/zare corresponding mean standard derivation at the zth cluster. fn,iand/n,iare the ithfrequency, and mode shape at model order, n, respectively; rfn,iand r/n,iare corre- sponding standard derivation, respectively. MAC representsthe modal assurance criteria (Pastor, Binda, & Har /C20carik, 2012).cis a weighting factor to reduce the effect of inaccur- ate mode shape on distance calculation c¼r~/zþr/n,i 2/C16/C17 : Equations (9) and (11) are similar; both consider the uncertainty of parameter estimates and the importance of mode shape difference. Figure 3 shows the flowchart of an improved clustering strategy. A more detailed introduction of self-adaptive clustering is referred to Cabboi et al. ( 2017). 3.2.3. Robust outlier detection The number of physical poles in the stabilization diagram has trends with the increase of model order, exhibiting vari- ability of modal estimates (Neu et al., 2017). The phenom- enon more frequently appears in the damping ratio because the damping ratio has a high scattered nature. Outlier detec- tion is applied to penalize undesirable modes in the final clusters for reducing identification variance from different measurements. In this study, robust outlier detection based on the minimum covariance determinant (MCD) is employed to identify outlying values from physical clusters as described in Hubert, Debruyne, and Rousseeuw (2018). A robust distance (RD) is defined as: RDxðÞ¼dðx,^lMCD,^RMCDÞ (12) where observation sample, x, is either frequency or damp- ing ratio in a physical cluster in our case. ^lMCD is the MCD estimates of location; ^RMCD is the MCD covariance estimate. Explicit derivation and introduction can be found in Hubert et al. ( 2018). A robust MCD estimator based on Eq. (12) is very powerful to flag outliers, as RD in Eq. (12) . It is not sensi- tive to diagnostic tools ’masking effect compared to statis- tical distance and Mahalanobis distance (Cerioli, 2010). Also, MCD has a high resistance to outliers and are more robust and efficient (Hubert et al., 2018). Furthermore, MCD has the advantage of requiring no user-defined threshold, like a box-plot method that needs to define limit bounds (Sarlo & Tarazaga, 2019). Robust outlier detection in this work can be done by the function ‘robustcov ’ in MATLAB. After outlier removal, the average frequency, damping ratio, and mode shape in each physical cluster are taken as a representative. For evaluating the quality of each identified cluster, uncertainty on the zthphysical clusters is quantified by Euclidean norm of uncertainty on modal parameters: rz¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð~r2 f,zþ~r2 /,zþ~r2 n,zÞq (13) where rzis the standard derivation of the zthclusters; ~rf,z, ~r/,z, and ~rn,zare the average values of standard derivations of all frequencies, damping ratios, and mode shapes in the zth clusters. 4. Application In this section, the performance of the proposed approach is validated by two field tests on the bridge, namely, the6 J. ZENG AND Y. HOON KIM Dowling Hall Footbridge located at Turfs University in the USA and the Z24 bridge benchmark located in Switzerland. The data are open sources, and many researchers used these data to test the algorithms in the research community. 4.1. Application 1: Steel frame footbridge Dowling Hall Footbridge is located at Tufts University, as shown in Figure 4(a) . The bridge is a two-span steel frame bridge, 144 ft (44 m) long and 12 ft (2.7 m) wide with a rein- forced concrete deck. Continuous health monitoring was designed and performed on Dowling Hall Footbridge from January 2010 to May 2010. The layout of eight accelerome- ters is shown in Figure 4(b) . More details of Dowling Hall Footbridge ’s information can be found in Moser and Moaveni ( 2011). In this study, the first six modal character- istics are used as baseline results obtained from the litera- ture (Moser & Moaveni, 2011) to evaluate the performance of the proposed approach. The acceleration data used in this study are obtained from vertical measurement under ambient excitation col- lected in the first week at 1:00 p.m. on January 7, 2010. The frequency range of interest is 0 –30 Hz. The sampling fre- quency is 128 Hz. Preparation parameters for SSI-cov/ref in this application are: i¼60, model order n¼40–150, refer- ence sensor ¼ð1, 2, 3, 4, 5, 6, 7, 8 Þ:4.1.1. Identification results The proposed approach described in Section 3 is utilized to analyze measured data. Figure 5(a –c)displays the modal iden- tification results at the pre-processing stage. The singularvalue spectrum (appeared in the curves in Figure 5 ) is plotted below the stabilization diagram. The standard derivation (6r) uncertainty bounds of the frequency are shown as hori- zontal bars. Figure 5(a) displays all possible physical modes remaining in the stabilization diagram after applying conven- tional validation criteria, e.g., damping ratio check and modal complexity check. Figure 5(b) depicts the stabilization dia- gram filtered by a supplementary uncertainty criterion. It isobserved only using conventional validation criteria, the sta-bilization diagram still looks busy, including lots of scattered poles, which are spurious modes. However, uncertainty criter- ion can eliminate as many spurious modes as possible com-pared to conventional validation criteria, which will speed uplater automated processes. The pre-processing stage ’s identifi- cation results demonstrate that the uncertainty criterion is more effective than conventional validation criteria. The clustering stage then starts with a calculated cluster- ing threshold using Eqs. (9) and(10) based on the remain- ing modes in Figure 5(b) . The proposed method ’s calculated threshold in this example is 0.022, while without the weight- ing factor, x, it is 0.0488. It implies the stricter threshold by Eqs. (9) and (10) that allows removing more spurious modes with keeping physical modes. Furthermore, the Figure 3. The flowchart of an improved self-adaptive clustering. Figure 4. Description of Application 1: (a) Dowling Hall Footbridge; (b) Sensor layout (Moser & Moaveni, 2011 ).STRUCTURE AND INFRASTRUCTURE ENGINEERING 7 updated threshold by improved clustering with Eq. (11) is 0.0086. Still, the original work (distance calculation without weighting factor) gives the updated threshold as 0.0304, indicating that the weighted distance tends to give a smaller value of the updated threshold.The identified modes are more consistent and stable. It may be attributed to the use of xcan improve the accuracy of measured mode shapes by distance calculation. Figure 5(c) shows the modal identification results after performing the improved self-adaptive clustering. It is observed that Figure 5. Identification results: (a) after conventional validation criteria; (b) after uncertainty criterion (c) after improved self-adaptive clusterin g. Figure 6. Damping ratio vs frequency: (a) before outlier detection; (b) after outlier detection.8 J. ZENG AND Y. HOON KIM clustering procedures remove spurious modes, the stabiliza- tion diagram is clarified with only remaining stable modes (vertical alignments). The first six modes in the reported work (Moaveni & Behmanesh, 2012) are used as a baseline for comparison, marked as M 1to M 6(a total of six clusters) inFigure 6 . It is noted that M 5and M 6are closely spaced modes, which are a common challenge in the OMA. The robust outlier detection is to remove outlying frequencies and damping ratios. As seen in Figure 6 (b) , damping ratios are tighter and more consistent. Table 1 presents identified frequencies and damping ratios along with the baseline data. Identified frequencies in this work agree well with those in the literature; the maximum relative difference (2.05%) is observed in the third mode. While larger variation is found in terms of damping ratio. It is because two tests were performed at a different time. Moaveniand Behmanesh ( 2012) reported baseline data, measured on April 4, 2009. In this study, measured data were collected on January 7, 2010. When considered the effect of environmental variables such as temperature, it is not surprising to have these differences. The frequency is less sensitive to environ- mental effects than the damping ratio. The uncertainty on modal parameters and physical clus- ters are also investigated in this example. The frequency and damping ratio in each mode are plotted as an open circleoverlapping the standard derivation ( r) error bar in Figure 7. And the uncertainty of modal frequencies is much smaller than those of damping ratios. It is often more difficult to accurately measure the damping ratio in practice. The uncertainty of identified physical clusters is also quantified byEq. (13) and shown in Table 2 . And the uncertainty of the sixth cluster is much larger than those of others,suggesting it is more challenging to identify the sixth cluster because this cluster contains weakly-excited and closelyspaced modes. Overall, the proposed approach successfully identifies six modes under ambient vibration, as shown in Figure 8 . The first six global mode shapes with corresponding uncertain-ties are presented; 62runcertainty bounds are plotted as blue dashed lines. Identified mode shapes have good agree-ment with those identified in the reported work (Moaveni &Behmanesh, 2012). Modes 3 and 4 are bending-torsional modes with evident rotational motion, while only verticaldeformation is found on other modes. In addition, uncer- tainty bounds for all modes are narrow, which concludes that the identification of mode shapes is accurate. For continuous SHM, it is crucial to track the change of modal parameters over time. In this example, the proposedapproach is applied to modal tracking with measured datacollected at every 1:00 p.m. from January 5 to February 28,in 2010 (total 55 datasets). The same procedures as the for-mer data analysis are applied for modal tracking. As showninFigure 9 , solid black lines indicate frequency estimates, and grey areas cover 62 standard derivations. All six modes are identified and tracked for all datasets by the proposedapproach. It is not surprising that frequencies varied overtime, mainly because of environmental change and ambientexcitation ’s randomness. The frequency at mode 6 has a relatively larger variation for two months, as this mode is not excited well and unstable to environmental change. Theresults illustrate the proposed approach can analyze massivedatasets with minimum human intervention. 4.1.2. Sensitivity analysis Two preparation parameters in SSI, e.g., maximum model order, n max, and time lag, i, significantly affect identification results as described in Section 2. The influence of nmaxandi is investigated to demonstrate the proposed approach isTable 1. Identification results. ModesFrequency (Hz) Damping ratio (%) BaselineaThe proposed approach BaselineaThe proposed approach 1st (M 1) 4.63 4.63 1.0 0.8 2nd (M 2) 6.07 6.04 0.6 0.1 3rd (M 3) 7.07 7.21 0.7 0.6 4th (M 4) 8.90 8.95 0.3 0.1 5th (M 5) 13.13 13.24 0.8 0.1 6th (M 6) 13.56 13.46 1.1 1.1 aNote : Moaveni and Behmanesh ( 2012 ). Figure 7. error bar of frequency (left, 62r) and damping ratio (right, 6r).Table 2. Uncertainty of physical clusters. No. cluster 123456 STD (%) 0.184 0.001 0.430 0.010 0.119 2.712 Note : STD denotes standard derivationSTRUCTURE AND INFRASTRUCTURE ENGINEERING 9 Figure 8. Identified mode shapes by the proposed approach and 62runcertainty bounds (blue dashed lines). Figure 9. Identified frequency of modes 1 –6 with 2-month data. Black solid lines: frequency estimates; grey shaded areas: 6two standard derivations.10 J. ZENG AND Y. HOON KIM robust to their choice. nmaxand iare varied from 70 to 160 and 30 to 120 in intervals of 10, respectively. As shown in Figure 10 , identified frequencies are almost invariant to a different choice of nmaxand i, suggesting the proposed approach is robust and not sensitive to these two prepar- ation parameters. It is very difficult to identify the best nmax and iin practice (Neu et al., 2017; Ubertini et al., 2013). Thus, insensitivity to them allows to more conveniently per-form automated OMA and continuous health monitoring. On the other hand, a different choice of uncertainty threshold (COV of frequency) is utilized to evaluate itseffect on identification results. The uncertainty threshold is varied from 1% to 5% in the interval of 1%. As shown in Figure 11 , the proposed approach yields the same frequen- cies regardless of COV thresholds, indicating a COVthreshold can be safely chosen in the range of 1% /C20COV /C205%: 4.2. Application 2: Concrete box girder bridge The proposed approach is also applied to the Z24 bridge benchmark to validate its performance. The Z24 bridgewas built in 1963 and located in Switzerland, serving to connect Koppigen with Utzenstorf and crossing over the A1 highway (see Figure 12(a) ). It is a post-tensioned con- crete box-girder bridge with a main span of 100 ft (30 m)and two sides span of 46 ft (14 m). Detailed introduction of the Z24 bridge can be found in Maeck and De Roeck (2003). The Z24 bridge was demolished at the end of 1998. Before the complete demolition, a short-term pro-gressive damage test was implemented on the bridge to investigate the effect of simulated damage on the safety of the bridge. A total of 17 different damage scenarios were designed under full forced and ambient excitation (Reynders & Roeck, 2009). In this work, acceleration response data from the scenario of No.8 for the new reference condition underambient excitation is used to assess the proposed approach. A total of 291 DOFs were measured (see Figure 12(b) ). Dueto the limited number of sensors, only at most 33 DOFs were measured for each set-up. Therefore, nine measure- ment set-ups were recorded with most 33 sensors to have full location coverage of the whole bridge, containing fivereference sensors that are common to each set-up and 28moving sensors whose location changes with different set- ups. For the No. 5 set-up, only 22 moving sensors were used. Samples of 65,536 data were recorded at each set-upat a 100 Hz sampling rate. For each dataset, preparation parameters in SSI-cov/ref are defined as: time lag is i¼50, model order ranges from 2 to 120, to create stabilization diagrams. Reference sensorsare selected as No. 29-33 (for set-up No. 5, as No. 23-27). After the stabilization diagram is created, the proposed approach in Section 3 is applied to automatically interpretthe stabilization diagram. 4.2.1. Identification results Nine stabilization diagrams corresponding to each set-up are created; results of the fifth set-up are only presented inFigure 13 due to space limitation in this article. The singular Figure 10. Frequencies at different parameters: (a) model order range sensitivity (fixed i¼60); (b) time lag range sensitivity (fixed nmax¼100). Figure 11. Frequencies at different COV threshold: ( i¼60;nmax¼100).STRUCTURE AND INFRASTRUCTURE ENGINEERING 11 value spectrum (appeared in curves in Figure 13 ) is also plotted below the stabilization diagram. 6r(standard deriv- ation) uncertainty bounds of frequency are plotted as hori-zontal bars. Figure 13(a) displays modal identification results using conventional validation criteria, many scatteredpoles which are definitely spurious modes, still retain in the stabilization diagram. However, the uncertainty criterion can remove most spurious modes, demonstrating that theuncertainty criterion is more effective than conventional val- idation criteria (see Figure 13(b) ). Figure 12. Description of the Z24 bridge: (a) Front and Top view (Maeck & De Roeck, 2003 ); (b) sensor layout (D €ohler et al., 2013 ). Figure 13. Pre-processing stage for No. 5 set-up: (a) after validation criteria; (b) after uncertainty criterion. Note: xis weighting factor in Eq. (9) . Figure 14. The clustering stage: (a) calculated clustering threshold; (b) after improved self-adaptive clustering for No. 5 set-up. Note : c is the weighting factor in Eq. (11)12 J. ZENG AND Y. HOON KIM Based on the remaining poles in the stabilization diagram after the pre-processing stage, the clustering threshold foreach measurement set-up is calculated using Eqs. (9) and (10). As shown in Figure 14(a) , all threshold values are sig- nificantly reduced compared to those calculated withoutweighting factor, as weighting factor can offset the effect of mode shape difference. Manual clustering thresholds in commercial OMA software are usually below 0.06 (Neu etal.,2017). The threshold derived from the newly proposed method is closer to the one from manual analysis, indicating proposed method ’s rationality and feasibility in practice. Mode clustering described in Section 3.2.2 is then imple- mented to group physical modes. The number of scattered poles is greatly removed by the proposed approach, onlyremaining obvious vertical alignments in the stabilization diagram in Figure 14(b) . In addition, an improved self-adaptive clustering that considers the weighting factor of cinEq. (11) tends to give a smaller updated threshold, implying identified modal parameters are more stable and consistent with each other(seeFigure 15 ). The use of the weighting factor can improve the performance of clustering. Because only the first six modes are present in baseline work (Reynders et al., 2012), the first six clusters are presented in Figure 14(b) , marked as P 1–P6:Robust outlier detection is used to identify outly- ing modes (See Figure 16 ). Finally, the average of modal parameters in each cluster is selected as representative. Table 3 shows the sample mean and sample standard derivation of frequency and damping ratio over nine meas-urement set-ups obtained from the proposed approach. Thestandard deviation in Table 3 represents the setup-to-setup sample statistics among all set-ups. The calculation of sam- ple standard derivation in Table 3 only considers the envir- onmental change among different set-ups rather than uncertainty sources summarized in Section 3.1.2. It is seen from Table 3 that the damping ratio has more significant variability than frequency, implying it is more challenging to identify damping ratio in practice, as the damping ratio is sensitive to environmental change. Overall, the proposedapproach ’s identified frequencies and damping ratios arealmost identical to those from baseline work, demonstrating low demand for human intervention. Uncertainties of modal parameters airing from assump- tions made in SSI, such as linear, stationary structuralbehavior, white noise, etc., are also studied. Figure 17 depicts the variability of frequency and damping ratio frommodes 1 to 6 across nine measurement set-ups, respectively, with open circles representing the parameter estimates and error bars covering 62rstandard derivations. Both frequen- cies and damping ratios change over time, while the damp-ing ratios have larger uncertainties. The negative dampingratio is immaterial in Figure 17(b) , such as mode 5 at No. 4 set-up and mode 6 at No. 1 set-up, merely because of theGaussian distribution approximation and the larger standardderivation. Mode 6 has relatively larger uncertainty sincethe mode is not excited well. Table 4 presents the average of standard derivation for each cluster over nine measurement set-ups using Eq. (13) . As expected, the sixth cluster has the highest uncertainty,implying it is relatively harder to identify this cluster, whichis also reflected in Figure 14(b) that the sixth vertical align- ments from the left form at a very ambiguous peak.Generally speaking, quantities of identified frequency anddamping ratio are consistent from one to another set-upnumbers, suggesting robust and fair performance on modal analysis. The global mode shapes are directly assembled from a local one in a single dataset by multiplying by a scaling fac-tor so that their common DOFs, at the location of