Operational modal analysis for estimating the dynamic properties of a stadium structure during a football game

During a football game, the ambient vibrations at the roof of a football stadium were recorded. A very large data set consisting of 4 hours of data, sampled at 80 Hz, is available. By a data reduction procedure, the complete data set could be analysed at once in a very short time. The data set was also split in shorter segments corresponding to certain events before, during and after the game to investigate the influence of varying operational conditions on the dynamic properties. As the structural vibrations were caused by unmeasurable wind and crowd excitation, Operational Modal Analysis has to be applied to find the dynamic characteristics of the structure. The new operational PolyMAX parameter estimation method is used and compared with Stochastic Subspace Identification. Stochastic Subspace Identification requires the correlations between the responses as primary data, whereas PolyMAX operates on spectra or half spectra (i.e. the Fourier transform of the positive time lags of the correlation functions). The main advantage of PolyMAX is that it yields extremely clear stabilisation diagrams, making an automation of the parameter identification process rather straightforward. This enables a continuous monitoring of the dynamic properties of a structure.


Introduction
For more than a year, the University of Sheffield monitored the Midland Road Stand at the Bradford and Bingley Stadium, home of Bradford City Football Club (Fig. 1, Left).Twelve accelerometers were used to instrument this stand.They were located at 6 frame locations at the intersection between the column at the back of the seating deck and the cantilever roof beam (Fig. 1, Right).Additionally, two video cameras were used to monitor the crowd for correlating its distribution and behaviour with the measured accelerations [1,2].The key aspects of stadia dynamics that it is hoped to learn about from these data are: -How does the presence of a crowd affect the modal properties of a stadium (i.e.what is the nature of crowdstructure interaction)?-Does a crowd introduce significant non-linearity into the human-structure dynamic system?-What forces do crowds apply to a stadium and how can these be modelled for future designs?
The University of Sheffield is interested in the innovative techniques that researchers around the world can apply to the data in an effort to understand more fully the nature of crowd excitation and interaction on stadia.Hereto, a part of the continuous monitoring data was made available through the university web site [3].Participants had the opportunity to disseminate their results during a special session that is being organised at IMAC 23, the International Modal Analysis Conference, Orlando (FL), 2005 [4].
Originally, this paper was part of that session and focused on the first key question raised above, by introducing the innovative PolyMAX method, which accurately and objectively identifies the modal parameters of a structure excited by ambient excitation [5].Upon request by the editor of this Special Issue of the Journal of Shock and Vibration, the paper has been expanded to include also a discussion on Stochastic Subspace Identification and a brief review of other Operational Modal Analysis methods often encountered in literature.These elements are contained in Section 2. Special emphasis is put on the PolyMAX method, including a discussion on the pre-processing of the raw operational time data and the post-processing of the poles and operational reference factors to obtain the mode shapes.Section 3 discusses the application of Stochastic Subspace Identification and PolyMAX to the stadium vibration data.In particular the efficiency of the methods are highlighted: in a very short time the modes of a structure are identified, even from very large datasets.The main advantage of PolyMAX is that it yields extremely clear stabilisation diagrams, making an automation of the parameter identification process rather straightforward.This enables a continuous monitoring of the dynamic properties of a structure, which is also demonstrated in Section 3.

Operational Modal Analysis
Before the monitoring system was installed, shaker tests were performed to learn more about the road stand [1].As the purpose of this study is to find out how the stadium during a game interacts with the crowd, shaker tests cannot be performed and only ambient excitation can be used.The ambient excitation came from the people occupying the stand, the traffic on the road behind the stand and the wind.It is practically impossible to measure this ambient excitation and, by consequence, the structural response is the only information available to identify the modal parameters of the structure.In this case one speaks of Operational Modal Analysis (OMA).Typical for the output-only case is that the lack of knowledge of the input is justified by the assumption that the input does not contain any information; or in other words, the input is white noise.The theoretical assumption of white noise turns out to be not too strict in practical applications.As long as the (unknown) input spectrum is quite flat, OMA methods will work fine.An overview and comparison of operational modal parameter estimation methods can be found in [6].This section will introduce the new Operational PolyMAX method, and position it with respect to earlier approaches discussed in [6].

