# Cross-Validation in Stochastic Subspace Identification

Simon Marwitz, Tom Lahmer, Volkmar Zabel

► To cite this version:

Simon Marwitz, Tom Lahmer, Volkmar Zabel. Cross-Validation in Stochastic Subspace Identification. 11th International Operational Modal Analysis Conference, Institut National de Recherche en Informatique et en Automatique [Inria], May 2025, Rennes, France. hal-05398773

HAL Id: hal-05398773

https://hal.science/hal-05398773v1

Submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cross-Validation in Stochastic Subspace Identification

Simon Marwitz<fnref n="1" />, Tom Lahmer<fnref n="2" />, Volkmar Zabel<fnref n="3" />

## ABSTRACT

Operational Modal Analysis (OMA) often relies on user-specified parameters, such as model order and input sequence lengths. These are typically selected based on expert knowledge. While stabilization diagrams effectively determine the model order, they are unsuitable for selecting other sensitive parameters. The existing criteria for the assessment of parameter sets in OMA fall short of meeting key requirements, including the independence from system dynamics, and applicability across different OMA methods. This study proposes a common validation criterion based on input reconstruction and cross-validation techniques. The covariance- and data-driven stochastic subspace identification methods, are extended for the purposes of multi-block estimation and synthesis. Moreover, a modally decomposed formulation is selected to facilitate an additional assessment of the identified modal parameters.

The proposed method is verified through the use of synthetic and experimental vibration signals, as well as a full-factorial parameter study, which collectively demonstrate the robustness and consistency of the proposed criterion.

The key objective of the developed methodology is to enhance confidence in identified modal parameters, thereby facilitating more reliable assessments in e. g. structural health monitoring or life-cycle analyses.

Keywords: Cross Validation, Synthetization, Reconstruction, Uncertainty, Stochastic Subspace Identification, Quality

## 1. INTRODUCTION

Operational Modal Analysis (OMA) applications frequently depend on user-specified parameter sets, such as the order of identified models and the length of input sequences, among other factors. These parameters are often selected based on expert knowledge and experience. Although the model order can be determined using a stabilization diagram, this method is not sufficient for selecting other sensitive parameters.

<footnote n="1">¹ Bauhaus-Universität Weimar, simon.jakob.marwitz@uni-weimar.de</footnote>

<footnote n="2">² Bauhaus-Universität Weimar, tom.lahmer@uni-weimar.de</footnote>

<footnote n="3">³ Universität Rostock, volkmar.zabel@uni-rostock.de</footnote>
Criteria for the assessment of a chosen set of parameters for OMA exist under the name of modal indicators, mode validation criteria, etc. However, they are typically

* specific to a single OMA method, e. g. the Consistent Mode Indicator [1] for the ERA, or

* specific to particular structures, e. g. the Modal Phase Colinearity [1] is only applicable to proportionally damped structures with properly spaced modes,

* not normalized, e. g. the Mean Phase Deviation [2] may take any value in [0, 180]°, or

* suffer from overfitting and / or sample bias, e. g. the Modal Contributions [3] use all data for training and validation.

The proposed criterion is derived from the identified systems in various OMA methods, assuming a direct correspondence to the modal parameters. The reconstruction error [4, 5] is a common quality measure based on reconstructed input signals, synthesized by the identified system. To mitigate overfitting and/or sample bias if the same input signals are used for identification and validation, cross-validation methods [6], a class of well-known techniques in statistics, are employed.

In this study, synthetization methods for two stochastic state-space models are extended to accommodate multi-block estimation and cross-validation. The method is verified on synthetic signals and demonstrated with experimental vibration measurements from a guyed tower. A full-factorial parameter study demonstrates the unbiased behavior and numerical robustness of the developed criterion.

## 2. CROSS-VALIDATION (CV) TECHNIQUES

