A complex ground motion uncertainty quantification method and stochastic simulation system
By screening and classifying ground motion records, constructing a joint probability distribution using an evolved power spectral density model and a Copula function, and implementing a non-Gaussian transformation using a higher-order Hermite polynomial, the accuracy and efficiency issues in simulating complex ground motions in existing technologies have been resolved, resulting in an efficient stochastic simulation and visualization tool.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- DALIAN UNIV OF TECH
- Filing Date
- 2026-01-06
- Publication Date
- 2026-06-02
AI Technical Summary
Existing technologies cannot accurately capture the physical characteristics of earthquakes when simulating complex earthquakes. They ignore nonlinearity and asymmetric correlation, have low computational efficiency, are cumbersome, and are difficult to promote and apply in engineering practice.
By screening measured ground motion records, mainshock-aftershock sequence type and near-fault pulse ground motion are identified. The joint probability distribution of model parameters is constructed using the evolved power spectral density model and Copula function. Non-Gaussian transformation is achieved using a high-order Hermite polynomial model. Non-Gaussian ground motion acceleration time history is generated by combining the spectral representation of intrinsic orthogonal decomposition and integrated into a visualization simulation platform.
It improves the realism and physical meaning of stochastic simulation results, enhances computational efficiency, lowers the technical threshold, and provides structural engineers with a convenient and efficient tool for seismic performance assessment.
Smart Images

