The use of recursive subspace-based identification methods is analyzed in the estimation of the most significant vibration frequencies along an elevated railroad segment of México City Metro Line 12 using ambient vibration measurements recorded from 2012, when the line was opened, to 2018. Due to the railroad characteristics, the use of high-order models and the systematic tuning of the methods are required to achieve low uncertainty in frequency estimation. A frequency history is generated using these high-order models in order to check for variations along the seven years where important events took place: two low- and one high-intensity earthquakes, paving, and construction of sidewalks and planters around the sensor station. Results are consistent for all methods under analysis in the identified frequencies, suggesting that the system has preserved its structural health. To produce independent results, spectral analysis was performed showing that the associated frequency history is again consistent with that generated with recursive subspace-based identification methods. Overall, results indicate that these subspace methods are suitable for frequency monitoring in the studied system offering, in the case of recursive N4SID, important advantages in terms of low computational cost, real-time implementation, and smaller uncertainty.

1. Introduction

Parameter identification in structural models is an important topic in civil engineering applications such as structural health monitoring, damage detection, and vibration control systems. Appropriate model identification allows the supervision of the dynamic properties of a structural system and to observe their variation in the presence of external perturbations, which may suggest structural damage. Related literature about identification methods is abundant. Ambient noise information or structure earthquake response is used in some methods to get modal properties [13]. Recently, identification methods based on subspaces, as those reported in [46], have been used in determining and monitoring modal properties in both buildings [7, 8] and bridges [9, 10]. A review and comparison with classical identification methods can be found in [11, 12]. Some critical aspects to be considered when applying these methods are their correct tuning, pursuing a balance between the identification accuracy and the demand for computational processing, and the presence of spurious modes that are mainly caused by noisy signals, coupling effects, and nonlinearities presented in the system [13, 14].

The goal of this work is to evaluate the effectiveness of two recursive subspace identification methods: numerical algorithms for subspace state space system identification (N4SID) and multivariable output error state space (MOESP) in a structural health application. Both methods are used to monitor the two most significant radial frequencies in a curved segment of an elevated railway in México City Metro Line 12, opened in 2012, where all the mentioned critical aspects are involved. To achieve this task, ambient vibration records from 2012 to 2018 are studied using the subspace methods previously mentioned and spectral analysis. Recursive formulations of both methods are investigated to assess their incorporation into a warning system [15] that, based on various structural performance indicators, including the frequency variation, provides information about the state of the structure before, during, and after seismic events take place. In the studied system, analysis is limited to the use of ambient vibration records.

This paper differs from other papers in the literature of structural health monitoring in several points [112] and references therein. In the first place, the nonlinear and coupling effects of the railroad segment require the use of high-order models, and if low uncertainty is desired, that has a direct impact on the complexity of the calculations and therefore in the feasibility of real-time implementation of structural health monitoring systems. In the second place, the online use of ambient vibration records implies a challenge in the methods tuning and calibration. For this purpose, extensive analyses with a long series of ambient noise records were performed to find systematic procedures to tune the online identification methods in order to reduce their convergence time and to minimize their uncertainty. To the authors’ knowledge, there is no similar analysis available in the literature of the two subspace identification methods under investigation for their online implementation using a limited number of measurements. Finally, a comparison with well-known spectral analysis tools is also included that increases the confidence of the obtained results.

Content is organized as follows: Section 2 describes the curved segment of the railway studied and the available instrumentation. In Section 3, recursive subspace methods used in the identification process are described. Section 4 explains the spectral and subspace methods used and their tuning. Section 5 shows and analyses the results, while Section 6 presents the conclusion and directions for future work.

2. Case Study and Mathematical Model

The subspace identification methods reported in [5, 6] were developed for the identification of linear systems. These methods construct a discrete state space mathematical model of the following form:where the system’s state corresponds to displacement and velocity signals evaluated at time . On the contrary, and represent the system input, the soil acceleration, the output, and the acceleration measured at specific points, respectively, evaluated at time . A structure with degrees of freedom, inputs, and outputs has states, and the state space matrix dimensions are , , , and .

Located in a zone of soft soil, in an area of high seismic risk, the studied system corresponds to a curved segment of a railway in México City Metro Line 12. The structural system is composed of several precast prestressed concrete elements, which were assembled on-site, and consists of one girder supported by two cantilever columns separated by 30 m with a height of 9.25 m connected to four 37 m piles of 1.2 m diameter. In this segment of the viaduct, the support of the girder consists of two pot bearings at the ends, which allow movements of rotation and horizontal translation. Figure 1 shows the schematic diagram of the structure and location of the sensors in both tangential (L) and radial (T) directions.