The reconstruction error of any system identification method can be estimated by comparing the identification input sequences with those synthesized by the identified system. However, it is recognized that training and evaluating a model on the same data can lead to over-optimistic validation estimates. CV methods [6] try to overcome these limitations. These methods divide the available data into a training set, which is used for system identification, and a validation set, against which the reconstructed system output is validated. The input and output of the system are not necessarily required to be timestep signals; they may also be correlation functions, spectral densities, or other derived quantities, depending on the identification method used.

The underlying idea of CV is to divide the available input data into $n_{bl}$ blocks and subsequently assemble them in a set $D = \{d_1, ..., d_{n_{bl}}\}$. A portion of these blocks is allocated for training a model (the training sample $D_t \subset D$), while the remainder is used for evaluating its performance (the validation sample $D_v = D_t^c$). A measure of fit $\gamma(M_{D_t}; d_i)$ quantifies the discrepancy between the output of an identified model $M_{D_t}$ and a validation sample $d_i \in D_v$. One such measure is developed in Sect. 4. The loss function is defined as the expected value of the validation measure

$$ \mathcal{L}(M_{D_t}, D_v) = \mathbb{E}[\gamma(M_{D_t}, D_v)]. \quad (1) $$

Different estimators exist depending on the employed method for assembling the data blocks into training and validation sets.

The simple Hold-Out (HO) estimator consists of a single ($b = 1$) training and validation set. The validation results are then averaged over the validation samples

$$ \hat{\mathcal{L}}_{\text{HO}}(M_{D_t}, D_v) = \frac{1}{n_v} \sum_{d_i \in D_v} \gamma(M_{D_t}, d_i). \quad (2) $$

However, the validation results are dependent on the selection of the training and validation sets, which introduces sample bias. To circumvent this, CV methods assemble the training and validation sets in $b$ distinct ways and average the validation measures. The Leave-One-Out and Leave-$p$-Out validators are examples of this category. The former employs all but one set for training and a single validation set, repeating this process across each set as the validation set, i. e. $b = n_{bl}$. It should be noted that both of
these methods are exhaustive, and as such, they may become computationally expensive if the training and/or validation procedure is so. A balanced choice that is not exhaustive but still reduces sample bias is the k-fold CV method. It involves the grouping of the data sets into k approximately equally sized validation subsets which are then successively used for training and validation. The CV estimator is subsequently averaged over each of these k groups.

$$ \hat{\mathcal{L}}_k(\mathcal{M}; D) = \frac{1}{k} \sum_{j=1}^{k} \hat{\mathcal{L}}_{\text{HO}}(\mathcal{M}_{D_{t,j}}; D_{v,j}). \quad (3) $$

### 3. MULTI-BLOCK STOCHASTIC SUBSPACE IDENTIFICATION

The Stochastic Subspace Identification (SSI) technique is a well-known method in operational modal analysis for identification of discrete, stochastic state-space models from measured signals. For a detailed description of the algorithms used in this study the reader is referred to [7–9]. The discrete stochastic state space system is defined by

$$ \begin{aligned} \mathbf{x}[n+1] &= \mathbf{A}\mathbf{x}[n] + \mathbf{w}[n], & \text{the state equation, and} \\ \mathbf{y}[n] &= \mathbf{C}\mathbf{x}[n] + \mathbf{v}[n], & \text{the output equation,} \end{aligned} \quad (4) $$

where $\mathbf{x}(t)$ is the state vector, $\mathbf{A}$ is the state matrix, $\mathbf{C}$ is the output matrix, and the stochastic white noise vectors $\mathbf{w}[n]$ and $\mathbf{v}[n]$ account for the unknown input terms. The eigendecomposition $\mathbf{A}\mathbf{\Psi} = \mathbf{\Psi}\mathbf{\Lambda}$, is the basis for the estimation of modal frequencies $f_j$, modal damping ratios $\zeta_j$, and mode shapes $\phi_j$, as well as the modal decoupling of the state space system.

The required extensions to apply multi-block cross-validation techniques with the SSI are given in the following. In the covariance-driven case, a block-modified form of the unbiased correlation function estimator between the channels $i_l = 1 \dots n_l$ and reference channels $i_r = 1 \dots n_r$ of the input signal $\mathbf{y}[n]$ is