Figure CN122133302A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of earthquake engineering, specifically to a method for quantifying the uncertainty of complex ground motions and a stochastic simulation system. Background Technology
[0002] In performance-based earthquake engineering design and assessment, nonlinear dynamic time history analysis is a core tool for evaluating the seismic performance of structures. This analytical method requires a large number of acceleration time history samples that accurately reflect the seismic hazard of the target site as input. Due to the scarcity of strong earthquake records, generating artificial ground motion time histories through stochastic simulation methods has become a key technology to meet engineering requirements.
[0003] Existing methods for generating artificial ground motions typically aim to match the target response spectrum and consider the non-stationary characteristics of ground motion intensity and spectrum changes over time. However, as research progresses, these methods still have significant shortcomings in accurately quantifying and simulating the uncertainties of complex ground motions.
[0004] First, existing techniques for statistical modeling earthquake motion often treat earthquake records with different physical causes (such as ordinary far-field ground motions, near-fault pulse ground motions, and mainshock-aftershock sequences) as a single stochastic process. This results in models failing to accurately capture the unique physical characteristics of specific types of earthquakes. Even when statistical analysis is performed on model parameters, simple linear correlation coefficient matrices are mostly used to describe the dependencies between parameters. This approach ignores the actual, complex nonlinear and asymmetric correlations, reducing the accuracy of the physical meaning of stochastic models.
[0005] Secondly, current simulation methods need improvement in terms of computational efficiency and statistical reproducibility. Traditional Monte Carlo stochastic simulations require generating massive amounts of samples to stably cover the probability space of parameters, resulting in high computational costs. Furthermore, most methods assume that ground motion processes follow a Gaussian distribution, while extensive experimental data show that ground motion acceleration time histories generally exhibit non-Gaussian characteristics such as sharp peaks and thick tails. This characteristic is crucial for assessing the cumulative damage effects on structures, but existing methods often ignore this or introduce non-Gaussian characteristics in complex and inefficient ways.
[0006] Finally, although the academic community has developed a variety of advanced stochastic simulation algorithms, these algorithms are often cumbersome, involving multiple independent professional stages such as data processing, parameter identification, probabilistic modeling, and randomization. This results in a fragmented simulation process, requiring a strong theoretical background and programming skills for implementation, which constitutes a high technical barrier for ordinary engineering technicians and limits the promotion and application of advanced simulation technologies in engineering practice. Summary of the Invention
[0007] To address the shortcomings of existing technologies, this invention provides a method for quantifying the uncertainty of complex ground motions and a stochastic simulation system. This solves the problems of existing technologies in simulating complex ground motions, such as distortion of physical characteristics and parameter correlations due to model simplification, low efficiency of stochastic simulation and inability to reproduce key non-Gaussian statistical features, and cumbersome and fragmented overall processes that are difficult to promote and apply in engineering practice.
[0008] To achieve the above objectives, the present invention provides a method for quantifying the uncertainty of complex seismic motions, the method comprising: Measured ground motion records are screened from the earthquake engineering database to identify mainshock-aftershock sequence-type ground motions and near-fault pulse ground motions, and the measured ground motion records are grouped according to the mean shear wave velocity. An evolutionary power spectral density model was adopted as a unified theoretical framework, and model parameters were identified based on measured ground motion records. Statistical analysis is performed on the identified model parameters to fit the optimal marginal probability distribution, and the correlation structure between the model parameters is constructed based on Copula function theory to establish a joint probability distribution model of the model parameters; A higher-order Hermite polynomial model for implementing non-Gaussian transformation is established, and the mapping relationship between the autocorrelation function of the standard Gaussian process and the autocorrelation function of the non-Gaussian process is established through the higher-order Hermite polynomial model. Random parameter samples are generated from the joint probability distribution model using a representative point set selection strategy. The random parameter samples are substituted into the evolution power spectral density model to determine the target evolution power spectral density. A Gaussian random process is generated by the spectral representation method based on intrinsic orthogonal decomposition. The Gaussian random process is then converted into a non-Gaussian ground motion acceleration time history sample using the mapping relationship. The response spectrum of the non-Gaussian ground motion acceleration time history sample is calculated and compared with the response spectrum of the measured ground motion record to verify the engineering applicability of the random simulation sample. The above steps are then integrated into a visual random ground motion simulation platform.
[0009] Furthermore, the evolved power spectral density model consists of a stationary power spectral density function and a time-frequency modulation function; the stationary power spectral density function is represented by the Clough-Penzien spectrum and includes parameters such as the site soil damping ratio, bedrock damping ratio, site soil dominant circular frequency, and bedrock dominant circular frequency; the time-frequency modulation function includes time parameters that control the attenuation of ground motion energy and the arrival time of the peak.
[0010] Furthermore, the process of identifying the parameters of the evolving power spectral density model based on measured ground motion records includes: decomposing the acceleration time history into wavelet packet coefficients of different times and frequencies using wavelet packet transform; reconstructing the wavelet packet coefficients into sub-signals of different frequency bands using inverse wavelet packet transform; calculating the energy of the sub-signals according to the definition of evolving power spectral density to obtain the energy distribution function estimated based on wavelet packet transform; and fitting the measured ground motion record using a genetic algorithm based on the principle of equal energy in the frequency domain with the goal of approximating the energy distribution function, thereby identifying the model parameters in the evolving power spectral density model.
[0011] Furthermore, the process of fitting the optimal marginal probability distribution of the model parameters includes: using the KS test method and the Akaike information content criterion to screen the optimal marginal probability distribution type of the model parameters from the candidate distribution models; the process of constructing the correlation structure between the model parameters includes: The joint probability distribution is described using a Copula function, and the optimal Copula function type is determined by the Akaike information content criterion. The correlation structure is then constructed based on the optimal Copula function type. When the measured ground motion record is a mainshock-aftershock sequence type ground motion, a high-dimensional correlation between the parameters of the mainshock evolution power spectral density model is established using a Copula structure, and the cross-correlation between the mainshock model parameters and the aftershock model parameters is established using a binary Copula function.
[0012] Furthermore, the process of establishing a high-order Hermite polynomial model for realizing the non-Gaussian transformation includes: setting the skewness and kurtosis statistical characteristics of the target non-Gaussian ground motion acceleration time history; calculating the Hermite shape coefficients based on the skewness and kurtosis statistical characteristics using linear moments; and constructing a second-order linear system containing the Hermite shape coefficients and polynomial coefficients. By solving the second-order linear system, the transformation relationship is determined, and the analytical mapping model between the autocorrelation function of the standard Gaussian process and the autocorrelation function of the non-Gaussian process is derived.
[0013] Furthermore, the specific steps for generating the Gaussian random process include: constructing a cross-spectral density matrix using the target evolutionary power spectral density generated by the evolutionary power spectral density model; performing eigenorthogonal decomposition on the cross-spectral density matrix to obtain eigenvalues and eigenvectors; and combining the eigenvalues, the eigenvectors, and the random phase angle to synthesize the Gaussian random process using spectral representation.
[0014] Furthermore, when the measured ground motion record is a mainshock-aftershock sequence type ground motion, the evolution power spectral density model includes the mainshock evolution power spectral density model and the aftershock evolution power spectral density model. The mainshock evolution power spectral density model and the aftershock evolution power spectral density model correspond to the independently set mainshock duration and aftershock duration, respectively, and a time interval is set between the mainshock and the aftershock; the model parameters include the dominant circular frequency of the site soil, the site soil damping ratio, the spectral intensity factor, and the time-frequency modulation parameters for the mainshock and the aftershock, respectively.
[0015] Furthermore, when the measured ground motion record is a near-fault pulse ground motion, the method includes decomposing the near-fault pulse ground motion into low-frequency pulse components and high-frequency residual components for separate simulation: using a pulse identification method based on continuous wavelet transform to extract the velocity time history of the strongest pulse direction; The velocity-time history of the strongest pulse direction is fitted using the Gabor pulse function model to identify pulse model parameters including peak pulse velocity, pulse occurrence time, pulse period, pulse half-wave number, and pulse phase angle, which are then used to generate the low-frequency pulse component. The residual acceleration time history is obtained by subtracting the fitted pulse component from the original record. The parameters of the residual acceleration time history are identified using the evolved power spectral density model to generate the high-frequency residual component. The generated low-frequency pulse component and the high-frequency residual component are superimposed to synthesize near-fault pulse ground motion.
[0016] A second aspect of the present invention provides a stochastic simulation system for quantifying the uncertainty of complex ground motion, the system comprising: The data filtering module is used to filter measured ground motion records from the database, identify non-stationary ground motions, mainshock-aftershock sequence-type ground motions, and near-fault pulse ground motions, and group them accordingly. The parameter identification module is used to identify the parameters of the evolving power spectral density model based on measured ground motion records, and to identify the pulse function model parameters for near-fault pulse ground motions. The probability modeling module is used to fit the optimal marginal probability distribution of the parameters and construct the correlation structure using the Copula function to establish the joint probability distribution model of the model parameters; The non-Gaussian model building module is used to establish higher-order Hermite multinomial models and autocorrelation function mapping relationships; The random simulation module is used to generate random parameter samples from the joint probability distribution model using a representative point set selection strategy, and to generate non-Gaussian ground motion acceleration time history samples using a spectral representation method based on intrinsic orthogonal decomposition. The visualization and interaction module is used to provide a user interface and display simulation results and verification comparison charts.
[0017] This invention provides a method for quantifying the uncertainty of complex seismic motions and a stochastic simulation system. It offers the following advantages: 1. This invention, through quantitative classification and identification of measured ground motion records in step S1, separates complex ground motion phenomena from ordinary records. Furthermore, it utilizes high-dimensional Copula theory to establish joint probability distribution models for the model parameters of different types of ground motions. This avoids the feature smoothing problem caused by mixing all ground motion records for modeling, and accurately captures and quantifies the unique physical characteristics of various ground motions and the complex nonlinear correlations between their parameters, thereby improving the realism and accuracy of the physical meaning of the stochastic simulation results.
[0018] 2. In step S4, this invention introduces a nonlinear transformation model based on high-order Hermite polynomials and estimates the model coefficients using the L-moment method, enabling the synthesized ground motion acceleration time history to reproduce the non-Gaussian statistical characteristics (such as spikes and skewness) commonly found in real records. Simultaneously, in step S5, the synthesis process employs a spectral representation method based on intrinsic orthogonal decomposition (POD) and a representative point set selection strategy. Compared to traditional Monte Carlo stochastic simulation methods, this approach reduces the required number of samples and computation time while maintaining simulation accuracy, thus improving the computational efficiency of generating large-scale stochastic ground motion samples.
[0019] 3. In step S7, this invention integrates the entire process from data screening, parameter identification, probabilistic modeling to stochastic synthesis and verification into a visualized stochastic simulation platform. This platform encapsulates complex stochastic dynamics theory and numerical algorithms within a graphical user interface. Users can complete stochastic simulations of specific sites and seismic motion types simply by configuring parameters, and can intuitively view and export analysis results that meet engineering application standards. This lowers the barrier to entry for advanced stochastic simulation technology, providing structural engineers with a convenient and efficient tool for seismic performance assessment and risk analysis. Attached Figure Description
[0020] Figure 1 This is a schematic diagram of the complex seismic uncertainty quantification process of the present invention; Figure 2 EPSD model parameters of the present invention Edge diagram; Figure 3 EPSD model parameters of the present invention Edge diagram; Figure 4 EPSD model parameters of the present invention Edge diagram; Figure 5 This is a schematic diagram of the system architecture of the present invention. Detailed Implementation
[0021] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0022] See attached document Figure 1 This invention provides a method for quantifying the uncertainty of complex ground motion, comprising the following steps: S1. Screen the measured ground motion records from the earthquake engineering center database, identify non-stationary ground motions, mainshock-aftershock sequence-type ground motions and near-fault pulse ground motions in the measured ground motion records, and group the measured ground motion records according to the average shear wave velocity. S2 adopts the evolutionary power spectral density model as a unified theoretical framework and identifies the model parameters of the evolutionary power spectral density model based on measured ground motion records; S3. Statistical analysis is performed on the identified model parameters to fit the optimal marginal probability distribution of the model parameters, and the correlation structure between the model parameters is constructed based on the Copula function theory, thereby establishing a joint probability distribution model of the model parameters. S4. Establish a high-order Hermite polynomial model to realize the non-Gaussian transform, and establish the mapping relationship between the autocorrelation function of the standard Gaussian process and the autocorrelation function of the non-Gaussian process through the high-order Hermite polynomial model. S5. Random parameter samples are generated from the joint probability distribution model using a representative point set selection strategy. The random parameter samples are substituted into the evolution power spectral density model to determine the target evolution power spectral density. A Gaussian random process is generated by the spectral representation method based on intrinsic orthogonal decomposition. The Gaussian random process is then converted into a non-Gaussian ground motion acceleration time history sample using a mapping relationship. S6. Calculate the response spectrum of the non-Gaussian ground motion acceleration time history sample and compare it with the response spectrum of the measured ground motion record to verify the engineering applicability of the random simulation sample. S7 is a visualization-based stochastic seismic simulation platform developed using Matlab GUI, integrating the workflow into the visualization-based stochastic seismic simulation platform.
[0023] The core algorithm principles and implementation process involved in steps S1 to S7 above will be explained in detail below with specific mathematical models and examples.
[0024] In step S1, the strategy for screening, identifying, and grouping ground motion data specifically includes the following sub-steps: Step S11: Acquisition and Preliminary Screening of Ground Motion Data. The NGA-West2 strong earthquake database from the Pacific Earthquake Engineering Center (PEER) is used as the data source. Records from stations installed on dams, bridges, or other structures are discarded, retaining only ground motion records from the free field or reference site to eliminate interference from structural interactions on ground motion characteristics. Data preprocessing, including denoising, baseline correction, and filtering, can be achieved using conventional signal processing techniques known to those skilled in the art and will not be elaborated upon here.
[0025] Step S12: Quantitative Identification of Near-Fault Pulse Ground Motion. For the selected records, the pulse characteristic determination logic is executed. The identification of pulse ground motion does not rely on a single subjective human judgment, but uses quantitative energy indicators.
[0026] Specifically, a continuous wavelet transform is performed on the velocity time history of each ground motion record, and the signal component corresponding to the wavelet coefficient with the highest energy is extracted as a potential pulse. The ratio of the total energy of the original record to the energy of the extracted potential pulse component is calculated, and combined with the peak ground velocity (PGV) index for determination. When the proportion of the extracted pulse component energy to the total energy exceeds a preset threshold (e.g., 30%), and the PGV is greater than a preset limit (e.g., 30 cm / s), the record is marked as a near-fault pulse ground motion. Records that fail to meet the above quantitative indicators are temporarily marked as non-pulse ground motions.
[0027] Step S13: Association Identification of Mainshock-Aftershock Sequences. Establish an association mapping between mainshocks and aftershocks in the database. Perform spatiotemporal matching based on station ID and event ID. The identification logic is as follows: multiple earthquake events recorded at the same station that are temporally continuous and belong to the same geological tectonic activity. The event with the largest magnitude and earliest occurrence is defined as the mainshock, and subsequent events are defined as aftershocks. For each pair of mainshock-aftershock sequences, extract the corresponding mainshock acceleration time history and aftershock acceleration time history, and record the time interval between them to form a mainshock-aftershock sequence dataset. Records that do not meet the sequence characteristics and the pulse characteristics in step S12 are classified as ordinary non-stationary ground motions.
[0028] Step S14: Site grouping based on shear wave velocity. To quantify the influence of site conditions on the seismic motion spectrum characteristics, the average shear wave velocity within 30 meters underground at the station location is used for grouping. The identified ground motion records (ordinary non-stationary and mainshock / aftershock sequences) were grouped. The grouping criteria followed the differences in engineering geological characteristics, dividing the measured ground motion records into statistical units with similar geological backgrounds. This provides a physical basis for subsequent steps involving statistical characteristic analysis of model parameters and fitting of marginal probability distributions under different site conditions. All filtered, identified, and grouped ground motion acceleration time history data and their metadata (including...) , , All parameters (e.g., 'etc.') are formatted and stored as input data for the subsequent parameter recognition module.
[0029] See attached document Figure 2 - Appendix Figure 4 As a specific embodiment of the present invention, after screening and grouping in step S1, the number of measured non-pulse and non-sequence ground motions corresponding to each site level is shown in Table 1.
[0030] Table 1. Number of measured ground motions corresponding to different field levels. ; In step S2, the parameter identification of the evolved power spectral density model based on wavelet packet transform specifically includes the following sub-steps: Step S21: Constructing the Evolving Power Spectral Density Model. The evolving power spectral density model is used as a unified physical model to describe the stochastic process of ground motion. This model defines the one-sided power spectral density (PSD) function of non-stationary ground motions. Time-frequency modulation function The one-sided power spectral density function of steady ground motion The product of these two products is expressed mathematically as follows: ; in, It is a time-frequency modulation function; and The unilateral PSDs are respectively for steady and non-steady ground motions; Represents a time variable; This represents the circumfrequency variable.
[0031] To accurately reflect the frequency domain characteristics of earthquake motion, the stationary power spectral density function... The Clough-Penzien spectral model was selected, which comprehensively considers the bedrock filtering effect and the amplification effect of local soil layers. Its specific expression is as follows: ; in, and These are the damping ratios of the site soil and bedrock, respectively; and These are the dominant circular frequencies of the site soil and bedrock, respectively. The spectral intensity factor is for non-stationary ground motion.
[0032] Spectral intensity factor Peak ground acceleration Peak factor And there are analytical relationships between the site soil parameters: ; in, Peak factor; This represents the average peak ground acceleration (PGA).
[0033] To describe the non-stationary changes in earthquake intensity and frequency over time, a time-frequency modulation function is used. The following exponential function form is adopted: ; in, These are parameters that control the attenuation of seismic energy and the time of peak arrival. It is an exponential function; The normalized time parameter, i.e., the time point when the time-frequency modulation function reaches its peak, is calculated as follows: ; in, It is the natural logarithm.
[0034] Based on the above definition, the theoretical energy distribution function of the evolution power spectral density model is... Through the Integrating in the time domain yields: ; Step S22: Estimation of Measured Energy Distribution Based on Wavelet Packet Transform. The measured ground motion records are processed using wavelet packet transform to extract their energy distribution characteristics in the time-frequency domain. The measured acceleration time history is then processed... Perform wavelet packet decomposition and calculate the th Layer Wavelet packet coefficients of each node : ; in, express Scale and translation parameters Next Wavelet packet coefficients corresponding to the frequency band; This represents integration operations over the entire time domain; The wavelet basis function is used; in this embodiment, the Meyer wavelet is employed.
[0035] By using inverse wavelet packet transform, the decomposed coefficients are reconstructed into a series of sub-signals with different frequency bands. : ; in, For the frequency range in Sub-signals between.
[0036] Based on the principle of energy conservation and using Parseval's theorem, calculate the... Energy of each frequency band sub-signal : ; Further calculation of the time-varying power spectral density of the sub-signal : ; in, The duration representing the acceleration time history. For the first The bandwidth of each frequency band.
[0037] Finally, the frequency domain energy distribution function estimated based on measured records is obtained. : ; Step S23: Model Parameter Identification and Optimization. Construct the objective function for parameter identification, aiming to minimize the theoretical energy distribution function. Compared with the measured energy distribution function The difference between them. The objective function is defined as the mean square error of both in the full frequency domain: ; in, It represents an infinitesimal frequency interval.
[0038] A genetic algorithm is used to optimize the objective function described above. Algorithm parameters such as population size, crossover probability, and mutation probability are set, and the optimal model parameter vector is obtained after multiple generations of iterative convergence. This completes the parameter identification of a single ground motion record.
[0039] Step S24: Special Processing of Near-Fault Pulse Ground Motion Records. For ground motion records identified as near-fault pulse type in Step S1, a decomposition strategy is adopted. Velocity pulse components are extracted using continuous wavelet transform. To accurately describe the pulse characteristics, a Gabor pulse function model is used to quantitatively fit the extracted velocity pulses. The velocity time history expression of this model is... for: ; in, These represent the peak pulse velocity, pulse occurrence time, pulse period, pulse half-wave number, and pulse phase angle, respectively.
[0040] To identify the above optimal parameter set Construct the following least-squares objective function, and minimize the extracted pulse. With theoretical models Solve by the energy difference between them: ; The low-frequency pulse component obtained from the fitting is subtracted from the original acceleration record to obtain the high-frequency residual component. For this high-frequency residual component, steps S21 to S23 are repeated to identify its corresponding model parameters using the evolved power spectral density model.
[0041] The low-frequency pulse component obtained from the fitting is subtracted from the original acceleration record to obtain the high-frequency residual component. For this high-frequency residual component, steps S21 to S23 are repeated to identify its corresponding model parameters using the evolved power spectral density model. Finally, the near-fault pulse ground motion is described by a set of pulse parameters and a set of EPSD model parameters.
[0042] Step S25: Parameter Identification of Mainshock-Aftershock Sequence Ground Motion. For mainshock-aftershock sequence ground motion, steps S21 to S23 are performed independently for both the mainshock and aftershock records. The identification results include the EPSD model parameter sets for the mainshock and the EPSD model parameter sets for the aftershocks, as well as the time interval between them. The mainshock and aftershocks have the same model structure, but the parameter values reflect their respective independent spectrum and non-stationary characteristics.
[0043] In step S3, the joint probability distribution of the model parameters is modeled based on Copula theory, which specifically includes the following sub-steps: Step S31: Marginal Probability Distribution Fitting of Model Parameters. For the massive model parameter sample set identified in Step S2, it is necessary to determine the statistical characteristics of each model parameter. For each parameter component in the sample set (e.g., ... , , , , , Statistical inference is performed separately for each distribution. A candidate probability distribution model library is constructed, including normal, log-normal, Gamma, Weibull, and generalized extreme value distributions. The parameters of each candidate distribution are estimated using maximum likelihood estimation (MLE), and the values of the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) are calculated. Based on the principle of minimizing AIC or BIC, the marginal probability density function with the best fit is selected from the candidate library. and the corresponding cumulative distribution function ,in Representing the Each model parameter The value of this parameter is given. Conventional methods for testing the statistical characteristics of parameters (such as the KS test) are well-known techniques in this field and will not be elaborated upon here.
[0044] See attached document Figure 2 - Appendix Figure 4 As a specific implementation example, the optimal marginal probability distribution fitting results of some EPSD model parameters for the selected non-stationary ground motions are shown in Table 2.
[0045] Table 2 Optimal marginal probability distribution of EPSD model parameters for non-stationary ground motion. ; For the special ground motion types identified in step S1, the same marginal probability distribution fitting and subsequent Copula modeling are performed. For example, for near-fault pulse ground motions, statistical modeling is performed on the pulse model parameters and residual component EPSD model parameters, and an example of the optimal marginal probability distribution fitting results is shown in Table 3.
[0046] Table 3. Parameter distribution types and distribution coefficients of the near-fault pulse ground motion stochastic model ; Step S32: Construction of Joint Probability Distribution Model Based on Copula Functions Due to the complex nonlinear correlations among the parameters of the seismic motion model (e.g., the correlation between site frequency and damping ratio), relying solely on marginal distributions cannot accurately describe the joint probability behavior of the parameter set. Copula function theory is introduced, and based on Sklar's theorem, multiple marginal distribution functions are connected through Copula functions to construct a joint distribution. For any two model parameters... and Its joint cumulative distribution function Represented as: ; in, Let Copula be the cumulative distribution function. and Parameters and The marginal cumulative distribution function.
[0047] Correspondingly, the joint probability density function Its analytical form is obtained by taking the partial derivative of the joint cumulative distribution function: ; in, The copula density function reflects the correlation structure among variables; and These are the marginal probability density functions.
[0048] Regarding the selection of Copula functions, a function library is established that includes elliptic Copulas (such as Gaussian and t-Copula) and Archimedean Copulas (such as Clayton, Gumbel, and Frank). The Euclidean distance or AIC value between empirical Copulas and various theoretical Copulas is calculated, and the Copula function form with the best fit and its related parameters (such as correlation coefficient) are selected. Or generate meta-parameters This allows for the quantification of the dependency structure between parameters.
[0049] See attached document Figure 2 - Appendix Figure 4 Step S33: High-dimensional Copula modeling of the mainshock-aftershock sequence. For mainshock-aftershock sequence-type ground motions, the random variables involved have high dimensionality (mainshock parameter set + aftershock parameter set + time interval, exceeding 10 dimensions), and there is asymmetric tail correlation between the variables. A Copula (D-vine Copula) structure is used to decompose and model this high-dimensional dependency.
[0050] As a specific embodiment, the number of measured main and aftershock sequences after being screened, identified and grouped in step S1 is shown in Table 4.
[0051] Table 4. Number of measured mainshock and aftershock sequences corresponding to different field levels ; For the above mainshock and aftershock sequence data, first execute step S31 to identify the optimal marginal probability distribution of each parameter of the mainshock (M) and aftershock (A), and the results are shown in Table 5.
[0052] Table 5. Optimal marginal probability distribution of EPSD model parameters for the mainshock and aftershock sequence. ; Subsequently, the correlation structure of the aforementioned high-dimensional parameters is modeled based on Copula. Taking a 4-dimensional random variable as an example, the joint probability density function based on the D-Vine structure is used. The total decomposition form is shown below: ; in, This represents the joint probability distribution function of a 4-dimensional random variable. and For random variables The cumulative distribution function (CDF) and probability density function (PDF). In the given Regarding the situation The conditional distribution function. for and The edge Copula density function between them. Represents a conditional random variable and The bivariate Copula density function between.
[0053] By constructing such a tree structure, the optimal vine structure is determined using the maximum spanning tree algorithm, and the parameters of each pair-copula are estimated in hierarchical order. An example of parameter information for a specific D-Vine Copula model is shown in Table 6.
[0054] Table 6. Parameter information for D-vine Copula ; Step S34: Derivation of the Conditional Probability Function for Random Sampling To support the generation of random parameters in the subsequent step S5, it is necessary to derive the Copula-based conditional cumulative distribution function (i.e., the h-function). For known variables... Variables under conditions Conditional cumulative distribution function The calculation formula is as follows: ; This formula establishes a mathematical pathway for deriving the conditional probability of another related variable from the marginal probability of a known variable, forming the basis for implementing the Rosenblatt transform or inverse transform sampling algorithm. For high-dimensional Vine Copula, the calculation of this conditional probability function is achieved by recursively calling the aforementioned binary conditional probability formula.
[0055] In step S4, the non-Gaussian mapping relation based on the higher-order Hermite polynomial is constructed, specifically including the following sub-steps: Step S41: Constructing a Non-Gaussian Transformation Model. To accurately simulate the prevalent non-Gaussian properties (i.e., the phenomenon of skewness and kurtosis deviating from a normal distribution) in seismic acceleration records, a memoryless nonlinear transformation model based on high-order Hermite polynomials is established. This model transforms the underlying standard Gaussian stochastic process... Mapped to a non-Gaussian stochastic process with a target marginal probability distribution Considering that seismic ground motion records often exhibit strong soft-hardening characteristics, a third-order polynomial alone is insufficient to describe the characteristics of higher-order moments. Therefore, a fifth-order Hermite polynomial model is adopted, the specific mathematical expression of which is: ; in, It is a stochastic process with a symmetric non-Gaussian marginal probability distribution function. This is the proportionality coefficient. This represents a standardized Gaussian random process. and This represents the Hermite shape factor, which can be obtained through linear moments.
[0056] scaling factor With coefficient , The following analytic relationship exists between them: ; Step S42: Solving for coefficients based on L-moments To determine the coefficients in the above model and This study utilizes L-moments estimation to replace the traditional central moment method, improving robustness to outliers in ground motion data. The first six L-moments are calculated based on measured ground motion data. The calculation of L-moments involves probability-weighted moments; the specific L-moments... , and By analyzing the observation data and its cumulative distribution function Points earned: ; ; ; Based on the calculated L-moment, a dimensionless L-skewness is defined. and L-kurtosis : ; ; Establish , With Hermite polynomial coefficients , The system of nonlinear mapping equations between them. The approximate linear system of equations obtained through least-squares fitting is as follows, used to directly solve for the coefficients: ; ; Step S43: Mapping Transformation of the Autocorrelation Function. Since the nonlinear transformation alters the second-order statistical characteristics (power spectral density or autocorrelation function) of the stochastic process, to ensure that the generated non-Gaussian ground motion has the correct target power spectral density, it is necessary to establish the non-Gaussian process autocorrelation function. Autocorrelation function of underlying Gaussian process The mapping relationship between them. Based on the orthogonality of Hermite polynomials, the analytical relation between them is derived: ; in, It represents the mathematical expectation.
[0057] Expanding the above equation, we obtain the autocorrelation function of the underlying Gaussian process. The fifth-degree polynomial equation, i.e., the mapping function : ; ; ; in, and These are the autocorrelation functions of a standard non-Gaussian process and a Gaussian process, respectively. and Represents the coefficients of a polynomial.
[0058] To simplify calculations, intermediate variables Defined as: ; By solving the above fifth-order equations, the autocorrelation function of the target non-Gaussian ground motion can be obtained. The autocorrelation function required to back-calculate the underlying Gaussian process The revised version This will serve as the input target for subsequent Gaussian process simulations, thereby ensuring that the final non-Gaussian ground motion after formula transformation not only satisfies the target edge distribution but also accurately satisfies the target power spectral density characteristics.
[0059] In step S5, the random ground motion synthesis combining the POD spectral representation and mapping transformation specifically includes the following sub-steps: Step S51: Random Parameter Generation Based on Representative Point Sets To maximize the coverage of the probability distribution characteristics of the parameter space with a limited number of samples, the traditional Monte Carlo random sampling method is abandoned, and a representative point set selection strategy based on iterative screening and rearrangement is adopted. Based on the joint probability distribution model established in Step S3, the multidimensional parameter space is discretized. Using F-discrepancy as an evaluation index of the uniformity of the point set distribution, this deviation is minimized to determine the optimal parameter combination sample set.
[0060] Specifically, the conditional cumulative distribution function derived in step S34 is transformed using the Rosenblatt inverse transform. This is applied to low-bias sequences (such as Sobol or Halton sequences) generated in hypercube space, thereby mapping them to ground motion model parameter vectors that conform to actual physical correlations. ,in , This represents the total number of simulated samples.
[0061] Step S52: Determining the target evolution power spectral density for each generated set of model parameter vectors. Substitute the values into the evolutionary power spectral density model in step S21 to calculate and determine the corresponding target evolutionary power spectral density function. This function uniquely determines the first... The second-order statistical characteristics of simulated ground motion in the time-frequency domain.
[0062] Step S53: Efficient Gaussian Process Simulation Based on POD To address the low computational efficiency of traditional spectral representation methods when simulating large numbers of samples, Proper Orthogonal Decomposition (POD) technology is introduced. The power spectral density function of the target is then analyzed. Discretize the matrix to construct a time-frequency matrix. Perform eigenvalue decomposition (or singular value decomposition) on this matrix to extract several principal eigenvalues. and its corresponding orthogonal eigenvectors (i.e., intrinsic modes). ,in , To truncate the number of modes, dimensionality reduction is achieved by truncating modes corresponding to high-order small eigenvalues. The Gaussian random process is then reconstructed using the dimensionality-reduced modes. The formula for the spectral representation based on POD is as follows: ; in, This indicates the operation of taking the real part; These are mutually independent random phase angle variables that follow a standard normal distribution; This represents the frequency step size. It is the characteristic function; The number of main modes to retain.
[0063] The Fast Fourier Transform (FFT) algorithm is used to accelerate the summation operation in the above formula, thereby efficiently generating standard Gaussian random process samples with the target power spectrum characteristics. .
[0064] Step S54: Synthesis of non-Gaussian ground motion time histories. The standard Gaussian random process samples generated in step S53 are synthesized. As input, before the transformation, the target power spectrum of the input Gaussian process is pre-corrected according to the autocorrelation function mapping relationship derived in step S43 to compensate for the spectral distortion caused by the nonlinear transformation. Then, the pre-corrected Gaussian process is substituted into the fifth-order Hermite polynomial nonlinear transformation model defined in step S41 to obtain the final synthesized non-Gaussian ground motion acceleration time history. After this transformation, the output time series It has: (1) Evolutionary power spectral density characteristics consistent with measured records; (2) The statistical characteristics of the model parameters conform to the joint probability distribution of the target; (3) True non-Gaussian (skewness and kurtosis) statistical properties.
[0065] For the simulation of near-fault pulse ground motion, the above-mentioned high-frequency residual components are generated. Based on this, pulse parameters are randomly generated according to the probability distribution of pulse parameters identified in step S24, low-frequency velocity pulses are constructed and differentially obtained to obtain acceleration pulses, and finally the two are superimposed in the time domain to synthesize a complete near-fault pulse ground motion time history.
[0066] In step S6, the reaction spectrum verification and applicability assessment of the simulated samples specifically includes the following sub-steps: Step S61: Seismic Response Calculation of a Single-Degree-of-Freedom System To evaluate the engineering characteristics of the generated non-Gaussian ground motion acceleration time history samples, the corresponding elastic response spectra are calculated. A series of systems with different natural periods are established. (or angular frequency) A single-degree-of-freedom (SDOF) linear elastic system. For each generated non-Gaussian ground motion acceleration time history sample... Using this as the input excitation, solve the differential equations of motion for the SDOF system: ; in, , and These represent the displacement, velocity, and acceleration responses of the SDOF system relative to the ground, respectively. The damping ratio is 0.05 (i.e., 5%), which is the standard engineering damping ratio in this embodiment.
[0067] For solving the aforementioned second-order linear differential equation, those skilled in the art can use the Newmark-β method, Wilson-θ method, or Duhamel integration method for numerical integration. The specific numerical iteration schemes are well-known techniques in the field and will not be elaborated upon here. The absolute acceleration response peak, i.e., the spectral acceleration, is obtained by solving the equation. : ; Traverse the preset period range (e.g., 0.01 seconds to 10 seconds) to obtain the complete acceleration response spectrum curve.
[0068] Step S62: Statistical Characteristic Analysis of Response Spectrum Considering that seismic response spectra typically follow a log-normal distribution, logarithmic domain statistical analysis is performed on the response spectra of both the simulated sample set and the measured record set. For those containing... The simulated sample set of time histories, at any period point At this point, calculate its logarithmic mean spectral value. and logarithmic standard deviation : ; ; The logarithmic mean was converted to a geometric mean response spectrum using an exponential transformation, and the corresponding plus or minus one standard deviation was calculated. Confidence intervals are used to characterize the dispersion of the simulated samples. Similarly, the same statistical calculations are performed on the measured ground motion record set selected in step S1 to obtain the target statistical spectrum curve.
[0069] Step S63: Consistency Verification of Simulated and Target Spectra. The statistical response spectrum (mean curve and standard deviation envelope) of the simulated sample is compared with the measured target statistical response spectrum in the same coordinate system. A relative error index is defined over the entire period. : ; in, The number of discrete periodic points is used. The ability of the evolved power spectral density model and the joint probability distribution model of parameters in steps S2 and S3 to capture the intensity and spectral characteristics of ground motion is verified by checking the consistency of the mean spectrum; the ability of the Copula correlation structure model to reproduce the random variability of ground motion is verified by checking the consistency of the standard deviation spectrum. If the relative error index is less than a preset threshold (e.g., 5%), and the simulated spectrum curve fluctuates within the confidence interval of the target spectrum curve, the simulated sample is deemed to have passed the response spectrum consistency test.
[0070] Step S64: Verification of Non-Gaussian Statistical Characteristics In addition to verifying the frequency domain and response spectrum characteristics, the non-Gaussian characteristics of the generated samples also need to be verified in the time domain. Calculate the skewness and kurtosis coefficients of the acceleration amplitudes of all simulated samples, and obtain their statistical mean. Compare this statistical mean with the target skewness and target kurtosis set in Step S4. Skewness coefficient and kurtosis coefficient The calculation formula is: ; ; in, These are the values at discrete points in the time history; The mean; This represents the total number of sampling points for this time period.
[0071] This step aims to verify the effectiveness of the nonlinear mapping and autocorrelation function correction strategy based on high-order Hermite polynomials in step S4, ensuring that the simulated ground motion truly reflects the peaks and asymmetric characteristics of the actual record that deviate from the Gaussian distribution, thereby ensuring that the cumulative damage effect of the structure can be accurately induced in the subsequent nonlinear dynamic time history analysis of the structure.
[0072] In step S7, the system integration and implementation of the visualization stochastic simulation platform specifically includes the following sub-steps: Step S71: Platform Architecture and Development Environment Setup. A visual human-computer interaction interface (GUI) is built based on the Matlab App Designer or GUIDE development environment. The platform architecture follows the Model-View-Controller (MVC) design pattern, decoupling the data processing logic (Model), interface display (View), and user interaction response (Controller). The model layer encapsulates all the core algorithms involved in steps S1 to S6, including data filtering, parameter identification, Copula probabilistic modeling, Hermite nonlinear mapping, and POD randomization algorithm. To ensure code execution efficiency and cross-platform compatibility, those skilled in the art can also use programming languages and development environments such as Python (based on PyQt or Tkinter libraries), C++ (based on the Qt framework), or C# (based on the WPF framework) to implement the above functions. These different programming implementations constitute equivalent embodiments of the present invention.
[0073] Step S72: Implementation of the parameter configuration module. Parameter input controls are deployed on the GUI interface to receive simulation conditions set by the user. Specifically, this includes: a drop-down menu for selecting the site category (corresponding to the four site groupings in step S1), and a function for setting the number of simulation samples. The module includes numerical input boxes and text editing boxes for setting target response spectrum characteristic parameters (such as magnitude, fault distance, and other conditional parameters). The parameter configuration module also includes validity verification logic; when the user-input parameters exceed a preset physical reasonable range (e.g., damping ratio less than 0 or greater than 1), an error message mechanism is triggered. This module converts the user-input physical parameters into a data structure recognized internally by the algorithm, which then serves as the input parameters for the subsequent computational kernel.
[0074] Step S73: The computation kernel is encapsulated and the design callback function is called in response to the start simulation command from the interface. This callback function triggers the execution of the background computation kernel. The computation kernel calls the functional sub-functions sequentially according to the data flow order: first, the probability model sub-function is called, which generates the model parameter vector from the joint probability distribution model established in step S3 according to the representative point set selection strategy in step S51; then, the evolution power spectrum generation sub-function is called, which calculates the target spectrum using the formula in step S2; next, the POD spectrum representation sub-function (step S53) is called to generate the Gaussian process; finally, the non-Gaussian mapping sub-function (step S54) is called to output the non-Gaussian ground motion acceleration time history. To improve the computational response speed, the computation kernel adopts matrix-based operations and enables multi-threaded parallel computation mode.
[0075] Step S74: The multi-dimensional visualization module platform integrates a drawing engine to present calculation results graphically in real time. The visualization content covers three dimensions: time domain, frequency domain, and statistical domain. Time-domain visualization: Plotting the generated acceleration time-history curves Velocity time history and displacement time history; Frequency domain visualization: Plotting the generated three-dimensional surface of the evolving power spectral density. And power spectrum slices at specific times, which intuitively reflect the time-frequency non-stationary characteristics of earthquake motion; Statistical domain display: Plot a comparison between the response spectrum (mean and standard deviation envelope) of the simulated sample and the target response spectrum, as well as statistical histograms of kurtosis and skewness of the simulated sample.
[0076] The drawing control supports interactive operations, including zooming, panning, data point snapping, and legend switching.
[0077] Step S75: Data Interaction and Export Interface. Establish a data persistence interface to support the export of simulated seismic motion data and corresponding metadata to a common engineering data format. Export formats include, but are not limited to, text files (.txt, .dat) suitable for finite element analysis software (such as OpenSees, ABAQUS, ANSYS), comma-separated value files (.csv) suitable for Excel processing, and binary data files (.mat) suitable for subsequent processing in Matlab. The exported content includes not only acceleration time series data but also the corresponding sampling frequency, duration, peak acceleration, and generated response spectrum data, thus achieving seamless integration with the structural dynamics analysis stage.
[0078] See attached document Figure 5The present invention provides a complex seismic motion uncertainty quantification stochastic simulation system, which includes: a data filtering module, a parameter identification module, a probability modeling module, a non-Gaussian model construction module, a stochastic simulation module, and a visualization interaction module.
[0079] The data filtering module is used to filter and group measured ground motion records from the Pacific Earthquake Engineering Center database. The module executes mainshock and aftershock sequence identification logic and near-fault pulse identification logic, and completes site grouping based on set shear wave velocity thresholds.
[0080] The parameter identification module is used to identify model parameters of the evolving power spectral density model based on measured ground motion records. For near-fault pulse ground motions, the module is also used to identify pulse function model parameters. The module integrates wavelet packet transform and genetic optimization algorithms to retrieve model parameters from measured acceleration time histories.
[0081] The probabilistic modeling module is used to fit the optimal marginal probability distribution of model parameters and construct the correlation structure between parameters using Copula functions, thus establishing a joint probability distribution model of the model parameters. The probabilistic modeling module includes a probability distribution model library and a Copula function library, and determines the optimal model type using information criteria.
[0082] The non-Gaussian model building module is used to establish higher-order Hermite polynomial models and autocorrelation function mappings. The module receives the set target skewness and kurtosis, solves for the Hermite polynomial coefficients, and establishes the transformation relationship between standard Gaussian and non-Gaussian processes.
[0083] The stochastic simulation module generates random parameter samples from the joint probability distribution model using a representative point set selection strategy, and generates non-Gaussian ground motion acceleration time history samples using a spectral representation based on eigenorthogonal decomposition. The stochastic simulation module generates random seed parameters based on the probability model and synthesizes ground motion time histories using fast Fourier transform techniques.
[0084] The visualization and interaction module provides a user interface and displays simulation results and verification comparison charts. Users can select the type of ground motion, set site conditions, input the number of simulations, and view the generated time-history curves, power spectral density curves, and response spectrum comparison charts through this module.
Claims
1. A method for quantifying the uncertainty of complex seismic motion, characterized in that, The method includes the following steps: S1: Screen measured ground motion records from the earthquake engineering center database, identify non-stationary ground motions, mainshock-aftershock sequence-type ground motions and near-fault pulse ground motions in the measured ground motion records, and group the measured ground motion records according to the average shear wave velocity; S2: The evolutionary power spectral density model is adopted as a unified theoretical framework, and the model parameters of the evolutionary power spectral density model are identified based on the measured ground motion records. S3: Perform statistical analysis on the identified model parameters, fit the optimal marginal probability distribution of the model parameters, and construct the correlation structure between the model parameters based on Copula function theory, thereby establishing a joint probability distribution model of the model parameters; S4: Establish a high-order Hermite polynomial model for implementing non-Gaussian transformation, and establish a mapping relationship between the autocorrelation function of the standard Gaussian process and the autocorrelation function of the non-Gaussian process through the high-order Hermite polynomial model. S5: Generate random parameter samples from the joint probability distribution model using a representative point set selection strategy, substitute the random parameter samples into the evolution power spectral density model to determine the target evolution power spectral density, generate a Gaussian random process using a spectral representation method based on intrinsic orthogonal decomposition, and convert the Gaussian random process into a non-Gaussian ground motion acceleration time history sample using the mapping relationship. S6: Calculate the response spectrum of the non-Gaussian ground motion acceleration time history sample and compare it with the response spectrum of the measured ground motion record to verify the engineering applicability of the random simulation sample; S7: Develop a visual seismic motion stochastic simulation platform based on Matlab GUI, and integrate the process into the visual seismic motion stochastic simulation platform.
2. The method for quantifying the uncertainty of complex seismic motion according to claim 1, characterized in that, In step S2, the evolved power spectral density model is constructed as follows: The stationary power spectral density function is represented by the Clough-Penzien spectrum, which includes parameters such as site soil damping ratio, bedrock damping ratio, site soil dominant circular frequency, and bedrock dominant circular frequency. The time-frequency modulation function includes time parameters that control the attenuation of ground motion energy and the arrival time of the peak.
3. The method for quantifying the uncertainty of complex seismic motion according to claim 1, characterized in that, In step S2, the specific process of identifying the model parameters of the evolved power spectral density model based on the measured ground motion records includes: The acceleration time history of the measured ground motion record is decomposed into wavelet packet coefficients of different times and frequencies using wavelet packet transform; The wavelet packet coefficients are reconstructed into sub-signals of different frequency bands using inverse wavelet packet transform; The energy of the sub-signal is calculated according to the definition of the evolved power spectral density, thereby obtaining the energy distribution function corresponding to the evolved power spectral density estimated based on wavelet packet transform; Based on the principle of equal energy in the frequency domain, a genetic algorithm is used to fit the measured ground motion record with the energy distribution function as the target, thereby identifying the model parameters in the evolved power spectral density model.
4. The method for quantifying the uncertainty of complex seismic motion according to claim 1, characterized in that, In step S3, the specific process of fitting the optimal marginal probability distribution of the model parameters includes: The optimal marginal probability distribution type of the model parameters is selected from the candidate distribution models using the KS test method and the Akaike information content criterion.
5. The method for quantifying the uncertainty of complex seismic motion according to claim 1, characterized in that, In step S3, the specific process of constructing the correlation structure between the model parameters based on Copula function theory includes: The joint probability distribution of the model parameters is described using a Copula function, and the optimal Copula function type is determined from the candidate Copula function types using the Akaike information content criterion. The correlation structure is then constructed based on the optimal Copula function type. When the measured ground motion record is a mainshock-aftershock sequence type ground motion, a high-dimensional correlation between the parameters of the mainshock evolution power spectral density model is established using a Copula structure, and the cross-correlation between the mainshock model parameters and the aftershock model parameters is established using a binary Copula function.
6. The method for quantifying the uncertainty of complex seismic motion according to claim 1, characterized in that, In step S4, the specific process of establishing a high-order Hermite polynomial model for implementing the non-Gaussian transform includes: Calculate the statistical characteristics of skewness and kurtosis of the measured ground motion acceleration time history; Based on the aforementioned skewness and kurtosis statistical characteristics, the Hermite shape factor is calculated using linear moments; Construct a second-order linear system containing the Hermite shape coefficients and polynomial coefficients; By solving the second-order linear system, the transformation relationship between the standard Gaussian process and the target non-Gaussian process is determined, and the analytical mapping model between the autocorrelation function of the standard Gaussian process and the autocorrelation function of the non-Gaussian process is derived.
7. The method for quantifying the uncertainty of complex seismic motion according to claim 1, characterized in that, In step S5, the spectral representation based on eigenorthogonal decomposition generates a Gaussian random process by embedding a fast Fourier transform. The specific steps for generating the Gaussian random process include: The cross-spectral density matrix is constructed using the target evolutionary power spectral density generated by the evolutionary power spectral density model. The cross-spectral density matrix is subjected to eigenorthogonal decomposition to obtain eigenvalues and eigenvectors; The Gaussian random process is synthesized by combining the eigenvalues, the eigenvectors, and the random phase angle using the spectral representation.
8. The method for quantifying the uncertainty of complex seismic motion according to claim 1, characterized in that, When the measured ground motion record is a mainshock-aftershock sequence type ground motion, the evolution power spectral density model includes the mainshock evolution power spectral density model and the aftershock evolution power spectral density model. The mainshock evolution power spectral density model and the aftershock evolution power spectral density model correspond to the independently set mainshock duration and aftershock duration, respectively, and a time interval is set between the mainshock and the aftershock. The model parameters include the dominant circular frequency of the site soil for both the mainshock and aftershock, the site soil damping ratio, the spectral intensity factor, and the time-frequency modulation parameters.
9. The method for quantifying the uncertainty of complex seismic motion according to claim 1, characterized in that, When the measured ground motion record is a near-fault pulse ground motion, step S2 further includes the step of decomposing the near-fault pulse ground motion into low-frequency pulse components and high-frequency residual components for separate simulation; The steps of the separate simulations include: extracting the velocity time history of the strongest pulse direction using a pulse recognition method based on continuous wavelet transform; The velocity-time history of the strongest pulse direction is fitted using the Gabor pulse function model to identify pulse model parameters including peak pulse velocity, pulse occurrence time, pulse period, pulse half-wave number, and pulse phase angle, which are then used to generate the low-frequency pulse component. The residual acceleration time history is obtained by subtracting the fitted pulse component from the original record, and the parameters of the residual acceleration time history are identified using the evolved power spectral density model to generate the high-frequency residual component. The generated low-frequency pulse component and the high-frequency residual component are superimposed to synthesize near-fault pulse ground motion.
10. A stochastic simulation system for quantifying the uncertainty of complex ground motions, applied to the method for quantifying the uncertainty of complex ground motions as described in claims 1-9, characterized in that, include: The data filtering module is used to filter and group measured ground motion records from the Pacific Earthquake Engineering Center database; The parameter identification module is used to identify the model parameters of the evolution power spectral density model based on the measured ground motion records, and for near-fault pulse ground motion, it is also used to identify the pulse function model parameters. The probability modeling module is used to fit the optimal marginal probability distribution of the model parameters and use the Copula function to construct the correlation structure between the parameters, thereby establishing a joint probability distribution model of the model parameters. The non-Gaussian model building module is used to establish higher-order Hermite multinomial models and autocorrelation function mapping relationships; The random simulation module is used to generate random parameter samples from the joint probability distribution model using a representative point set selection strategy, and to generate non-Gaussian ground motion acceleration time history samples in conjunction with the spectral representation based on intrinsic orthogonal decomposition. The visualization and interaction module is used to provide a user interface and display simulation results and verification comparison charts.