There is previous work in this railway segment, based on offline data processing. In [16], static and dynamic tests are reported, while in [17, 18], the effect of some large-scale earthquakes over the system properties was analyzed. In [15], a system was developed to monitor the health of the structure based on the analysis of short time events and long-term effects. The first conclusion of these studies was that the structure health has been preserved, as relative displacements between girders and columns, tangential and radial rotations at both foundation and capital levels, and vibration frequencies in the radial direction of the curve did not report significant differences from 2012 to 2017. The second conclusion was about the importance of monitoring frequencies in the radial T direction as when a train passes at high speed, larger displacements were registered in the T direction when compared with those in the L direction. Based on these conclusions, in this work, frequency monitoring is performed only in the radial T direction where acceleration signals used correspond to A5CT–A6CT (columns 5 and 6) and AET (soil) (Figure 1).

In this railroad system, there are difficulties in identifying frequencies that arise from two sources: the characteristic of the measured signals and the system itself. Regarding the first point, the measured signals, the pass of trains, and the maintenance work were occasionally registered at the top of the column, while traffic was registered at the bottom. In addition, the input sensor is located over soft soil where the dominant frequencies range from 0.9 to 1.1 Hz. On the contrary, based on field tests performed in 2014 [16], in this curved segment of the railway, nonlinear behavior and coupling effects have been observed, and both tangential and radial frequencies lie in similar intervals 2.1–2.3 Hz and 2.0–2.2 Hz, phenomena that can be related to the curved system shape. Moreover, the identification also can be affected by the presence of nonstructural elements.

3. Identification Methods

In this section, recursive subspace methods used in the identification process are described. The first step of the N4SID method is to sort the input and output signals to build the following Hankel matrices:where is the number of block rows used in its construction and is the number of data used in the identification process. Then, a relation between the input and output signals is proposed:where matrix projects the row space from one matrix over the complement row space of the other. After this calculation, the observability matrix is recovered through a singular value decomposition (SVD) of matrix , allowing to extract the pair (, ) and from here the modal information of the system. This method uses information windows from a finite set of data and provides only one identified model. In [5], a systematic and recursive process is shown where the observability matrix is updated with each new input-output pair data. Here, is the input/output compressed matrix, whose dimension is constant. Additionally, corresponds to the forgetting factor with . Then, with an adequate tuning and by updating matrix , the system matrixes are identified recursively. With a slight difference in the treatment of the data, the objective of both MOESP and N4SID methods is to recursively update matrix . For the description of MOESP, readers can refer to [4, 5]. In both methods, a singular value decomposition is performed each iteration over matrix :such that the system order can be determined by removing the elements with the lowest values in the diagonal matrix . Stabilization diagrams are a useful tool to choose the identified order model; consistency of modal properties is analyzed as a function of the identified order model. So, matrix has the singular values associated with . Finally, the column space corresponds to the image space of the extended observability matrix , which can be derived from . The system matrix is calculated as follows:where the above matrix product is performed by removing the last rows of the first term and the first rows of the second term. Superscript † represents the Moore–Penrose pseudoinverse. Frequency and damping ratios are derived from the eigenvalues in matrix :where is the -eigenvalue in and represents the sampling time. Mode shapes are derived from the eigenvectors in matrix and the observability matrix . In the studied system, the use of mode shapes is limited because these are observed in a single point, the upper part of the columns.

Differences between N4SID and MOESP subspace identification methods are described, for example, in Kim and Lynch [19] who showed that the N4SID method uses an oblique projection, whereas the MOESP method uses an orthogonal one. The advantage of using orthogonal projections is in the capability to deal with ill-conditioned Hankel matrices produced by poor quality input/output signals obtained from experimental tests performed in noisy environments, as it is mentioned, for instance, by Hachicha et al. [20]. However, the MOESP method demands a high computational effort, the orthogonal projection involves a LQ decomposition whose computational processing increases exponentially as the order of the identified model increases. Studies as those of Chen and Loh [21] focus on improving the MOESP method by reducing its demand for computational resources. The use of improved subspace identification methods, most of them derived from the original ones, was left for future work because the goal in this paper is to evaluate its use as monitoring systems and later integration into an alert system for the studied system.