$$ \hat{R}_{i_l i_r, i_{bl}}[m] = \frac{1}{N_{bl} - m} \sum_{n=0}^{N_{bl} - m - 1} y_{i_l}[i_{bl}N_{bl} + n] y_{i_r}[i_{bl}N_{bl} + m], \quad (5) $$

where $i_{bl} = 1 \dots n_{bl}$ is the block number, $N_{bl}$ is the block length in timesteps, and $m = 1 \dots M$ is the lag step of the correlation function. For short block lengths $N_{bl}$ and high lags $m$, large confidence intervals may occur [10]. These can be circumvented by allowing the second input signal $y_{i_r}[n+m]$ to overlap into the consecutive block by $m$ timesteps, at the cost of reduced block independence. The identification of $\mathbf{A}$, $\mathbf{C}$ and $\mathbf{G}_{\text{ref}} = \mathbb{E}[\mathbf{x}[n+1] \mathbf{y}_{\text{ref}}[n]^\top]$ then proceeds as shown in [7] using an averaged correlation function estimate over a (training) subset of these $n_{bl}$ correlation function estimates.

In the data-driven case, the projection matrix [11] must be estimated from multiple blocks of signals. Döhler [9] has developed a method to compute independent and identically distributed estimates of reduced projection matrices. This method is extended here for full projection matrices, which are required for input reconstruction.

A block Hankel matrix of the so-called past and future signals is split into $n_{bl}$ blocks along its block-columns as shown in [9] such that

$$ \begin{bmatrix} \mathbf{Y}_{-,1} & \mathbf{Y}_{-,2} & \cdots & \mathbf{Y}_{-,n_{bl}} \\ \mathbf{Y}_{+,1} & \mathbf{Y}_{+,2} & \cdots & \mathbf{Y}_{+,n_{bl}} \end{bmatrix} = \begin{bmatrix} \mathbf{L}_1 & \mathbf{L}_2 & \cdots & \mathbf{L}_{n_{bl}} \end{bmatrix} \begin{bmatrix} \mathbf{Q}_1 & 0 & \cdots & 0 \\ 0 & \mathbf{Q}_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathbf{Q}_{n_{bl}} \end{bmatrix}, \quad (6) $$

where $\mathbf{L}_j$, $\mathbf{Q}_j$ are the individual LQ decompositions of each block.

Another thin LQ decomposition of the column-stacked triangular matrices is defined

$$ [\mathbf{L}_1 \ \mathbf{L}_2 \ \cdots \ \mathbf{L}_{n_{\text{bl}}}] = \tilde{\mathbf{L}}\tilde{\mathbf{Q}} \quad (7) $$
where uniqueness is achieved by restraining the diagonal of $\tilde{L}$ to positive values, i. e. by multiplication of the rows / columns of $\tilde{L}, \tilde{Q}$ with sgn(diag($\tilde{L}$)). Döhler [9] has shown, that

$$L = \tilde{L}, \quad Q = \tilde{Q} \begin{bmatrix} Q_1 & 0 & \cdots & 0 \\ 0 & Q_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & Q_{n_{\text{bl}}} \end{bmatrix} \qquad (8)$$

can be assumed. The identification then proceeds as detailed in [7] using only a (training) subset of these $n_{\text{bl}}$ blocks in Eq (6). The system $A, C$ is finally identified by solving the following relation [7]

$$\begin{bmatrix} N^{-1/2} & \hat{X}_{q+1} & \cdots & y[q+N] \end{bmatrix} = \begin{bmatrix} A \\ C \end{bmatrix} \hat{X}_q + \begin{bmatrix} \rho_w \\ \rho_v \end{bmatrix}, \qquad (9)$$

in a least-squares sense. The residuals yield the noise correlations $\Sigma_{ww}, \Sigma_{wv}, \Sigma_{vv}$, which are required for the construction of the Kálmán gain in the next section.

## 4. (MODAL) RECONSTRUCTION OF INPUT SEQUENCES

