## 1 Introduction  
The effective health monitoring of critical infrastructure, such as buildings and bridges, is of paramount importance to ensure their longevity, provide high-quality service over their life cycle, and mitigate the potential loss of human lives and economic assets resulting from structural failures, especially during seismic events [1, 2]. The utilization of measured modal parameters, encompassing natural frequencies, modal damping ratios, and mode shapes, has become a widely adopted and straightforward approach in structural health monitoring [3, 4, 5]. These parameters not only facilitate damage detection, but also find applications in tasks such as updating finite element models to match real behavior and implementing vibration control strategies. The seismic response monitoring of critical infrastructure in earthquake-prone regions poses significant challenges. The accurate identification of modal parameters from measured seismic responses is essential for understanding structural behavior and its potential vulnerabilities. A substantial body of the literature has explored the modal parameter identification of structures under seismic excitation [6, 7, 8]. In this context, two primary categories of approaches, time series methods and subspace methods, have been frequently employed in the time domain or time-frequency domain. The time series models include autoregressive models with exogenous input (ARX) and autoregressive moving average models with exogenous input (ARMAX). For instance, [PERSON] and [PERSON] [9] and [PERSON] and [PERSON] [10] employed ARX models to extract modal parameters of buildings through least-squares estimations of the coefficient matrices of ARX models in the time domain. In the time-frequency domain, [PERSON] et al. [11] and [PERSON] and [PERSON] [12] introduced discrete and continuous wavelet transforms into ARX models, respectively. Using an ARMAX model, [PERSON] and [PERSON] [10] and [PERSON] and [PERSON] [13] identified modal parameters of buildings through the estimation of coefficient matrices of ARMAX models using the maximum likelihood method and the differential iteration method in the time domain, respectively, while [PERSON] et al. [14; 15] determined the modal parameters of a structure with incompletely measured base excitations through an iterative multi-stage estimation of coefficient matrices. [PERSON] et al. [16] proposed recursive ARMAX models for identifying instantaneous modal parameters of a time-variant system. Various approaches based on time-varying ARX [17; 18; 19] have been employed to determine instantaneous modal parameters of time-variant or nonlinear systems. Notably, ARMAX modeling is more time consuming than ARX modeling due to the involvement of unmeasured white noise. Subspace methods are classified as deterministic, stochastic, and combined deterministic-stochastic subspace identification techniques [20; 21]. Stochastic subspace identification (SSI) approaches, which include covariance-driven [22; 23] and data-driven [24; 25] SSI algorithms according to the type of data matrices used in the algorithms, have been often employed to process output-only stationary responses, like ambient vibrations. Comprehensive literature reviews on SSI approaches in terms of methodology and engineering applications can be found in the review papers by [PERSON] et al. [7] and [PERSON] et al. [8]. Although earthquake responses, theoretically, do not fit stochastic state-space models, [PERSON] et al. [26] applied a recursive covariance-driven SSI approach to identify the natural frequencies and mode shapes of a tower from its earthquake responses pre-processed using recursive singular spectrum analysis. [PERSON] and [PERSON] [27] presented an improved data-driven SSI approach to process the earthquake responses of structures.