The first phase of this evaluation consists in tuning both subspace methods, and this was done by means of a 1200 s ambient vibration record from 2013, with a sampling time of 5 ms, whose acceleration signals and identified frequencies, based on average power spectrum ratios with 40 s time windows, are shown in Figure 2. For that purpose, using stabilization diagrams allows checking for the minimum order where identified frequencies are consistent. Using this ambient vibration record from 2013, Figure 3 shows the first and second identified system frequencies with nonrecursive N4SID and both recursive N4SID and MOESP methods for models from order 20 to 100 with increments of two, where the dotted line refers to the spectral values reported in Figure 2. As mentioned, the nonrecursive N4SID method provides only one identified model and recursive methods provide a frequency history at 0.5 s intervals. The frequency shown in Figure 3 corresponds to the average frequency reported in the recursive subspace identification methods. It can be observed that the identified frequencies are consistent for a model order of 50 and higher; lower order models report inconsistent results in both nonrecursive and recursive methods. This finding shows that the task of identifying the principal system frequencies cannot be underestimated and it is not a trivial task for the system under investigation.

The second step was related to determining Hankel matrix dimensions. The literature suggests that the number of rows to be higher than the order model , meanwhile the data number should be higher than . On the contrary, the MOESP method can provide accurate results with less information but demanding high computational effort. Considering this, Hankel matrices were built as follows:(i)N4SID: the number of rows was fixed to and the data number to (60 s). Condensed data matrix has a dimension of . System frequencies are identified every 0.5 s.(ii)MOESP: the number of rows was fixed to and the data number to (30 s). Condensed data matrix has a dimension of . System frequencies are identified every 0.5 s.

As an extra step, spurious modes were removed. The presence of such modes occurs either by sensor problems as bias and offset values or by external perturbations as the passage of trains in the elevated viaduct, as well as cars and trucks close to the sensor station. Marchesiello et al. [14] proposed a criterion to remove spurious modes by using identified mode shapes and the modal mass. In the studied system, the elevated railway section, mode shapes are observed at one point, the upper part of the railway, limiting the use of that and similar criteria. Additionally, identified mode shapes may result in complex values that can be associated with nonlinear behavior, coupling effects, and low signal-to-noise ratio signals [22, 23]. Spurious modes were removed by inspecting eigenvalues. First, pairs of complex conjugated eigenvalues were preserved removing the nonpaired ones, and later, eigenvalues with higher damping values (20%) were removed. Although limited, these steps allow removing most of the spurious modes caused by sensor errors or by external perturbations.

Figure 4 shows the identified system frequencies using the acceleration records of AET and A5CT sensors, where it is clear that the MOESP identification method provides results with lower variations although with a processing time of 200 min, compared with less than 60 s used by N4SID. The simulations were performed in MATLAB software in a workstation Windows 10 OS, 3.06 GHz CPU and 16 GB RAM. Apparently, the advantage of using the MOESP method seems to be minimal; the variation in the identified frequencies is less than 1% when recursive N4SID is used.

4. Spectral Analysis

Spectral analysis represents a comparison point in regard to the results based on subspace identification methods, and in this section, spectral analysis applied to 300 s records is discussed. First, using the input and output average Fourier spectra, in intervals of 10 s, spectral ratios are determined for both columns with respect to the ground acceleration signal: A5CT/AET and A6CT/AET. Then, these ratios are analyzed around the frequency ranges of interest. Considering as an example, for the first frequency of the system, the spectral analysis is summarized as follows:(1)Magnitudes of spectral ratios, whose frequencies lie within 2 to 3 Hz range, are normalized with respect to the highest value. This provides spectral ratios with values smaller than or equal to one.(2)Normalized spectral ratios and their corresponding frequency values are sorted at user-defined intervals, for instance, from 0.75 to 1. Then, the most representative frequencies of the system are registered.(3)Average frequency and standard deviation values are calculated. Then, frequency values outside the average and standard deviation range are removed and these statistical variables are calculated again.

This process is executed again inspecting the spectral ratios whose frequencies lie in the interval from 5 to 6 Hz, where the second system’s frequency is located. This spectral analysis was performed for each one of the ambient vibration records from 2012 to 2018 in both columns. The above process provides an estimate of the first and second most significant frequencies of the system in a fast way; each identification demands less than 0.1 s of computational processing time. The goal is to obtain a baseline to compare the use of both recursive subspace identification methods. In order to illustrate this procedure, spectral identification was performed considering the same test signal of 2013. In Table 1, average frequencies and standard deviation are reported and compared with those ones identified by using N4SID and MOESP subspace identification methods. It can be observed that the standard deviation is higher when using spectral analysis, but it is lower than 5% of the average frequency values. Comparing with respect to the values obtained by means of power average power spectrum ratios, the difference is below 5% using both spectral analysis and subspace identification methods.