The states and corresponding system matrices are unique up to within a similarity transform [11]. Thus, a modally decoupled identified stochastic forward state-space model can be constructed as follows

$$\begin{align} x_m[n+1] &= \Lambda x_m[k] + w_m[n] \notag \\ y[n] &= \Phi x_m[n] + v_m[n]. \tag{10} \end{align}$$

Similarly, a forward innovation model, which is obtained by applying a Kálmán filter to the forward model, can be transformed. By eliminating the innovation sequences, the modally decoupled forward innovation model is expressed as

$$\begin{align} \hat{x}_m[n+1] &= (\Lambda - K_m C_m) \hat{x}_m[n] + K_m y[n] \tag{11} \\ y[n] &= \tilde{y}_m[n] + \tilde{y}_e[n]. \notag \end{align}$$

In both cases the modal states and system matrices are

$$x_m[n] = \Psi^{-1}x[n], \quad A = \Psi\Lambda\Psi^{-1}, \quad C_m = C\Psi, \text{ and } K_m = \Psi^{-1}K. \qquad (12)$$

**Modal Decomposition of Correlation Functions** Given a modally decoupled identified stochastic forward state-space model, Reynders [4] has demonstrated the synthetization of correlation functions

$$\tilde{R}_m[m] = \sum_{j=1}^{n_m} \left( \lambda_j^{m-1} \phi_j g_{m,j}^{\top} + \lambda_{j^*}^{m-1} \phi_{j^*} g_{m,j^*}^{\top} \right), \qquad (13)$$

where $n_m$ the total number of modes, $j = 1... n_m$ is the mode number and, by slight abuse of notation, $j^*$ the index of its complex conjugate counterpart. The modal participation vector $g_{m,j}$ is the $j^{\text{th}}$ row of the modally decomposed next-state output covariance matrix, that is obtained by solving $G_{\text{ref}} = \Psi G_m$ [12]. An example will be introduced in Sec. 5. and is shown in Fig. 1.

**Modal Decomposition of Kálmán-Filter State Sequences** The modal decoupling of the Kálmán filter states of the input signal allows the estimation of each mode's contribution to the signal as shown by [3]. Under stationary conditions the Kálmán gain can be formulated in terms of the noise covariance matrices [3, 13] and the forward state covariance matrix follows from the solution of the Ricatti equation.

As the output matrix $C_m$ is also modally decoupled, the reconstructed signal in Eq. (11) can be obtained from the Kálmán filter states as

$$\tilde{y}_m[n] = \sum_{j=1}^{n_m} (\hat{x}_{m,j}[n] C_{m,j} + \hat{x}_{m,j^*}[n] C_{m,j^*}) \qquad (14)$$

An example will be introduced in Sec. 5. and is shown in Fig. 2.
<figure><img src="image_2.png" /><figcaption>**Figure 1: Verification Example:** Synthetic correlation functions, modally decomposed by an identified, order-28 state-space model. Only the auto- and cross-correlations between the two reference channels ("201 x", "201 y") are shown. Time lags up to $\tau = 40$ s are displayed and $\delta_j$ indicates the relative contribution of mode j to the input (measured) correlation function.</figcaption></figure>

**Validation Measure and Modal Contribution** A validation measure is developed based on the reconstructed input sequences

$$d[l] = \sum_{j=1}^{n_m} \tilde{d}_j[l] + d^{\epsilon}[l] , \qquad (15)$$

where $\tilde{d}_j[l]$ corresponds to either Eq. (13), (14) or an equivalent formulation for other OMA methods. Here, $l = 1...L$ and $L = \{M, N\}$ is the sequence index (correlation / time), and $d^{\epsilon}[l]$ is an additional error term resulting from an insufficient fit of the identified system. $d[l]$ consists of $p = \{n_1, n_1 n_r\}$ rows that correspond to the signal channels or channel-pairs of the correlation function.

A validation measure of the identified system is then obtained from the normalized correlation coefficient of the input sequences with their synthesized counterpart, averaged over the rows and modes