Operational PolyMAX method
This section discusses how the PolyMAX method, which is very successful in classical modal analysis [7], can be used for Operational Modal Analysis.

Output-only frequency-domain model
Frequency-domain Operational Modal Analysis methods, such as PolyMAX, require output spectra as primary data.In this subsection, it will be shown that, under the assumption of white noise input, output spectra can be modelled in a very similar way as Frequency Response Functions (FRFs).It is well known that the modal decomposition of an FRF matrix [H(ω)] ∈ C l×m is [8,9]: where l is the number of outputs; m is the number of inputs; n is the number of complex conjugated mode pairs; • * is the complex conjugate of a matrix; • H is the complex conjugate transpose (Hermitian) of a matrix; {v i } ∈ C l are the mode shapes; < l T i >∈ C m are the modal participation factors and λ i are the poles, which are related to the eigenfrequencies ω i and damping ratios ξ i as follows: The input spectra [S uu ] ∈ C m×m and output spectra [S yy (ω)] ∈ C l×l of a system represented by the FRF matrix [H(ω)] are related as: In case of operational data the output spectra are the only available information.The deterministic knowledge of the input is replaced by the assumption that the input is white noise.A property of white noise is that it has a constant power spectrum.Hence [S uu ] in Eq. ( 3) is independent of the frequency ω.The modal decomposition of the output spectrum matrix is obtained by inserting Eq. ( 1) into Eq.( 3) and converting to the partial fraction form [10][11][12]: where < g i >∈ C l are the so-called operational reference factors, which replace the modal participation factors in case of output-only data.Their physical interpretation is less obvious as they are a function of all modal parameters of the system and the constant input spectrum matrix.Note that the order of the power spectrum model is twice the order of the FRF model.The goal of OMA is to identify the right hand side terms of Eq. ( 4) based on measured output data pre-processed into output spectra (Section 2.1.2).