Once tuned both recursive subspace identification methods and having performed the spectral analysis, the next step was the identification of system frequencies by using ambient vibration records from 2012, when the Metro Line was opened, to 2018. The identification process was limited to the use of 300 s ambient vibration records with 200 samples per second, where frequency and damping values were identified at 0.5 s intervals when subspace identification methods were used. In the next section, results are presented when the first and second system frequencies are identified using both subspace identification methods and the spectral analysis previously mentioned. The results are organized by column, each one of them located at the extreme of the curved railway section.

5. Identification Results

While ambient vibration records from 2012 to 2014 were registered on specific dates, ambient vibration tests were carried out on a scheduled basis every day from May 2015 to February 2018; except for December 2015 and from July to August 2016. Measurements were obtained always at 2:00 a.m. because the Metro Line 12 is out of use at this time and there is low variation of ambient temperature, although the passage of trains and maintenance work were registered very occasionally. The mass of the maintenance trains does not affect significantly the frequency of the structure because it is much smaller than that of service trains. About the temperature, registered from 2012 to 2018 at 2:00 a.m., the average and standard deviation values are 14.3 ± 2.5 and 16 ± 3°C in columns 5 and 6, respectively. At night, in México City, there is low variation in the temperature during the seasons. So, a correlation was not observed between the frequencies and the temperature. After removing measurements either ill-recorded or registered at different times and trying to preserve consistency and homogeneity in the monitoring of system frequencies, the number of records by column was 728. It is important to mention that due to the high demand for computer processing, only 160 ambient vibration records by column were selected with a symmetrical separation in the MOESP identification case. Figure 5 shows the average identified frequencies from 2015 to 2018; global mean and standard deviation values are marked. It is to note that results are consistent using both subspace identification methods being slightly higher the variation in MOESP case, but below 5%.

To analyze the results more efficiently, these are first summarized by month from 2012 to 2018 where spectral analysis and N4SID identification results are compared. After this, results are organized considering a set of events that took place, including one high-intensity earthquake in September 2017 and two low-intensity earthquakes one in September 2017 and another in February 2018. Results based on the MOESP identification method are introduced at this step where a set of representative ambient vibration records were considered before and after the reported events.

5.1. Frequency History Analysis

Starting with the first column (A5CT), Figure 6 shows the average identified frequencies by month from 2012 to 2018. Standard deviation values were calculated based on the average identified frequencies per day throughout each month. These standard deviation values are below 3% of the average frequency ones, using both N4SID and spectral analysis. Although there is a difference in tracking the frequencies by using N4SID or spectral analysis, this difference remains within confidence thresholds and the path followed by the frequencies is comparable, as it can be observed in Figure 6. Both system frequencies present a significant increase from August 2012 to May 2015, the first one from 2.16 to 2.28 Hz (5.56%) and the second one from 5.31 to 5.65 Hz (6.40%). Subsequent variations in monitored frequencies are minimal, until September 2017 when a high-intensity earthquake of 7.1 Mw took place at 18:14 UTC with epicenter at 110 km south of the site of the instrumented elevated railroad segment of Metro Line 12 in México City; previous earthquake, of September 7th, did not affect the system frequencies. The most affected system frequency was the first one, reporting a reduction from 2.32 to 2.25 Hz (3.03%) to subsequently recover gradually. Later, a 7.2 Mw earthquake occurred on February 16, 2018, at 23:39 UTC in Oaxaca State, where no significant changes in system frequencies were observed.

Similar results were drawn in the second column (A6CT), as it is shown in Figure 7 from August 2012 to May 2015. The first frequency reveals an increase from 2.17 to 2.29 Hz (5.80%) and the second one from 5.37 to 5.61 Hz (4.42%). In the second column case, there are two aspects to be mentioned: first, the difference by using N4SID or spectral analysis is notoriously lower and, second, the effect of the September 2017 earthquake on the first system frequency is slightly more noticeable, and it decreases from 2.33 to 2.26 Hz (3.01%), although once again the preevent frequency value is gradually recovered. Finally, except for the increase reported from August 2012 to May 2015, the second system frequency did not report significant variations in any column, even when the mentioned earthquakes took place. Therefore, both columns reported consistent and comparable results suggesting N4SID and spectral analysis are suitable for monitoring the most significant vibration frequencies.