$$\delta = \sum_{j=1}^{n_m} \delta_j, \quad \text{where} \quad \delta_j = \frac{1}{p} \sum_{p} \frac{\sigma_{d\tilde{d},p}^2}{\sqrt{\sigma_{d\tilde{d},p}^2 \sigma_{dd,p}^2}}. \qquad (16)$$

Here, $\delta_j$ is the contribution of each identified mode to the input data, which is also known as the modal participation factor [5] or the modal contribution [3]. The (co)-variances for each row $p$ of the input and synthesized sequences, respectively, are defined as

$$\sigma_{d\tilde{d},p}^2 = \sum_{l=1}^{L} \tilde{d}_p[l] \tilde{d}_p[l], \quad \sigma_{dd,p}^2 = \sum_{l=1}^{L} d_p[l] d_p[l], \quad \sigma_{d\tilde{d},p}^2 = \sum_{l=1}^{L} d_p[l] \tilde{d}_{j,p}[l] \qquad (17)$$

assuming zero mean sequences.

In order to adapt this measure to the cross-validation approach, the input sequences and their synthesized counterpart are split into $n_{bl}$ blocks along the running index $l$. Subsequently, $\delta$ is computed on the validation subset only, according to the chosen validation estimator, as defined in Sec. 2..
<figure><img src="image_3.png" /><figcaption>**Figure 2: Application Example:** Experimental vibration response signals, modally decomposed by an identified, order-30 state-space model. Only the signals at the two reference channels ("S23 x", "S23 y") are shown. A representative time window of $t = 175$ s is displayed and $\delta_j$ indicates the relative contribution of mode $j$ to the measured signal.</figcaption></figure>

# 5. NUMERICAL VERIFICATION AND EXPERIMENTAL APPLICATION

A guyed steel transmission mast with a height of 200 m has been selected to demonstrate the developed method, as shown in Figure 3. A *verification example* using synthetic acceleration response signals and an *application example* using experimental acceleration response measurements are shown.

The analysis of the acceleration vibration response measurements carried out in 2019 was presented in [14]. Figure 4 shows an extract of the identified modal characteristics.

A simplified numerical model is developed in this study, modeling the mast as a continuous beam with elastic supports, representing the pre-stressed guy cables, and a tuned mass damper, which introduces non-proportional damping; thus the modal parameters must be interpreted as non-classical [15]. Further details with respect to the modeling are described in [14, 16].

Synthetic acceleration response time histories are generated by an inverse Fourier transform of the frequency response function, utilizing generalized modal coefficients [5], multiplied with spectral stochastic wind excitation [17]. The signals are virtually sampled at 5 locations and two directions and varying sample rates and durations. Additional white Gaussian measurement noise is added s. t. SNR = -10 dB.

**Modal Synthetization of Input Sequences**
The following example demonstrates the modal reconstruction of correlation functions and response signals from identified state-space models. The size of the inputs is chosen<fnref n="1" /> as $M = 2048$ and $n_{ord} = 28$ in the covariance-driven case and $p = q = 1024$ and $n_{ord} = 30$ in the data driven case, where $p$ and $q$ are the number of block-rows and -columns of the Block-Hankel matrix and $n_{ord}$ is the order of the identified stochastic state space model. At a duration $T = 3600$ s and sample rate $f_s = 14$ Hz the signal has a total of $N = 50400$ timesteps.

Fig. 1 shows the modal reconstruction of correlation functions for selected physical modes of the veri-

<footnote n="1">¹For the sake of a smooth display of synthesized signals, a much longer correlation function than required for identification purposes was chosen.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Geometry and dimensions of the antenna mast (Photograph: Wolfgang Greiner)</figcaption></figure>

<figure><img src="image_5.png" /></figure>

<table><caption>**Figure 4:** Part of the identified modal parameters of the mast (Note the rotated coordinate systems)</caption><thead><tr><th>Φ<sub>j</sub></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>f<sub>d,j</sub> [Hz]</td><td>0.196</td><td>0.192</td><td>0.471</td><td>0.501</td><td>1.024</td><td>1.043</td><td>1.326</td><td>1.336</td></tr><tr><td>ζ<sub>j</sub> [%]</td><td>0.466</td><td>0.479</td><td>0.615</td><td>0.365</td><td>0.47</td><td>0.25</td><td>0.29</td><td>0.23</td></tr></tbody></table>

