## 1 Introduction  
Many engineering problems are closely related to nonlinear vibrations, some of which can cause irreversible effects on systems. [PERSON], [PERSON], and others introduced nonlinear energy sinks to achieve vibration reduction [1; 2; 3; 4]. However, most vibrations are unavoidable, making the analysis of system vibrations and the prediction of their behavior significantly important for practical engineering evaluation. Meanwhile, multi-scale dynamical systems and the phenomenon of bursting oscillations induced by the multi-scale nature of systems are common in various fields of natural science and engineering. For example, in mechanics, bistable asymmetric laminated composite plates exhibit bursting oscillations under excitation, leading to jump behaviors [5]; in neuroscience, the transmission between nerve cells and endocrine cells at synapses also exhibits bursting oscillations [6]. Therefore, the study of bursting oscillations in this context becomes important and meaningful. A system is said to exhibit bursting oscillations if it has periodic orbits whose time series (at least in one state variable or some coordinate scales) alternate between rapid oscillations and near-stationary behaviors [7; 8]. Currently, the main analysis frameworks for bursting oscillations are based on mathematical modeling and dynamical systems theory. The establishment of models often relies on understanding physical processes, and stability analysis, bifurcation theory, and multi-scale analysis are important components of theoretical research on bursting oscillations. Due to the nonlinear nature of bursting oscillations, theoretical analysis often requires numerical simulations to aid in simulating the dynamics of bursting oscillations, verifying theoretical predictions, and exploring complex states in state space [9; 10]. [PERSON] first proposed the fast-slow analysis method in 1985 [11], which usually involves local linearization and stability analysis at equilibrium points, followed by the division of the system into fast and slow subsystems based on time scales for further study. Although the fast-slow analysis method is a powerful tool for analyzing bursting oscillations, it faces challenges when dealing with highly nonlinear, strongly coupled systems or systems lacking distinct time-scale separation. Meanwhile, in practical applications, commonly used numerical verification methods for bursting oscillations include phase diagram plotting and time history curve plotting. Comprehensive analysis of the complexity of bursting oscillations in practical applications often requires the combination of various numerical methods based on specific research problems. Choosing appropriate numerical methods helps enhance the understanding of the mechanisms of bursting oscillations and improve the efficiency of theoretical research. In recent years, the rapid development of data-driven techniques related to dynamical systems cannot be overlooked. Unlike local linearization around fixed points in stability analysis theory, many data-driven methods focus on finding global linear approximations of systems [12; 13; 14]. This idea of global approximation originates from Koopman operator theory, whose core concept is to transform originally nonlinear dynamical systems into infinite-dimensional linear systems by introducing the Koopman operator, allowing the use of linear system tools for analysis. Although global linearization is very attractive, the extremely high dimensionality of the Koopman operator makes numerical processing difficult. Thanks to the introduction of Takens delay embedding theory [15] and Koopman Mode Decomposition (KMD) theory, as well as the need for rapid modeling in many fields (e.g., rapid flow field reconstruction in industrial applications), many methods for studying the entire dynamical system using univariate measurement data and studying the underlying dynamical system through delay embedding matrices have emerged. In 2014, [PERSON] et al.