In order to analyze more easily the system frequency variations, a set of important events that took place between the years 2012 to 2018 are shown in Table 2. These events could affect the system frequencies. The first three events occurred due to human intervention: paving, sidewalks, and planters constructed around the field station and maintenance work consisting of ballast change in the railways of the curved segment. The last three events correspond to the previously mentioned earthquakes. Based on these six events, results were once again organized and summarized by sections in Table 3 and Figure 8 for column 5. The most notorious variations occurred in the first system frequency. Following the N4SID identification method case, the first system frequency starts at 2.18 Hz and shows an increase to 2.30 Hz after event . Next, there is a reduction to 2.22 Hz and subsequent recovery to 2.30 Hz, before and after events and , respectively. Then, the frequency value drops to 2.25 Hz when event , September 2017 earthquake, took place. Finally, the frequency is recovered apparently and unaffected by event , the earthquake occurred in February 2018. On the contrary, the second system frequency was only affected significantly by events and with an increase from 5.43 to 5.62 Hz, and later, the frequency remains around 5.60 Hz. In the MOESP identification case, the first system frequency increases from 2.08 to 2.27 Hz after event and recovers to 2.21 Hz after . Next, change is observed before and after event with a reduction from 2.30 to 2.22 Hz and then shows a gradual frequency recovering. In the second system frequency case, MOESP is consistent with N4SID results. Similar behavior can be observed when spectral analysis is used. Comparable and consistent results were reported for column 6, as it is shown in Table 4 and Figure 9.

5.2. Damping History Analysis

By using the 2013 ambient vibration record, identified damping ratios were as follows: versus () and versus () by using N4SID and MOESP identification methods, respectively. Figure 10 shows the average identified damping ratios and corresponding standard deviations from 2015 to 2018. N4SID standard deviation values are 41 and 26% for damping ratios and , respectively; similar variation results have been reported in [10, 13] when using subspace identification methods. In contrast, MOESP standard deviation values are 53 and 87% for damping ratios and , respectively. Results based on the N4SID method remain in a more stable region such that they can be considered as a comparison value when using other identification methods.

The difficulty in identifying damping accurately is well known. For example, Papagiannopoulos and Beskos [24] presented a review comparing damping identification methods, when buildings are subjected to earthquakes, based on wavelets, modal minimization, least squares, and frequency-domain transfer function. An important remark is about the complexity in identifying accurately damping when structural systems show either nonlinear behavior or interaction soil-structure. In the studied system, elements that affect the accuracy in the identification are the interaction and coupling modes, noise-signal ratio, and diverse sources of nonlinearity, such as the support conditions at the beam-column interface, rail-structure interaction, and soil-structure interaction effects. Based on these observations, the identification of damping for the studied system should be carried out using different identification methods and excitation signals to establish confidence intervals of damping.

6. Conclusions

Based on ambient vibration records from 2012 to 2018, the two most significant vibration frequencies and corresponding damping coefficients were monitored along two columns of an elevated railway portion in México City Metro Line 12. This was achieved by the use of two recursive subspace identification methods and spectral analysis. This system has nonlinear behavior and coupling effects that require the use of high-order models and a systematic tuning of the subspace identification methods, in order to reduce uncertainty. The use of high-order models increases the dimension of the involved Hankel matrices and the corresponding computational cost.

Results show consistency in tracking both system frequencies where variations caused by external events such as paving and sidewalks construction and the presence of strong earthquakes can be detected. From 2012 to the first third of 2014, both frequencies reported an increase below 10% due to the paving and sidewalks construction around the field station. After the maintenance work in mid-2015, ballast change in the railways, both frequencies reported a minimal decrease and remained constant until September 2017 when an earthquake of high intensity took place and where only the first frequency was affected with a reduction around of 5% with subsequent gradual recovery. Finally, the earthquake reported in February 2018 did not affect the system frequencies. On the contrary, even though the damping coefficients identified are not conclusive about their real values, these were consistent and remained constant from 2012 to 2018.

These results were consistent in both columns using both subspace identification methods and confirmed with the spectral analysis. It can be concluded that recursive subspace identification methods N4SID and MOESP are suitable as structural health monitoring systems. Moreover, even with high-order models, the N4SID method can be used online because of its low computational processing cost, reserving the MOESP method to validate the identification results. Part of the future work consists of using the recursive N4SID method identifying modal parameters and its variation when the pass of trains and seismic events are registered.

Data Availability

The use of the data in the paper requires permission from Universidad Nacional Autónoma de México.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The support of Mexico City’s government and the active involvement of the Institute of Engineering of the National Autonomous University of Mexico (IIUNAM) to design, deploy, and maintain the monitoring system of the railway is recognized. We appreciate the help of staff in charge of the operation and maintenance of the instrumentation of the structure presented. The participation of the following technicians and students is gratefully appreciated: Luis A. Aguilar, Mauricio Ayala, Carlos Huerta, Miguel A. Mendoza, Alejandro Mora, José Camargo, and Leonardo Ramirez.