fication (numerical) example. Similarly, Fig. 2 shows the modally reconstructed response signal of the application (measurement) example.

The examples demonstrate the effectiveness of the employed reconstruction methods in two scenarios. In the application example, the high overall reconstruction rate indicates low-noise signals. While noise significantly affects the response signal reconstruction, its impact is effectively reduced due to the averaging nature of correlation functions in their reconstruction. The examples also show that the modal contributions, can reliably be used as modal indicators, but are dependent on the excitation and experimental setup, highlighting the need for careful consideration of these factors.

**Full-Factorial Study of the Cross-Validation Measure** To demonstrate the developed cross-validation procedure, the total modal reconstruction indicator $\delta = \sum \delta_j$ from the previous example is used as the validation measure $\gamma$. The input signals are split into $n_{bl} = 40$ blocks and assembled into $k = 10$ randomized subsets of cardinality $n_t = 36$ for identification and $n_v = 4$ for validation. The loss estimate, which is used in the following two full-factorial parameter studies, follows as

$$ \hat{\mathcal{L}}_k = \frac{1}{k} \sum_k \delta_k. \qquad (18) $$

In the case of the experimental application example, the model order $n_{ord} = 4...40$ and the number of block-rows $q = 20...200$ are each sampled at 37 steps. The resulting variations of the loss estimate in Figure 5 indicate, that the validation measure increases with increasing model order and increasing time lags. However, a closer investigation reveals the occurrence of modal splitting at higher model orders. This is not represented in the reconstruction criterion and points to a possible direction of further research.

In the second study, the numerical verification example, the signal duration $T = 600...7200$ s and sample rate $f_s = 3.9...70$ Hz are varied similarly. A stochastic state space model of order $n_{ord} = 28$ and the
<figure><img src="image_6.png" /><figcaption>**Figure 5: Application Example:** Loss estimate $\hat{\mathcal{L}}_k$ for cross-validation of identified stochastic forward innovation state space models (SSI-data) at different model orders $n_{ord}$ and numbers of block-rows/-columns $p, q$.</figcaption></figure>

corresponding validation measure is identified from correlation functions with length $\tau_{max} = 8$ s. The loss estimate is displayed in Figure 6, which indicates, that the sample rate should be chosen about a factor four higher than the highest natural frequency and a longer signal duration is favorable.

The examples show, that the proposed method is consistent over a wide range of parameters and indicates its usefulness for purposes of optimization or uncertainty quantification.

# 6. CONCLUSIONS AND OUTLOOK

In summary, the present study proposes an improved validation criterion for selecting parameter sets in Operational Modal Analysis. It combines input reconstruction with cross-validation techniques to provide a more objective and reliable means for determining the influence of sensitive parameters in the measurement and analysis process.

Notably, the criterion is application- and method-agnostic, scalar and normalized, and robust with respect to overfitting and sample bias. It is designed to be independent of system and excitation dynamics. Extensions of two stochastic subspace identification methods have been developed to include multi-block estimation and reconstruction in a modally decomposed manner, for further use as a stabilization criterion of single modes. The method has been implemented within an object-oriented, open-source Python framework [18].

The effectiveness of this approach has been demonstrated through both numerical and experimental vibration signals. Its application to an experimental example, a 200-meter-tall guyed steel mast, demonstrates its applicability and highlights the critical need for consideration of noise and excitation conditions in the modal reconstruction process. A full-factorial study has further demonstrated the criterion's consistent performance over a broad range of parameters, indicating its usefulness for optimization and uncertainty quantification in modal analysis.
<figure><img src="image_7.png" /><figcaption>**Figure 6: Verification Example:** Loss estimate $\hat{\mathcal{L}}_k$ for cross-validation of stochastic forward state space models (SSI-cov) identified from signals with different durations $T$ and sample rates $f_s$.</figcaption></figure>