Pre-processing operational data
Power spectra are defined as the Fourier transform of the correlation sequences.The most popular non-parametric spectrum estimate is the so-called weighted averaged periodogram (also known as modified Welch's periodogram).Weighting means that the signal is weighted by one of the classical windows (Hanning, Hamming, . ..) to reduce leakage.Welch's method starts with computing the discrete Fourier transform (DFT) Y (ω) ∈ C l of the weighted outputs: where y k ∈ R l is the sampled output vector (the measurements); w k denotes the time window; k is the time instant; N the number of time samples and ∆t the sampling time.An unbiased estimate of the spectrum is the weighted periodogram: The variance of the estimate is reduced by splitting the signal in possibly overlapping blocks, computing the weighted periodogram of all blocks and taking the average: where P is the number of blocks and superindex b denotes the block index.Another non-parametric spectrum estimate is the so-called weighted correlogram.It will be shown that this estimate has some specific advantages in a modal analysis context.First the correlation sequence has to be estimated (i being the time lag): High-speed (FFT-based) implementations exist to compute the correlations as in Eq. ( 8); see for instance [13].The weighted correlogram is the DFT of the weighted estimated correlation sequence: where L is the maximum number of time lags at which the correlations are estimated.This number is typically much smaller than the number of data samples to avoid the greater statistical variance associated with the higher lags of the correlation estimates.In a modal analysis context, the weighted correlogram has the following advantages: -It is sufficient to compute the so-called half spectra which are obtained by using only the correlations having a positive time lag in Eq. ( 9): The relation between the half spectra Eq. ( 10) and the full spectra Eq. ( 9) is the following: It can be shown (see for instance [10,14]) that the modal decomposition of these half spectra only consists of the first two terms in Eq. ( 4): The advantage in modal analysis is that models of lower order can be fitted without affecting the quality.
-Under the white noise input assumption, the output correlations are equivalent to impulse response.So, just like in impact testing, it seems logical to apply an exponential window w k to the correlations before computing the DFT Eq. (10).The exponential window reduces the effect of leakage and the influence of the higher time lags, which have a larger variance.Moreover, the application of an exponential window to impulse responses or correlations is compatible with the modal model and the pole estimates can be corrected.This is not the case when a Hanning window is used: such a window always leads to biased damping estimates.
Interesting discussions on the estimation of spectra from the measured time histories in a modal analysis context can be found in [10,14,15].
The weighted correlogram approach is illustrated using the Bradford stadium vibration data (More details on the data can be found in Section 3). Figure 2 shows the time history of a typical acceleration channel.Figure 3 (Left) shows how thousands time data samples are reduced to 256 correlation samples using Eq.(8).The effect of the exponential window is visible.Figure 3 (Right) shows the half spectrum estimate according to Eq. (10).

PolyMAX modal parameter estimation
After having pre-processed output data into output spectra (Section 2.1.2),it is now the task to identify a modal model (Section 2.1.1).By comparing Eqs ( 12) with (1), it is clear that FRFs and half spectra can be parameterised in exactly the same way.By consequence, the same modal parameter estimation methods can be used in both cases.PolyMAX is such a method.It will be reviewed in the following.
PolyMAX originates from the so-called least-squares complex frequency-domain (LSCF) estimation method which was introduced to find initial values for the iterative maximum likelihood method [16].The LSCF method estimates a so-called common-denominator transfer function model [17].Quickly it was found that these "initial values" yielded already very accurate modal parameters with a very small computational effort [16,18,19].The most important advantage of the LSCF estimator over the available and widely applied parameter estimation techniques [9] is the fact that very clear stabilisation diagrams are obtained.
It was found that the identified common-denominator model closely fitted the measured frequency response function (FRF) data.However, when converting this model to a modal model Eq. ( 1) by reducing the residues to a rank-one matrix using the singular value decomposition (SVD), the quality of the fit decreased [18].Another feature of the common-denominator implementation is that the stabilisation diagram can only be constructed using pole information (eigenfrequencies and damping ratios).Neither participation factors nor mode shapes are available at first instance.The theoretically associated drawback is that closely spaced poles will erroneously show up as a single pole [20].These two reasons provided the motivation for developing "PolyMAX", which is a polyreference version of the LSCF method, using a so-called right matrix-fraction model.In this approach, also the participation factors are available when constructing the stabilisation diagram.The main benefits of the polyreference variant are the facts that the SVD step to decompose the residues can be avoided and that closely spaced poles can be separated.PolyMAX was introduced and validated using industrial data sets in [7,20].The present paper will indicate how PolyMAX can also be applied to operational data when appropriate pre-processing (Section 2.1.2) and post-processing (Section 2.1.4)is applied.
In the PolyMAX method, following so-called right matrix-fraction model is assumed to represent the measured half spectrum matrix: where [β r ] ∈ C l×m are the numerator matrix polynomial coefficients; [α r ] ∈ C m×m are the denominator matrix polynomial coefficients; p is the model order and m is in this case not the number of inputs as in Eq. ( 1), but the number of outputs selected as references.For estimating the modal parameters, it is indeed sufficient to compute only the cross spectra between all outputs l and a limited set of references m.From this point forward, S + yy (ω) is an l × m matrix instead of an l × l as above.Please note that a so-called z-domain model (i.e. a frequency-domain model that is derived from a discrete-time model) is used in Eq. ( 13), with z = exp(jω∆t).Equation ( 13) can be written down for all values ω of the frequency axis of the half spectra data.Basically, the unknown model coefficients [α r ] , [β r ] are then found as the least-squares solution of these equations (after linearisation).More details about this procedure can be found in [7,20].
Once the denominator coefficients [α r ] are determined, the poles λ i and operational reference factors < g i > are retrieved as the eigenvalues and eigenvectors of their companion matrix, which is rather classical [9,21].A p th order right matrix-fraction model yields pm poles.This procedure allows constructing a stabilisation diagram [22] for increasing model orders p and using stability criteria for eigenfrequencies, damping ratios and operational reference factors.An efficient construction of the PolyMAX stabilisation diagram is possible by formulating the least squares problem for the maximum model order p max .Considering submatrices of appropriate dimensions can then solve lower-order problems.Examples of PolyMAX stabilisation diagrams will be shown in Section 3.

Operational residuals
The interpretation of the stabilisation diagram yields a set of poles and corresponding operational reference factors.The mode shapes can then be found from Eq. ( 12): in which LR, U R ∈ R l×m , respectively the lower and upper residuals, have been introduced to model the influence of the out-of-band modes in the considered frequency band.The only unknowns in Eq. ( 14) are the mode shapes {v i } and the lower and upper residuals.They are readily obtained by solving Eq. ( 14) in a linear least-squares sense.This second step is commonly called least-squares frequency-domain (LSFD) method [9,21].The operational residuals LR, U R in Eq. ( 14) are different from the usual residuals in case of FRFs.They were determined by verifying the asymptotic behaviour of the output spectra of a single-degree-of-freedom (SDOF) mechanical system excited by white noise (Fig. 4).The formats of the residuals that have to be applied in the LSFD method are listed in Table 1 for FRFs Eq. ( 1), classical full output spectra Eq. ( 4) and half spectra Eq. (12).In general, they depend on the measured output quantity (displacement -velocity -acceleration).Deriving displacements to velocities (and similarly velocities to accelerations) implies multiplication by jω for FRFs and (jω) 2 for full spectra.Rather surprising is that this is not the case for half spectra: the residuals are independent of the output quantity.

Stochastic subspace identification 2.2.1. Stochastic state space model
A rather advanced yet popular OMA method is Stochastic Subspace Identification (SSI).In this method a socalled stochastic state space model Eq. ( 15) is identified from output correlations or directly from measured output data [23][24][25][26].It can be shown that a stochastic state space model is a good representation of a vibration structure excited by some unknown forces, which are assumed to be white noise signals (see for instance [12,26] for a detailed derivation): where y k ∈ R l is the sampled output vector (the measurements); x k ∈ R 2n is the discrete state vector; w k ∈ R 2n is the process noise, typically due to disturbances and modelling inaccuracies, but here also due to the unknown excitation of the structure; v k ∈ R l is the measurement noise, typically due to sensor inaccuracy, but here also due to the unknown excitation of the structure; k is the time instant; l is the number of outputs; 2n is the system order.The matrix A ∈ R 2n×2n is the state transition matrix that completely characterises the dynamics of the system by its eigenvalues; and C ∈ R l×2n is the output matrix that specifies how the internal states are transformed to the outside world.Note that the system order 2n in Eq. ( 15) is twice the number of complex conjugated pole pairs n in Eq. ( 1).
The only knowns in Eq. ( 15) are the output measurements y k .The challenge is now to determine the other parameters and in particular the system matrices A,C form which the modal parameters can be derived (Section 2.2.3).

Pre-processing operational data
For reasons of simplicity, only the correlation-driven variant of SSI will be discussed.Hence, the raw time measurements need to be pre-processed into correlations.In fact, this was already explained in Section 2.1.2as computing correlations was an intermediate step in obtaining the half spectra required in the PolyMAX approach.As stated there, correlations can be computed in a highly efficient way using the FFT implementation of Eq. ( 8).Note that there is no need to apply an exponential window in case of the SSI method.

SSI modal parameter estimation
Considering the definitions of the i th lag correlation matrix R i ∈ R l×l and the next state -output correlation matrix G ∈ R 2n×l : where E is the expected value operator, it is rather straightforward to prove that the output correlations decompose as (see for instance [25,27]): The correlation sequence R i can be estimated from measurement data as described in Section 2.2.2; so it remains to be solved how to decompose the correlations as to recover A, C, G. Hereto, a block Toeplitz matrix T 1|i ∈ R li×li is formed that consists of the correlation estimates (Definition: a Toeplitz matrix is constant across the diagonals): The second equality follows from applying property Eq.(17).The definitions of the extended observability matrix O i ∈ R li×2n and the reversed extended stochastic controllability matrix Γ i ∈ R 2n×li are obvious from Eq. (18).For li > 2n, and in case of an observable and controllable system, the rank of the li × li Toeplitz matrix equals 2n, since it is the product of a matrix with 2n columns and a matrix with 2n rows.The actual implementation of SSI consists of computing the Singular Value Decomposition (SVD) of T 1|i , truncating the SVD to the model order 2n, estimating O i and Γ i by splitting the SVD in two parts and finally estimating A, C, G from O i and Γ i in a least-squares sense.Implementation details can be found in the above-cited references.It is straightforward to slightly reformulate the SSI method, so that it only needs the correlations between all outputs and a sub-set of references [24].
Finally, the modal parameters can be computed from the system matrices A, C. The derivation starts with the eigenvalue decomposition of A: where Ψ ∈ C 2n×2n is the eigenvector matrix and Λ d ∈ C 2n×2n is a diagonal matrix containing the discrete-time eigenvalues µ i , which are related to the continuous-time eigenvalues λ i as: The eigenfrequencies and damping ratios are related to λ i as indicated in Eq. ( 2).Finally, the mode shapes V ∈ C 2n×2n are found as: In theory the system order 2n can be determined by inspecting the number of non-zero singular values of T 1|i Eq. ( 18).In practice however, the estimated Toeplitz matrix is affected by "noise" (modelling inaccuracies, measurement noise, computational noise, finite data length, . ..) leading to singular values that are all different from zero.Sometimes it is suggested to look at the "gap" between two successive singular values.The singular value where the maximal gap occurs yields the model order.This criterion should however not be applied too dogmatically, since in most  practical cases there is no gap.However, to obtain a good model for modal analysis applications, it is probably a better idea to construct a stabilisation diagram, by identifying a whole set of models with different order.The stabilisation diagram was already introduced in Section 2.1.3.In case of the SSI method, an efficient construction of the stabilisation diagram is achieved by computing the SVD of the Toeplitz matrix only once.Models of different order are then obtained by including a different number of singular values and vectors in the computation of O i and Γ i .An example of an SSI stabilisation diagram will be shown in Section 3.

Other methods
In this section, a brief overview of other OMA methods often encountered in literature will be given.An overview of names and acronyms can be found in Table 2.

Peak-picking
The simplest method to obtain rough estimates of the modal parameters of a structure subjected to ambient loading is the so-called peak-picking method.The method is named after the key step of the method: the identification of the eigenfrequencies as the peaks of power spectrum plots.A violation of the basic assumptions (low damping and well-separated frequencies) leads to erroneous results.In fact the method identifies operational deflection shapes (ODS) instead of mode shapes and for closely spaced modes such an ODS will be the superposition of multiple modes.Other disadvantages are that the selection of the eigenfrequencies can become a subjective task if the spectrum peaks are not very clear and that the method does not yield any damping estimates.
Evidently, inspecting the peaks in spectrum plots remains a very useful concept, for instance to confront the eigenfrequency estimates of the more advanced curve-fitters with.

The Complex Mode Indication Function (CMIF)
A more advanced frequency-domain OMA method consists of computing the SVD of the output spectrum matrix.Various literature references exist on the method, including the application to FRF data in classical modal analysis [28][29][30][31].Also various names are circulating: "method based upon the diagonalisation of the spectral density matrix," the Complex Mode Indication Function (CMIF), and the Frequency Domain Decomposition (FDD) method.
The method is based on the fact that the transfer function or spectrum matrix evaluated at a certain frequency is only determined by a few modes.The number of significantly contributing modes determines the rank of the spectrum matrix.The SVD is typically used for estimating the rank of a matrix: the number of non-zero singular values equals Table 3 Comparison of modal parameters from PolyMAX and Stochastic Subspace Identification (SSI).Eigenfrequencies (f ) and damping ratios (ξ) for both sets are given as well as their differences (∆).MP C stands for "mode phase colinearity" and expresses the complexity of a mode shape.MAC stands for "modal assurance criterion" and equals the squared correlation between two mode shapes.the rank [32].Therefore, plotting the singular values as a function of frequency yields the eigenfrequencies as local maxima.The method is also able to detect closely spaced modes: more than one singular value will reach a local maximum around the close eigenfrequencies.The first singular vector at resonance is an estimate of the mode shape at that frequency.In case of mode multiplicity at a resonance frequency, every singular vector corresponding to a non-zero singular value is considered as a mode shape estimate.
Extensions of the CMIF method are possible that do estimate eigenfrequencies and damping ratios differently as in the PP method.After applying the SVD to the spectrum matrix, this matrix is in fact decomposed in SDOF systems.To such a system, SDOF modal parameter estimation methods could be applied, extensively documented in the modal analysis literature [8,9].
The main disadvantage of the peak-picking method remains: also in the CMIF method the initial selection of eigenfrequencies (as peaks in singular value plots) is a rather subjective task.This is one of the drivers behind the development of more advanced parametric methods (such as PolyMAX, SSI, MLE), that fit models to the measured data and allow a more objective selection of model order and physical poles based on, for instance, stabilisation diagrams.

Maximum likelihood estimator (MLE)
Contrary to the PP method or the CMIF, which consider only one mode at a time, this method estimates the parameterised spectrum matrix as a whole.Maximum Likelihood (ML) identification is an optimisation-based method that estimates the parameters of a model by minimising an error norm.A discussion on the use of the ML estimator to identify parametric frequency-domain models can be found in [33,34].The ML method results in equations that are non-linear in the unknown parameters.This requires an iterative procedure.Therefore it is no surprise that often mentioned drawbacks of ML estimators are the high computational load and the fact that they are not suited to handle large amounts of data.During the last years attention has been paid to the optimisation of the ML method: the algorithm has been modified to keep the memory requirements as low as possible; and using an adapted parameterisation and fast signal processing techniques, an important reduction of the computation time was possible.It has been shown that ML identification is a robust method to find the modal parameters of a structure from a large and noisy data set [35,36].Originally intended for application to FRFs, the method was extended to use spectra as primary data, so that it also could be used in output-only cases [37].Very promising are also the "MLE-alike" methods that complement linear least-squares estimators with one additional iteration: a single iteration keeps the computation time low, but still improves greatly the accuracy of the modal parameter estimates as demonstrated in [38].

Stadium data analysis
The stadium structure was introduced in Section 1.The University of Sheffield made a 4-hours recording of stadium accelerations available [3].The data acquisition started about 1 hour before the football game between Bradford City and Sheffield United on 23 November 2002.The data include following events (Fig. 5): -Filling up of the stand These events represent different crowd-structure interaction scenarios and allow investigating how the presence of a crowd affects the modal properties of the stadium.The original data were sampled at 80 Hz.As we are mainly interested in the modal parameters below 10 Hz, the data were resampled to 20 Hz.By consequence, 288000 samples (20 Hz × 14400 s) are available for each of the 12 accelerometers.The "sliding" RMS values represented in Fig. 5 are based on a frame of 1000 samples and an overlap between the frames of 50%.

Analysis of the complete data set
Although it is expected that the different crowd configurations identified in Fig. 5 will influence the modal parameters of the stadium [39], a first analysis will be performed on the complete data set.This allows illustrating the efficiency of the PolyMAX operational modal parameter estimation process and getting an idea of the modes that can be expected from a more detailed analysis.Evidently, the risk of inconsistencies exits due to changing modal parameters during the data acquisition.
Two sensors are selected as reference sensors and, as explained in Section 2.1.2and illustrated in Figs 2 and 3, the 12 raw time histories are reduced to a 12 × 2 cross correlation Eq. ( 8) and cross spectrum matrix Eq. ( 10).An exponential window of 1% has been applied to the cross-correlations before computing the half spectra by a DFT.The number of time lags at which the correlations have been computed was 256.Afterwards, the PolyMAX method was applied to estimate the modal parameters.The whole process required only 4 s of computation time on a year-2003 PC (Fig. 6). Figure 7 compares the PolyMAX (Section 2.1) stabilisation diagram with the one obtained using SSI (Section 2.2).Although both methods are known to yield accurate modal parameter estimates, it is obvious that PolyMAX considerably facilitates the identification process; in the sense that it is much easier to pick the stable poles from the diagram.The symbols indicating the position (i.e.frequency versus model order) of the poles in Fig. 13.PolyMAX stabilisation diagrams for different events: Filling -1 st goal -sitting -leaving.a stabilisation diagram represent the "stability" of a pole with respect to the closest pole one model order below.Stability bounds are defined for frequency, damping and eigenvector values and the following symbols are used: "f" for frequency stabilisation, "d" for damping and frequency or "v" for eigenvector and frequency and "s" for full stabilisation.As typical stability criteria, following values are used: -1% for frequency stability, -5% for damping stability, -2% for eigenvector stability.
These "defaults" reflect the accuracy of the estimates that can be expected in a wide range of industrial modal analysis applications.
Section 3.2 will compare PolyMAX and SSI results in more detail.As an overall quality indicator of the parameter estimation process, the sum 1 of measured cross spectra is compared with the sum of spectra that are synthesised from the identified modal parameters using Eq. ( 14); see Fig. 8.The good correspondence indicates that all major dynamic properties have been extracted from the data.

Empty stadium
In this section, we will discuss in detail the OMA of the empty stadium.The data used here involve the last 2400 s of the records, as indicated in Fig. 2 and pre-processed in Fig. 3. Figure 9 shows the PolyMAX stabilisation diagram and the quality of the parameter estimation process by comparing the sum of measured and synthesised cross spectra.Again, the correspondence is excellent.It is noteworthy that not only the absolute levels of vibration are much smaller when comparing the spectra of the empty stadium (Fig. 9, Right) with the spectra computed from the complete data set (Fig. 8), but that also the relative contribution of the modes is changing.So, it can be expected that, depending on the selected event (Fig. 5), the modes will be better or less well excited and the mode extraction process will be easier or, on the contrary, more difficult.Some of the identified mode shapes are represented in Fig. 10.The correspondence between measured and synthesised half spectra (Fig. 9, right and Fig. 11) already confirms the accuracy of the PolyMAX estimates, but it still interesting to compare them with estimates from the SSI method.Table 3 contains this comparison.The eigenfrequency agreement is excellent, whereas the damping ratio differences fall within the typical uncertainty on these estimates.Both PolyMAX and SSI estimate modes with low complexity as indicated by the M P C values close to 100%.The correspondence between the mode shapes is very good (very high M AC values in last column).The complete matrix of M AC values is represented in Fig. 12.

Other events
Based on Fig. 5, certain events are isolated and an Operational Modal Analysis using the PolyMAX method was performed.The pre-processing and processing parameters were the same in all cases: correlations up to lag 256 were computed, an exponential window of 1% was applied, and the PolyMAX method was running up to model order 60 between 2 and 8 Hz.The number of correlation time lags determines the frequency resolution of the spectra and the exponential window strength their "quality" (see discussion in Section 2.1.2). Figure 13 shows the PolyMAX stabilisation diagrams for the events: filling up -1 st goal -sitting -leaving.Figure 14 compares the sum of measured and synthesised cross spectra for all these events.Again it can be concluded that PolyMAX is able to extract the main dynamic features from the data of all events.

Automatic modal analysis
As discussed earlier in this paper, the main advantage of PolyMAX is that it yields extremely clear stabilisation diagrams.This makes an automation of the parameter identification process rather straightforward and enables a continuous monitoring of the dynamic properties of a structure.An overview and discussion on automatic modal analysis techniques can be found in [22].
In this section, a method has been used which is able to autonomously interpret stabilisation diagrams based on a set of heuristic rules [12,40].Unlike in previous sections, no attempts were made to isolate certain events, but the data were simply cut in segments of equal length.Two different segmentations were applied: -Segments of 4096 samples (204.8 s @ 20 Hz) without overlap -Segments of 16384 samples (819.2 s @ 20 Hz) with 75% overlap, to obtain the same time resolution as in previous case.
Figure 15 shows the automatic PolyMAX results.From these diagrams following observations can be made: -The 3 lowest modes are easiest to track.
-The higher modes can be tracked when the stadium is filling up or when it is emptying or empty.
-During the football game, the higher modes are not well excited and/or not consistently estimated and/or not well tracked.-Enlarging the segment length yields more stable results (Bottom diagram of Fig. 15 vs. Top diagram) and an improved tracking.Of course, larger segments will make it more difficult to identify changes in modal parameters that only occur during a short time.
Figure 16 shows the tracking of mode 3 based on Fig. 15.During the game, when people are sometimes heavily interacting with the structure, the frequencies tend to decrease and the damping is increasing.The estimated mode shapes may also be change during the game, as evidenced by the drops in MAC values.The most stable conditions occur when the stadium is empty, filling up with people before the game or when it is emptying after the game.
Despite the fact that the analysis of individual events is quite successful (Section 3.3), tracking is difficult for some of the modes.Possible causes are that modes change too much at certain crowd-structure interaction events or that they may not be well excited in certain data segments.For instance, the crowd is mainly exciting the lower modes.

Conclusions
The paper discussed the use of the Operational Modal Analysis for monitoring a stadium structure during a football game.The Operational PolyMAX and Stochastic Subspace Identification method were discussed in some detail and an overview of other commonly used methods was given.The main advantage of PolyMAX is that it yields very clear stabilisation diagrams, easing dramatically the problem of selecting the model order and the structural system poles.The analysis results compare favourably with current best-of-class commercially available methods, without increasing the computational effort.
Although having a completely different mathematical background, both PolyMAX and Stochastic Subspace Identification approximately yield the same modal parameter estimates.The main difference is that PolyMAX is easier to automate.Such an automated processing was successfully applied to the whole data set split in different segments.Although the theoretical white noise excitation assumption is very questionable when people are interacting with the structure, the application of Operational Modal Analysis confirmed that the frequencies of the structure are decreasing when the number of people is increasing (added mass) and that the damping ratios are generally increasing.Both PolyMAX and the Stochastic Subspace Identification ("Op.Time MDOF") are implemented in the LMS Test.Lab Structures software [41].

Fig. 2 .
Fig. 2. Stadium ambient acceleration time data.(Top) 4 hours recording; (Bottom) selection of the last 2400 s for further processing.

Fig. 4 .
Fig. 4. SDOF response to white noise input.(Top) Full spectrum; (Bottom) half spectrum.(Left) displacement; (Middle) velocity; (Right) acceleration response.The upper and lower residuals can be represented by straight lines on a log-log scale.

Fig. 5 .Fig. 6 .
Fig. 5. Typical acceleration time history and "sliding" RMS value with an indication of the events during the game.

Fig. 12 .
Fig. 12. MAC values assessing the mode shape correlation between SSI and PolyMAX mode shapes.-People seated during the game -Half time -Emptying of the stand -Empty stand -5 goals (by the way, all scored by Sheffield United)

Fig. 15 .
Fig.15.Unusual stabilisation diagrams representing automatic modal analysis results: the y-axis is not the model order, but represents the time.Automatic PolyMAX analyses were performed on data segments of 204.8 s and no overlap (Top) and 819.2 s with 75% overlap (Bottom).The background image is the normalised 1 st singular value (CMIF) of the spectrum matrix computed for each segment.The curves on the Left represent the sliding RMS values (Fig.5).

Table 2
Overview of parameter estimation methods that are found in OMA literature.An acronym and a name designate the methods.Alternative names, which can be found in literature for basically the same methods, are also represented in the table