It is acknowledged that the method exhibits certain limitations, which offer directions for future research. Specifically, the utilization of other Cross-Validation methods, such as Leave-One-Out (LOO) or Leave-K-Out (LKO), has not been explored, despite being theoretically applicable. The application examples also highlight challenges in the context of modal splitting, indicating potential areas for future research. While the method has been successfully demonstrated for the SSI methods and additionally implemented for the pLSCF method, its extension to other OMA methods, though conceptually possible, remains to be developed and implemented.

## ACKNOWLEDGEMENTS

The financial support of the German Research Foundation (DFG) for the project *Assessment and Reduction of Uncertainties in Operational Modal Analysis* (RUN-OMA) is gratefully acknowledged.

## DATA AVAILABILITY STATEMENT

The implementations developed as part of this study are available under a GPL-3.0 license at https://github.com/pyOMA-dev/pyOMA [18]. The data that support the findings of this study are available upon request from the corresponding author, SM.
REFERENCES

[1] R. S. Pappa, K. B. Elliott, and A. Schenk. "A consistent-mode indicator for the eigensystem realization algorithm". In: *AIAA Dynamics Specialists Conference* (Apr. 16, 1992).

[2] E. Reynders, J. Houbrechts, and G. De Roeck. "Fully automated (operational) modal analysis". In: *Mechanical Systems and Signal Processing* 29 (2012), pp. 228–250.

[3] F. J. Cara et al. "Modal contribution and state space order selection in operational modal analysis". In: *Mechanical Systems and Signal Processing* 38.2 (2013), pp. 276–298.

[4] E. Reynders. "System identification methods for (operational) modal analysis: review and comparison". In: *Archives of Computational Methods in Engineering* (2012), pp. 51–124.

[5] R. Brincker and C. Ventura. *Introduction to Operational Modal Analysis*. Wiley, 2015.

[6] S. Arlot and A. Celisse. "A survey of cross-validation procedures for model selection". In: *Statistics Surveys* 4.0 (2010), pp. 40–79.

[7] B. Peeters and G. de Roeck. "Reference-Based Stochastic-Subspace Identification for Output-Only Modal Analysis". In: *Mechanical Systems and Signal Processing* (1999).

[8] M. Döhler and L. Mevel. "Fast multi-order computation of system matrices in subspace-based system identification". In: *Control Engineering Practice* 20.9 (2012), pp. 882–894.

[9] M. Döhler. "Subspace-based system identification and fault detection: Algorithms for large systems and application to structural vibration analysis". Theses. Université Rennes, 2011.

[10] P. Stoica. *Spectral analysis of signals*. Upper Saddle River, N.J: Prentice Hall, 2005.

[11] P. Van Overschee and B. De Moor. *Subspace Identification for Linear Systems: Theory, Implementation, Applications*. Kluwer Academic Publishers, 1996.

[12] B. Peeters and G. de Roeck. "System identification and damage detection in civil engineering". eng. PhD thesis. KU Leuven, 2000.

[13] R. H. Shumway and D. S. Stoffer. *Time Series Analysis and Its Applications: With R Examples*. Springer International Publishing, 2017.

[14] V. Zabel, S. Marwitz, and A. Habtemariam. "Bestimmung von modalen Parametern seilabgespannter Rohrmasten". In: *Berichte der Fachtagung Baustatik-Baupraxis*. Vol. 14. Baustatik Baupraxis e. V. Mar. 2020, pp. 519–526.

[15] E. Kausel. *Advanced Structural Dynamics*. Cambridge University Press, Aug. 2017.

[16] S. Marwitz, T.Lahmer, and V. Zabel. "Quantification of Polymorphic Uncertainties - A Quasi-Monte Carlo Approach". In: *ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering* (June 1, 2024).

[17] M. Clobes. *Identifikation und Simulation instationärer Übertragung der Windturbulenz im Zeitbereich*. Aachen, 2008.

[18] Simon Marwitz, Volkmar Zabel, and Andrei Udrea. *pyOMA*. 2025.