A method and apparatus for generating seismic waves based on target spectrum and Green's function

Abstract

This invention relates to the field of numerical simulation and computational technology for seismic wave propagation, and discloses a method and apparatus for generating seismic ground motions based on target spectra and Green's functions. It aims to address the difficulty in achieving a balance between physical consistency, computational efficiency, and spectral compatibility in existing methods. The scheme mainly includes: preprocessing environmental noise records and extracting empirical Green's functions through unsupervised clustering; constructing a spatially continuous propagation tensor by performing depth correction and gradient interpolation based on surface wave eigenfunctions; discretizing the kinematically finite fault model into sub-source units and generating long-period seismic ground motions including path effects using propagation tensor convolution; optimizing source parameters through Bayesian inversion using the GMPE target spectrum as a constraint; and finally generating high-frequency components and fusing them with the long-period waveform to obtain a broadband seismic ground motion time history. This invention achieves a unification of physical propagation mechanisms and statistical spectral constraints, significantly improving computational efficiency while ensuring spectral compatibility, and is particularly suitable for basin areas.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation and calculation technology for seismic wave propagation, and specifically to a method and apparatus for generating seismic ground motion based on target spectrum and Green's function. Background Technology

[0002] With the acceleration of urbanization, the seismic safety of major infrastructure, urban clusters, and energy hubs is becoming increasingly prominent. Accurate and reliable seismic motion input is the foundation for seismic design of engineering structures, seismic hazard analysis, and urban resilience assessment. However, in large sedimentary basin areas (such as the Sichuan Basin and the Guandong Basin), the basin effect significantly amplifies seismic waves, and the limited observational data and computational resources make it a challenging problem to quickly generate seismic motion time histories that both meet the target response spectrum requirements and accurately reflect the physical processes of wave propagation.

[0003] Currently, seismic motion simulation methods are mainly divided into two categories:

[0004] The first category is deterministic numerical simulation methods, such as the finite difference method, the finite element method, or the spectral element method. These methods, by establishing a detailed three-dimensional subsurface velocity structure model, can relatively completely characterize the propagation process of seismic waves in complex media, including basin effects and topographic effects. However, this method has extremely high computational resource requirements; a single simulation typically consumes [amount missing]. The core hourly rate is on the order of magnitude and highly dependent on high-resolution velocity models, making it difficult to promote and apply in areas lacking detailed geological data.

[0005] The second category is empirical or stochastic simulation methods, such as EXSIM and SMSIM. These methods are based on source spectrum models and empirical attenuation relationships, offering high computational efficiency and ease of operation, and are widely used in engineering practice. However, to simplify calculations, these methods typically involve significant simplification of the wave propagation process, making it difficult to accurately reproduce the phase characteristics and basin amplification effects of seismic waves, especially in the long-period range, where simulation results often deviate considerably from measured records.

[0006] Furthermore, traditional methods often employ empirical spectral scaling or amplitude adjustment techniques when dealing with the matching problem between ground motion and target response spectra. While these approaches can ensure that the response spectrum of the synthesized ground motion meets the specifications, they often disrupt the original phase relationship and wavefield consistency of the waveform, resulting in unclear physical meaning of the generated ground motion time history and making it difficult to use for refined structural nonlinear analysis.

[0007] The above problems are particularly prominent in areas lacking strong earthquake observation data. How to generate broadband ground motion time histories that both meet the target spectrum constraints and have physical consistency without relying on large-scale three-dimensional velocity models and a large amount of measured records has become a key problem that current ground motion simulation technology urgently needs to overcome. Summary of the Invention

[0008] This invention aims to address the problem that existing seismic motion simulation methods struggle to balance physical consistency, computational efficiency, and spectral compatibility, and are unable to rapidly generate broadband seismic motion time histories that both conform to the target response spectrum and accurately reflect the basin wave propagation characteristics under limited observation data. The invention proposes a seismic motion generation method and apparatus based on the target spectrum and Green's function.

[0009] The technical solution adopted by the present invention to solve the above-mentioned technical problems is as follows:

[0010] In a first aspect, the present invention provides a method for generating seismic motion based on a target spectrum and a Green's function, the method comprising:

[0011] Step 1: Preprocess the continuous environmental noise records of the target area, and evaluate the waveform similarity of the preprocessed noise records using an unsupervised clustering algorithm to extract the empirical Green's function. The empirical Green's function is used to characterize the wave propagation characteristics and basin amplification effect within the target area.

[0012] Step 2: Based on the surface wave eigenfunctions, perform depth correction on the extracted empirical Green's function to convert it from the surface response to the Green's function at the source depth; use the gradient approximation method to perform spatial interpolation on the depth-corrected Green's function to construct a spatially continuous propagation tensor covering the fault plane within the target area;

[0013] Step 3: Based on the kinematic finite fault model, the fault plane in the target area is discretized into multiple sub-source units. Using the spatial continuous propagation tensor, the source time function of each sub-source unit is convolved and superimposed with the corresponding Green's function to generate a long-period ground motion waveform that includes path and site effects.

[0014] Step 4: Using the target response spectrum predicted by the empirical ground motion prediction equation as a constraint, establish a posterior distribution model of the source parameters; with the goal of minimizing the residual between the synthetic spectrum of the long-period ground motion waveform and the target response spectrum, perform global optimization inversion of the source parameters using the sequential Monte Carlo annealing sampling method to obtain the optimal source parameters.

[0015] Step 5: Based on the optimal source parameters, a short-period high-frequency ground motion component matching the long-period ground motion waveform is generated using a high-frequency radiation model; the long-period ground motion waveform and the short-period high-frequency ground motion component are fused in the frequency domain using a smoothing weighting function to generate the final ground motion time history covering the target wideband range.

[0016] Furthermore, step 1 specifically includes:

[0017] Step 11: Perform segmentation, bandpass filtering, spectral whitening, and amplitude normalization on the continuous environmental noise record to obtain the preprocessed noise record. The formula for amplitude normalization is:

[0018] or ;

[0019] in, This indicates the original noise waveform at time [time]. amplitude, This represents the normalized amplitude. Indicates taking The absolute value, express The root mean square value;

[0020] Step 12: For the preprocessed noise records, calculate the cross-correlation function between each pair of stations to obtain the cross-correlation function set. The formula for calculating the cross-correlation function is:

[0021] ;

[0022] in, Indicates the first The station and the first Time shift between individual stations The cross-correlation function value at the location, Indicates the first Each station is at time Preprocessed noise recordings Indicates the first Each station is at time Preprocessed noise recording;

[0023] Step 13: For each pair of cross-correlation functions in the set of cross-correlation functions, calculate the dynamic time warped distance. The calculation formula is as follows:

[0024] ;

[0025] in, Indicates the first Article and No. The dynamic time-warped distance between the cross-correlation functions. This indicates taking the minimum value. This represents the optimal curved path. This represents the index pair of corresponding points in two cross-correlation functions. Represents corresponding points in two noise waveforms and Local distance metric Indicates the first Cross-correlation function at point amplitude, Indicates the first Cross-correlation function at point The amplitude;

[0026] Step 14: Construct the kernel function matrix based on the dynamic time warping distance, the expression of which is:

[0027] ;

[0028] in, This represents the kernel function value after mapping by kernel principal component analysis, used to measure the... Article and No. The similarity between cross-correlation functions, This represents the kernel function bandwidth, used to control the similarity decay scale. Represents the natural exponential function;

[0029] Step 15: After performing kernel principal component analysis to reduce the dimensionality of the kernel function matrix, clustering is performed using a Gaussian mixture model. The cluster number corresponding to the minimum Bayesian information criterion value is selected, and the cluster with the smallest variance is chosen. The center of this cluster is used as the empirical Green's function. The expression for the Bayesian information criterion is:

[0030] ;

[0031] in, Represents the Bayesian information criterion value. This represents the number of clusters in the Gaussian mixture model. This indicates the number of cross-correlation functions participating in the clustering. This represents the maximum likelihood function value of the Gaussian mixture model. It represents the natural logarithm.

[0032] Furthermore, in step 2, depth correction is performed on the extracted empirical Green's function based on the surface wave eigenfunctions, specifically including:

[0033] Based on the eigenfunctions of Love and Rayleigh waves, a depth correction factor is calculated to correct the empirical Green's function from the surface response to the focal depth. At this point, the depth-corrected Green's function is obtained, and its approximate relationship is:

[0034] ;

[0035] ;

[0036] in, and Representing depth Green's functions for the Love wave and Rayleigh wave at that location. and These represent the Green's functions at the Earth's surface corresponding to the Love wave and the Rayleigh wave, respectively. and These represent the eigenfunction values ​​of the Love wave and Rayleigh wave at the corresponding depths, respectively.

[0037] Spatial interpolation of the depth-corrected Green's function is performed using the gradient approximation method, specifically including:

[0038] For locations within the fault plane that are not directly observed, a first-order gradient approximation is used to spatially interpolate the depth-corrected Green's function, constructing a spatially continuous propagation tensor, the expression of which is:

[0039] ;

[0040] in, Indicates the interpolation point to be determined. Green's function value at that point, Represents a known point The depth-corrected Green's function at that point, Represents the gradient operator. Represent the Green's function at a given point The gradient vector at that point, This represents the displacement vector from the known point to the interpolation point to be determined;

[0041] The interpolation scaling factor is derived from the first-order gradient approximation formula and used to evaluate the accuracy of the interpolation result. Its formula is as follows:

[0042] ;

[0043] in, This represents the interpolation scaling factor. The closer the value is to 1, the higher the interpolation accuracy. and This represents the spatial position vector of two known points. Indicates starting from a known point to a known point Green's function value, and Indicates the component index of the Green's function. This represents its gradient vector.

[0044] Furthermore, step 3 specifically includes:

[0045] Step 31: Discretize the kinematic finite tomographic model into multiple sub-source units, each sub-source unit corresponding to a spatial location;

[0046] Step 32: In the frequency domain, using the spatially continuous propagation tensor, convolve the source moment tensor of each sub-source element with the Green's function corresponding to the location of that sub-source element and sum them to obtain the seismic motion spectrum at the observation point:

[0047] ;

[0048] in, Indicates the first observation point The seismic motion spectrum of the components, Represents angular frequency. Indicates the total number of sub-source units. Indicates the first Individual source units The spectrum of the Green's function to the observation point Indicates the first The source moment tensor spectrum of each sub-source element;

[0049] Step 33: Perform an inverse Fourier transform on the seismic motion spectrum to obtain the long-period seismic motion waveform. The transformation formula is as follows:

[0050] ;

[0051] in, This represents the long-period ground motion waveform in the time domain. Indicates time, This represents the inverse Fourier transform operator.

[0052] Furthermore, in step 4, the posterior distribution model is:

[0053] ;

[0054] in, The posterior probability density function representing the source parameters. This represents the source parameters to be inverted. Empirical Seismic Ground Motion Prediction Equation The predicted target response spectrum, Indicates direct proportion. This represents the natural exponential function. Indicates based on current source parameters The spectral function of the synthesized long-period seismic motion waveform. Indicates matrix transpose. Represents the covariance matrix of the observed data The inverse matrix, Represents the regularization parameter. Represents the source slip vector. Represents the covariance matrix of the sliding distribution The inverse matrix.

[0055] Furthermore, the source slip vector The prior distribution is constrained by the Von Kármán autocorrelation model, and its spatial covariance function is:

[0056] ;

[0057] in, Indicates distance The covariance of the slip between two points Indicates the spatial distance between two points. Represents the sliding variance. Indicates the sliding standard deviation. This represents the power exponent, used to control the roughness of the sliding distribution. Represents the gamma function operator. This indicates the relevant scale, used to control the attenuation distance of spatial correlation. Indicates the order is The modified Bessel function.

[0058] Furthermore, in step 4, the source parameters are globally optimized and inverted using the sequential Monte Carlo annealing sampling method, specifically including:

[0059] The posterior distribution samples of the source parameters were obtained using the sequential Monte Carlo annealing sampling method. An intermediate posterior distribution was constructed by introducing a temperature coefficient.

[0060] ;

[0061] in, Indicates the first The intermediate posterior distribution at each temperature step Indicates the first Temperature coefficient at each temperature step , This represents the prior distribution of the source parameters. Represents the likelihood function;

[0062] With temperature coefficient The sampling process increments from 0 to 1, starting from the prior distribution. Smooth transition to posterior distribution The final converged sample distribution corresponds to the region that minimizes the residual between the synthetic spectrum and the target response spectrum of the long-period ground motion waveform. The maximum a posteriori solution or the sample mean is selected from the converged sample set to obtain the optimal source parameters. .

[0063] Furthermore, in step 5, the high-frequency radiation model is: The spectral model, whose spectral expression is:

[0064] ;

[0065] in, This represents the high-frequency spectrum corresponding to short-period high-frequency ground motion components. Indicates frequency, This represents the scalar seismic moment determined based on the optimal source parameters. Represents angular frequency. Indicates the angular frequency;

[0066] The turning frequency The calculation formula is:

[0067] ;

[0068] in, Represents a constant. Used to adjust dimensions. Indicates the shear wave velocity of the medium. This indicates stress drop.

[0069] Further, in step 5, the long-period ground motion waveform and the short-period high-frequency ground motion component are fused in the frequency domain using a smoothing weighting function, specifically including:

[0070] Step 51: Within the preset transition frequency band, define a low-frequency weighting function and a high-frequency weighting function such that their sum is 1. The expression is as follows:

[0071] ;

[0072] ;

[0073] in, This represents the low-frequency weighting function. Represents a high-frequency weighting function. Indicates frequency, Represents the cosine function. Pi is a constant. Indicates the center frequency of the transition band. This represents half of the transition bandwidth;

[0074] Step 52: Using the low-frequency weighting function and the high-frequency weighting function, the low-frequency spectrum corresponding to the long-period ground motion waveform and the high-frequency spectrum corresponding to the short-period high-frequency ground motion component are weighted and superimposed to obtain a broadband spectrum. The superposition formula is as follows:

[0075] ;

[0076] in, Indicates broadband spectrum. This represents the low-frequency spectrum corresponding to long-period ground motion waveforms. This represents the high-frequency spectrum corresponding to the short-period high-frequency ground motion component.

[0077] Step 53: Perform an inverse Fourier transform on the broadband spectrum combined with the phase spectrum to obtain the final ground motion time history covering the target broadband range:

[0078] ;

[0079] in, Indicates the final earthquake time history, Indicates time, This represents the inverse Fourier transform operator. Represents the phase term in complex exponential form. Represents the natural constant. Represents the imaginary unit. Represents the phase spectrum, where frequency is... The function.

[0080] Secondly, the present invention provides a seismic motion generation device based on a target spectrum and a Green's function, for implementing the seismic motion generation method based on a target spectrum and a Green's function as described in the first aspect, the device comprising:

[0081] The Green's function extraction module is used to preprocess the continuous environmental noise records of the target area, and to evaluate the waveform similarity of the preprocessed noise records through an unsupervised clustering algorithm to extract the empirical Green's function. The empirical Green's function is used to characterize the wave propagation characteristics and basin amplification effect within the target area.

[0082] The Green's function correction and interpolation module is used to perform depth correction on the extracted empirical Green's function based on the surface wave eigenfunction, converting it from the surface response to the Green's function at the source depth; the gradient approximation method is used to perform spatial interpolation on the depth-corrected Green's function to construct a spatially continuous propagation tensor covering the fault plane in the target area.

[0083] The long-period ground motion generation module is used to discretize the fault plane in the target area into multiple sub-source units based on the kinematic finite fault model. Using the spatial continuous propagation tensor, the source time function of each sub-source unit is convolved and superimposed with the corresponding Green's function to generate a long-period ground motion waveform that includes path and site effects.

[0084] The source parameter inversion module is used to establish a posterior distribution model of the source parameters based on the target response spectrum predicted by the empirical ground motion prediction equation as a constraint; with the goal of minimizing the residual between the synthetic spectrum of the long-period ground motion waveform and the target response spectrum, the source parameters are globally optimized and inverted using the sequential Monte Carlo annealing sampling method to obtain the optimal source parameters.

[0085] The broadband fusion module is used to generate short-period high-frequency ground motion components that match the long-period ground motion waveform based on the optimal source parameters and using a high-frequency radiation model; and to fuse the long-period ground motion waveform and the short-period high-frequency ground motion components in the frequency domain through a smoothing weighting function to generate the final ground motion time history covering the target broadband range.

[0086] The beneficial effects of this invention are as follows: The ground motion generation method and apparatus based on the target spectrum and Green's function provided by this invention, by integrating the environmental noise Green's function with Bayesian inversion optimization, can achieve the unification of physical propagation mechanism and statistical spectrum constraints without relying on a large-scale three-dimensional velocity model. The generated ground motion time history not only realistically reflects the basin wave propagation characteristics but also accurately matches the target response spectrum, avoiding the problem of destroying waveform phase consistency in order to meet spectral requirements in traditional methods. At the same time, this invention improves computational efficiency compared with the traditional three-dimensional finite difference method, and can quickly realize regionally transferable ground motion synthesis under limited observation data conditions. It fills the technical gap in the existing technology where it is difficult to balance physical consistency, computational efficiency and spectral compatibility, and provides efficient and reliable technical support for seismic design of major projects and seismic hazard assessment of urban agglomerations. Attached Figure Description

[0087] Figure 1 A schematic flowchart of the seismic motion generation method based on the target spectrum and Green's function provided for the embodiment;

[0088] Figure 2 The diagram shows the source slip inversion results under Laplacian smoothing constraints provided in the embodiment; where (a) is the average slip distribution under Laplacian smoothing constraints, (b) is the maximum a posteriori probability slip distribution under Laplacian smoothing constraints, and (c) is the slip standard deviation distribution under Laplacian smoothing constraints.

[0089] Figure 3 The following is a schematic diagram of the source slip inversion results under the Von Kármán prior constraint provided in the example; where (a) is the average slip distribution under the Von Kármán prior constraint, (b) is the maximum a posteriori probability slip distribution under the Von Kármán prior constraint, and (c) is the slip standard deviation distribution under the Von Kármán prior constraint.

[0090] Figure 4The following is a schematic diagram of the spectral function error before Bayesian update provided in the embodiment; wherein, (a) is the distribution result of spectral function error with source distance under the condition of 1s before Bayesian update, (b) is the distribution result of spectral function error with source distance under the condition of 3s before Bayesian update, (c) is the distribution result of spectral function error with source distance under the condition of 5s before Bayesian update, and (d) is the distribution result of spectral function error with source distance under the condition of 10s before Bayesian update;

[0091] Figure 5 The following is a schematic diagram of the spectral function error after Bayesian update provided in the embodiment; wherein, (a) is the distribution result of spectral function error with source distance under the condition of Bayesian update period of 1s, (b) is the distribution result of spectral function error with source distance under the condition of Bayesian update period of 3s, (c) is the distribution result of spectral function error with source distance under the condition of Bayesian update period of 5s, and (d) is the distribution result of spectral function error with source distance under the condition of Bayesian update period of 10s;

[0092] Figure 6 The following is a schematic diagram of the spatial distribution of the empirical Green's function extracted from the basin region and the error of the response spectrum function under different periods, provided for the example; wherein, (a) is the empirical Green's function extracted from MESO stations in the basin region, (b) is the spatial distribution of the error of the response spectrum function in the basin region under the condition of period 1s, (c) is the spatial distribution of the error of the response spectrum function in the basin region under the condition of period 3s, and (d) is the spatial distribution of the error of the response spectrum function in the basin region under the condition of period 5s.

[0093] Figure 7 The diagram illustrates a comparison of the ground motion parameters of the method in this embodiment with the EXSIM method and the finite difference method, provided for the implementation example; wherein, (a) shows the comparison results of peak ground acceleration under different source distances, (b) shows the comparison results of peak ground velocity under different source distances, and (c) shows the comparison results of pseudospectral acceleration under different source distances.

[0094] Figure 8 A schematic diagram comparing the computational costs of the method in this embodiment with those of the EXSIM method and the finite difference method, provided for the purpose of this embodiment.

[0095] Figure 9 A schematic diagram of the seismic motion generation device based on the target spectrum and Green's function provided in the embodiment. Detailed Implementation

[0096] Existing seismic motion simulation methods rely excessively on high-precision three-dimensional velocity models, oversimplify wave propagation physics, and crudely handle target spectrum matching requirements. This makes it difficult to achieve a balance between physical consistency, computational efficiency, and spectral compatibility, hindering the rapid generation of broadband seismic time histories that both conform to the target response spectrum and accurately reflect basin wave propagation characteristics under limited observational data. Therefore, this invention proposes a technical solution.

[0097] In this invention, firstly, unsupervised clustering analysis is performed on environmental noise records of the target area to extract an empirical Green's function that truly reflects the regional wave propagation characteristics and basin amplification effect. Then, the empirical Green's function is deeply corrected based on surface wave eigenfunctions, and a spatially continuous propagation tensor covering the entire fault plane is constructed using the gradient approximation method, thereby achieving a physical description of the propagation path without relying on a fine velocity model. Next, the fault plane is discretized into multiple sub-source units, and the spatially continuous propagation tensor is convolved and superimposed to generate a long-period ground motion waveform that includes path and site effects. Furthermore, the target response spectrum predicted by the empirical ground motion prediction equation is used as a constraint, and the source parameters are globally optimized using Bayesian inversion and sequential Monte Carlo annealing sampling methods, so that the synthesized spectrum and the target spectrum are statistically consistent, thus avoiding the physical distortion caused by the forced waveform adjustment to meet spectral requirements in traditional methods. Finally, a high-frequency radiation model is used to generate short-period components, and the short-period components are fused with the long-period waveform through a frequency domain smoothing weight function to obtain the final ground motion time history covering a wide frequency band. This invention, through the implementation principle of combining physical drive and data constraints, fundamentally solves the problem of the difficulty in balancing physical consistency, computational efficiency and spectral compatibility in the prior art.

[0098] The technical solutions in this embodiment 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.

[0099] Figure 1 A flowchart illustrating a seismic motion generation method based on the target spectrum and Green's function is shown. Please refer to [link / reference]. Figure 1 The method includes the following steps:

[0100] Step 1: Preprocess the continuous environmental noise records of the target area, and evaluate the waveform similarity of the preprocessed noise records using an unsupervised clustering algorithm to extract the empirical Green's function. The empirical Green's function is used to characterize the wave propagation characteristics and basin amplification effect within the target area.

[0101] Specifically, this step aims to extract a stable and repeatable empirical Green's function from continuous environmental noise records in the target area using an unsupervised learning method. This function is used to characterize the propagation characteristics of seismic waves in the target area and the basin amplification effect.

[0102] In this embodiment, step 1 specifically includes steps 11 to 15:

[0103] Step 11: Perform segmentation, bandpass filtering, spectral whitening, and amplitude normalization on the continuous environmental noise record to obtain the preprocessed noise record. The formula for amplitude normalization is:

[0104] or ;

[0105] in, This indicates the original noise waveform at time [time]. amplitude, This represents the normalized amplitude. Indicates taking The absolute value, express The root mean square value.

[0106] In practical applications, the continuous environmental noise recorded by each station in the target area is first preprocessed. In this embodiment, the continuous environmental noise waveform is segmented into 30-minute time windows to ensure a moderate time series length, controllable computational load, and reduced interference from non-stationary noise sources. Then, bandpass filtering is applied to the noise waveform of each time window, with a frequency range set to 1-20 seconds, corresponding to the main surface wave period range of the basin. The filtered noise records are then subjected to spectral whitening and amplitude normalization to reduce the influence of near-field source energy differences and occasional high-energy events. Through the above preprocessing, standardized noise records suitable for subsequent analysis are obtained.

[0107] Step 12: For the preprocessed noise records, calculate the cross-correlation function between each pair of stations to obtain the cross-correlation function set. The formula for calculating the cross-correlation function is:

[0108] ;

[0109] in, Indicates the first The station and the first Time shift between individual stations The cross-correlation function value at the location, Indicates the first Each station is at time Preprocessed noise recordings Indicates the first Each station is at time The preprocessed noise record.

[0110] In practical applications, the cross-correlation function between each pair of stations is calculated for the preprocessed noise records to extract wave propagation information between stations. By segmenting and averaging continuous records of more than one year, a stable set of cross-correlation functions can be obtained.

[0111] Step 13: For each pair of cross-correlation functions in the set of cross-correlation functions, calculate the dynamic time warped distance. The calculation formula is as follows:

[0112] ;

[0113] in, Indicates the first Article and No. The dynamic time-warped distance between the cross-correlation functions. This indicates taking the minimum value. This represents the optimal curved path. This represents the index pair of corresponding points in two cross-correlation functions. Represents corresponding points in two noise waveforms and Local distance metric Indicates the first Cross-correlation function at point amplitude, Indicates the first Cross-correlation function at point The amplitude.

[0114] Specifically, to select the most stable and reliable waveforms from the set of cross-correlation functions, this step introduces a dynamic time warping algorithm to evaluate the similarity between each pair of cross-correlation functions. Dynamic time warping can effectively handle the slight scaling of waveforms on the time axis and more accurately measure the similarity between waveforms. The smaller the dynamic time warping distance, the more similar the waveform shapes of the two cross-correlation functions are.

[0115] Step 14: Construct the kernel function matrix based on the dynamic time warping distance, the expression of which is:

[0116] ;

[0117] in, This represents the kernel function value after mapping by kernel principal component analysis, used to measure the... Article and No. The similarity between cross-correlation functions, This represents the kernel function bandwidth, used to control the similarity decay scale. This represents the natural exponential function.

[0118] Specifically, the kernel function matrix is ​​used to map waveform similarity to a high-dimensional feature space. During the calculation of the kernel function value, the similarity between waveforms is transformed into numerical values ​​in the kernel matrix through kernel function transformation, laying the foundation for subsequent dimensionality reduction and clustering.

[0119] Step 15: After performing kernel principal component analysis to reduce the dimensionality of the kernel function matrix, clustering is performed using a Gaussian mixture model. The cluster number corresponding to the minimum Bayesian information criterion value is selected, and the cluster with the smallest variance is chosen. The center of this cluster is used as the empirical Green's function. The expression for the Bayesian information criterion is:

[0120] ;

[0121] in, Represents the Bayesian information criterion value. This represents the number of clusters in the Gaussian mixture model. This indicates the number of cross-correlation functions participating in the clustering. This represents the maximum likelihood function value of the Gaussian mixture model. It represents the natural logarithm.

[0122] In practical applications, kernel principal component analysis (KPCA) is used to reduce the dimensionality of the kernel function matrix, extracting the main feature components and reducing the data dimensionality while preserving the main differences between waveforms. After dimensionality reduction, a Gaussian mixture model is used to perform cluster analysis on the dimensionality-reduced feature vectors. The number of clusters is optimized using the Bayesian information criterion, selecting the number of clusters that minimizes the Bayesian information criterion value as the optimal number of clusters. The cluster with the smallest variance is then selected from the corresponding clusters, and the center of this cluster is used as the final extracted empirical Green's function.

[0123] Through the above steps, an empirical Green's function with clear physical meaning and high stability is extracted from the environmental noise record, providing a reliable input for subsequent depth correction, spatial interpolation, and ground motion synthesis.

[0124] Step 2: Based on the surface wave eigenfunctions, perform depth correction on the extracted empirical Green's function to convert it from the surface response to the Green's function at the source depth; use the gradient approximation method to perform spatial interpolation on the depth-corrected Green's function to construct a spatially continuous propagation tensor covering the fault plane within the target area.

[0125] This step aims to correct the empirical Green's function extracted in step 1 to the actual source depth, and to construct a spatially continuous propagation tensor covering the entire fault plane through spatial interpolation, so as to provide an accurate propagation operator for subsequent long-period ground motion synthesis.

[0126] In this embodiment, depth correction is performed on the extracted empirical Green's function based on the surface wave eigenfunction, specifically including:

[0127] Based on the eigenfunctions of Love and Rayleigh waves, a depth correction factor is calculated to correct the empirical Green's function from the surface response to the focal depth. At this point, the depth-corrected Green's function is obtained, and its approximate relationship is:

[0128] ;

[0129] ;

[0130] in, and Representing depth Green's functions for the Love wave and Rayleigh wave at that location. and These represent the Green's functions at the Earth's surface corresponding to the Love wave and the Rayleigh wave, respectively. and These represent the eigenfunction values ​​of the Love wave and Rayleigh wave at the corresponding depths, respectively.

[0131] Specifically, the empirical Green's function extracted in step 1 is based on surface station records and reflects the wave propagation characteristics between two points on the Earth's surface. However, the hypocenter of an actual earthquake is located at a certain depth underground. Therefore, the surface Green's function needs to be corrected to the hypocenter depth to accurately describe the wave propagation process from the underground hypocenter to the surface receiving point.

[0132] This embodiment focuses on depth correction for the main components of surface waves: Love waves and Rayleigh waves. Based on the eigenfunctions of Love waves and Rayleigh waves, depth correction factors are calculated, and the empirical Green's function is corrected from the surface response to the source depth. At this depth, the Green's function after depth correction is obtained. Here, the eigenfunctions reflect the relative magnitudes of surface wave amplitudes at different depths. The above formula approximates the Green's function at the source depth by scaling the Green's function observed at the surface according to the ratio of the eigenfunctions.

[0133] In this embodiment, the gradient approximation method is used to perform spatial interpolation on the depth-corrected Green's function, specifically including:

[0134] For locations within the fault plane that are not directly observed, a first-order gradient approximation is used to spatially interpolate the depth-corrected Green's function, constructing a spatially continuous propagation tensor, the expression of which is:

[0135] ;

[0136] in, Indicates the interpolation point to be determined. Green's function value at that point, Represents a known point The depth-corrected Green's function at that point, Represents the gradient operator. Represent the Green's function at a given point The gradient vector at that point, This represents the displacement vector from the known point to the interpolation point to be determined.

[0137] In practical applications, fault planes are continuous planar structures, while the empirical Green's function extracted from station pairs only corresponds to a limited number of station locations. To obtain the Green's function value at any location on the fault plane, spatial interpolation of the depth-corrected Green's function is required.

[0138] This embodiment employs the first-order gradient approximation method for spatial interpolation. For locations within the fault plane that are not directly observed, the Green's function values ​​at nearby unknown points are estimated using the known points' Green's function values ​​and their spatial rate of change. The physical meaning of the above formula is that within a small range near a known point, the change in the Green's function can be approximated as linear; therefore, the function value at an unknown point can be estimated by adding the change along the displacement direction to the known point's function value. This method allows for the construction of a spatially continuous propagation tensor covering the entire fault plane without the need for dense observations, providing Green's functions corresponding to arbitrary sub-source units for subsequent finite fault synthesis.

[0139] In this embodiment, to evaluate the reliability of the interpolation results, this step also introduces an interpolation scaling factor to quantify the interpolation accuracy. The estimation formula for the interpolation scaling factor can be derived from the first-order gradient approximation formula:

[0140] ;

[0141] in, This represents the interpolation scaling factor. The closer the value is to 1, the higher the interpolation accuracy. and This represents the spatial position vector of two known points. Indicates starting from a known point to a known point Green's function value, and Indicates the component index of the Green's function. This represents its gradient vector.

[0142] Interpolation scaling factor The physical meaning is the ratio of the interpolation result to the known point value. When When the value is close to 1, it indicates that the interpolation increment is relatively small compared to the original value, and the interpolation result has high reliability; when... A significant deviation from 1 indicates drastic spatial changes and substantial interpolation uncertainty. In practical applications, this can be determined based on... The value is used to judge the confidence level of the interpolation result, or to assign lower weights to unreliable regions of interpolation in subsequent inversion.

[0143] Through the aforementioned depth correction and spatial interpolation, this embodiment converts the empirical Green's function at the surface into a spatially continuous propagation tensor applicable to arbitrary source depths and fault plane locations, providing a propagation operator with clear physical meaning and spatially continuous coverage for the finite fault ground motion synthesis in step 3.

[0144] Step 3: Based on the kinematic finite fault model, the fault plane in the target area is discretized into multiple sub-source units. Using the spatial continuous propagation tensor, the source time function of each sub-source unit is convolved and superimposed with the corresponding Green's function to generate a long-period ground motion waveform that includes path and site effects.

[0145] This step aims to utilize the spatially continuous propagation tensor constructed in step 2, combined with a kinematically finite fault model, to generate long-period seismic motion waveforms that incorporate path and site effects through convolutional superposition. This process organically combines the source rupture process with the wave propagation process, achieving complete seismic motion synthesis from the source to the observation point on the Earth's surface.

[0146] In this embodiment, step 3 specifically includes steps 31 to 33:

[0147] Step 31: Discretize the kinematic finite tomographic model into multiple sub-source units, each sub-source unit corresponding to a spatial location.

[0148] In practical applications, the fault planes within the target area are first discretized based on a kinematic finite fault model. The kinematic finite fault model is a mathematical description of the earthquake rupture process. It treats the fault plane as a set of multiple source elements, each of which corresponds to a spatial location and has independent source parameters (such as slip, rupture time, source time function, etc.).

[0149] In this embodiment, the fault plane is discretized into several sub-source elements using a 2km×2km grid. Each sub-source element can be considered an independent sub-source, and its source parameters are set according to the kinematic model, including the source moment tensor, rupture time, and source time function. This discretization process allows the complex continuous rupture process to be decomposed into the superposition of multiple simple sub-sources, facilitating subsequent numerical calculations.

[0150] Step 32: In the frequency domain, using the spatially continuous propagation tensor, convolve the source moment tensor of each sub-source element with the Green's function corresponding to the location of that sub-source element and sum them to obtain the seismic motion spectrum at the observation point:

[0151] ;

[0152] in, Indicates the first observation point The seismic motion spectrum of the components, Represents angular frequency. Indicates the total number of sub-source units. Indicates the first Individual source units The spectrum of the Green's function to the observation point Indicates the first The source moment tensor spectrum of each sub-source element.

[0153] In practical applications, after the fault plane is discretized, the contribution of each sub-source unit to the surface observation point is calculated using the spatially continuous propagation tensor constructed in step 2, and then superimposed and synthesized in the frequency domain. The reason for choosing to perform the calculation in the frequency domain is that the convolution of the Green's function and the source time function is simplified to a product operation in the frequency domain, which can significantly improve computational efficiency.

[0154] The physical meaning of the above formula for calculating the seismic motion spectrum is that the seismic motion at the observation point on the surface is a linear superposition of seismic waves generated by all sub-source units of the fault plane at the observation point. The contribution of each sub-source unit is determined by its own source characteristics and propagation path characteristics. By summing the contributions of all sub-source units, the complete seismic motion spectrum can be obtained.

[0155] Step 33: Perform an inverse Fourier transform on the seismic motion spectrum to obtain the long-period seismic motion waveform. The transformation formula is as follows:

[0156] ;

[0157] in, This represents the long-period ground motion waveform in the time domain. Indicates time, This represents the inverse Fourier transform operator.

[0158] Specifically, step 32 yields the seismic motion spectrum in the frequency domain, while engineering applications require a time-domain waveform. Therefore, an inverse Fourier transform (IFT) is needed to convert the seismic motion spectrum back to the time domain. Through the IFT, the frequency domain information is restored to a time series, yielding the displacement change curve at each observation point over time. This waveform already contains information about the entire process from source rupture to wave propagation and then to the site response: source characteristics are... The path effect and basin amplification effect are manifested through the Green's function. The site effect is implicit in the surface response of the Green's function.

[0159] Through the above steps, this embodiment realizes finite fault ground motion synthesis based on physical processes. Compared with traditional empirical methods, this embodiment has the following advantages: First, the physical meaning is clear: the contribution of each sub-source unit is calculated based on the actual wave propagation process, rather than empirical attenuation relationships; second, spatial continuity: the spatially continuous propagation tensor constructed in step 2 ensures a smooth transition at any position on the fault plane; third, high computational efficiency: frequency domain calculation is used, avoiding the complex operations of time domain convolution; fourth, strong scalability: the size and number of sub-source units can be adjusted as needed to balance computational accuracy and efficiency.

[0160] Step 4: Using the target response spectrum predicted by the empirical ground motion prediction equation as a constraint, establish a posterior distribution model of the source parameters; with the goal of minimizing the residual between the synthetic spectrum of the long-period ground motion waveform and the target response spectrum, perform global optimization inversion of the source parameters using the sequential Monte Carlo annealing sampling method to obtain the optimal source parameters.

[0161] This step aims to globally optimize the source parameters using a Bayesian inversion framework, constrained by the target response spectrum predicted by the empirical ground motion prediction equation. This ensures that the synthetic spectrum of the long-period ground motion waveform generated in step 3 is statistically consistent with the target response spectrum, thus obtaining the optimal source parameters. This process organically combines physical simulation with statistical constraints, solving the physical distortion problem caused by traditional methods that forcibly adjust the waveform to meet the target spectrum.

[0162] This step transforms the problem of source parameter inversion into a Bayesian inference problem. In the Bayesian framework, source parameters are treated as random variables, the distribution of which is determined by prior information (such as geological structure knowledge) and observational data (i.e., the target response spectrum).

[0163] Specifically, based on the empirical earthquake motion prediction equation ( The predicted target response spectrum For observation data, the source parameter vector Establish the posterior distribution model, which is as follows:

[0164] ;

[0165] in, The posterior probability density function representing the source parameters. This represents the source parameters to be inverted. Empirical Seismic Ground Motion Prediction Equation The predicted target response spectrum, Indicates direct proportion. This represents the natural exponential function. Indicates based on current source parameters The spectral function of the synthesized long-period seismic motion waveform. Indicates matrix transpose. Represents the covariance matrix of the observed data The inverse matrix, Represents the regularization parameter. Represents the source slip vector. Represents the covariance matrix of the sliding distribution The inverse matrix.

[0166] The expression consists of two parts: one is the likelihood term. This reflects the degree of fit between the synthesized spectrum and the target spectrum. The smaller this term is, the closer the synthesized spectrum is to the target spectrum; the second is the prior term. This reflects the physical constraints on the source slip distribution. This term penalizes the smoothness of the slip vector to prevent unreasonable and drastic changes in the inversion results.

[0167] By maximizing the posterior probability density function This allows us to obtain source parameter estimates that most closely approximate the synthetic spectrum with the target response spectrum while satisfying physical constraints.

[0168] To ensure the physical validity of the inverted source slip distribution, this embodiment employs the Von Kármán autocorrelation model to constrain the prior distribution of the source slip vector. The Von Kármán model can describe the spatial correlation and roughness characteristics of earthquake rupture surface slip, and is more physically meaningful than the traditional Laplacian smoothing constraint. Its spatial covariance function is:

[0169] ;

[0170] in, Indicates distance The covariance of the slip between two points Indicates the spatial distance between two points. Represents the sliding variance. Indicates the sliding standard deviation. This represents the power exponent, used to control the roughness of the sliding distribution. Represents the gamma function operator. This indicates the relevant scale, used to control the attenuation distance of spatial correlation. Indicates the order is The modified Bessel function.

[0171] To illustrate the impact of different prior constraints on the source slip inversion results, this embodiment further compares the source slip distribution and its uncertainties under traditional Laplacian smoothing constraints and Von Kármán autocorrelation prior constraints. Figure 2 and Figure 3 As shown, compared with the traditional Laplacian smoothing constraint, the inversion results using the Von Kármán prior constraint can more effectively focus the main slip region, generate a more spatially correlated and physically reasonable fracture distribution, and significantly reduce the uncertainty of the inversion.

[0172] Due to the posterior distribution Typically characterized by high dimensionality and nonlinearity, these parameters are difficult to solve directly using analytical methods. This embodiment employs a sequential Monte Carlo annealing sampling method, which introduces a temperature coefficient to gradually transition from a priori distribution to a posterior distribution, achieving efficient sampling of the parameter space.

[0173] In practical applications, firstly, the sequential Monte Carlo annealing sampling method is used to obtain the posterior distribution samples of the source parameters. By introducing a temperature coefficient, an intermediate posterior distribution is constructed.

[0174] ;

[0175] in, Indicates the first The intermediate posterior distribution at each temperature step Indicates the sampling step index. Indicates the current temperature coefficient. , This represents the prior distribution of the source parameters. Represents the likelihood function;

[0176] In the above construction method, when When the distance is close to 0, the intermediate posterior distribution approximates the prior distribution. Sampling is easy; with As the likelihood function gradually increases, its influence gradually strengthens; when When the value reaches 1, the intermediate distribution is the target posterior distribution. .

[0177] The iterative process of sampling is as follows:

[0178] Initialization: Several sliding field samples are generated using Von Kármán priors as initial particles;

[0179] Gradually increase the temperature coefficient: Set a temperature coefficient sequence;

[0180] Weight update: At each temperature step, calculate the importance weight of each particle. Particles with higher weights are those that are more in line with the target distribution at the current temperature.

[0181] Resampling: Particles are resampled according to their weights, eliminating particles with low weights and replicating particles with high weights;

[0182] Metropolis update: Apply random perturbations to resampled particles to maintain particle diversity;

[0183] Convergence judgment: Repeat the above process until the temperature coefficient equals 1, and obtain the final posterior sample set.

[0184] With temperature coefficient The sampling process increments from 0 to 1, starting from the prior distribution. Smooth transition to posterior distribution The final converged sample distribution corresponds to the region that minimizes the residual between the synthetic spectrum and the target response spectrum of the long-period ground motion waveform. The maximum a posteriori solution or the sample mean is selected from the converged sample set to obtain the optimal source parameters. .

[0185] After obtaining the optimal source parameters, the ground motion response spectrum within the basin area is calculated using the extracted empirical Green's function, and the error distribution between the calculated spectrum and the target response spectrum is further analyzed. To verify the improvement effect of Bayesian update on the target response spectrum constraint, this embodiment compares the error distribution of different periodic spectral functions before and after the Bayesian update, such as... Figure 4 and Figure 5 As shown. Figure 4 and Figure 5 In the graph, the horizontal axis represents the source distance, and the vertical axis represents the spectral function error; different color depths or symbols indicate the error distribution under different periodic conditions. Figure 5 It can be seen that after Bayesian update, the overall error of the spectral function under different periods is more concentrated near zero, and the dispersion is reduced, indicating that the residual between the synthesized spectrum and the target reaction spectrum is effectively reduced.

[0186] To verify the applicability of the method in this embodiment at the regional scale, Figure 6 This paper presents an empirical Green's function extracted from a basin region and the spatial distribution of errors in the response spectrum function under different periods. Figure 6 As can be seen in (a), the extracted empirical Green's function has relatively stable propagation characteristics, and is derived from... Figure 6 As can be seen from (b) to (d) in the figure, the overall error of the response spectrum function in the basin area is relatively small and the spatial distribution is relatively balanced under different periods, indicating that the method of this embodiment can achieve target spectrum matching in the basin area well.

[0187] Using the aforementioned Bayesian inversion framework, this embodiment combines physical simulation with statistical constraints to achieve global optimization of source parameters. This process not only ensures the statistical consistency between the synthetic ground motion and the target response spectrum, but also guarantees the rationality of the source parameters through physical prior constraints, laying a solid foundation for subsequent high-frequency component generation and broadband fusion.

[0188] Step 5: Based on the optimal source parameters, a short-period high-frequency ground motion component matching the long-period ground motion waveform is generated using a high-frequency radiation model; the long-period ground motion waveform and the short-period high-frequency ground motion component are fused in the frequency domain using a smoothing weighting function to generate the final ground motion time history covering the target wideband range.

[0189] This step aims to fuse the long-period ground motion waveform generated in step 3 with the high-frequency component generated based on the optimal source parameters obtained in step 4, to obtain the final ground motion time history covering a wide frequency band of 0.1-10Hz. This process solves the problem that a single method cannot accurately simulate both the propagation characteristics of long-period waves and the radiation characteristics of short-period high frequencies simultaneously.

[0190] Specifically, the long-period ground motion waveform obtained in step 3 (usually corresponding to a period range of 1-5 seconds) can accurately reflect the propagation path effect and basin amplification characteristics of seismic waves, but it is limited by the frequency band range of the empirical Green's function and is difficult to include high-frequency components (usually corresponding to a period range of 0.1-1 seconds). To obtain complete broadband ground motion, high-frequency components need to be added.

[0191] This embodiment adopts The spectral model (Brune model) generates high-frequency components. This model is a classic seismological model describing the source spectrum and can reasonably characterize the energy attenuation characteristics of high-frequency ground motions. Its spectral expression is:

[0192] ;

[0193] in, This represents the high-frequency spectrum corresponding to short-period high-frequency ground motion components. Frequency is expressed in Hz. This represents the scalar seismic moment determined based on the optimal source parameters, in N·m. Represents angular frequency. This indicates the angular frequency.

[0194] The physical meaning of this formula is: in the low-frequency range, the spectral amplitude increases with... Increase; in the high-frequency range, the spectral amplitude tends to be constant. Angular frequency. It is the dividing point between high and low frequency behavior, determined by the source scale, and the calculation formula is:

[0195] ;

[0196] in, Represents a constant. Used to adjust dimensions. This represents the shear wave velocity of the medium, expressed in m / s. This represents the stress drop, expressed in MPa.

[0197] This formula shows that the angular frequency is directly proportional to the medium shear wave velocity and inversely proportional to the cube root of the ratio of seismic moment to stress drop. The larger the source scale (i.e., The larger the stress drop, the lower the angular frequency and the relatively weaker the high-frequency components; the greater the stress drop, the higher the angular frequency and the relatively stronger the high-frequency components.

[0198] The phase of the high-frequency waveform is generated using a random phase generation method, and the amplitude is as described above. This allows us to determine and obtain short-period high-frequency ground motion components that match the long-period waveform.

[0199] In this embodiment, the long-period ground motion waveform and the short-period high-frequency ground motion component are fused in the frequency domain using a smoothing weighting function, specifically including:

[0200] Step 51: Within the preset transition frequency band, define a low-frequency weighting function and a high-frequency weighting function such that their sum is 1. The expression is as follows:

[0201] ;

[0202] ;

[0203] in, This represents the low-frequency weighting function. Represents a high-frequency weighting function. Indicates frequency, Represents the cosine function. Pi is a constant. Indicates the center frequency of the transition band. This represents half of the transition bandwidth;

[0204] Step 52: Using the low-frequency weighting function and the high-frequency weighting function, the low-frequency spectrum corresponding to the long-period ground motion waveform and the high-frequency spectrum corresponding to the short-period high-frequency ground motion component are weighted and superimposed to obtain a broadband spectrum. The superposition formula is as follows:

[0205] ;

[0206] in, Indicates broadband spectrum. This represents the low-frequency spectrum corresponding to long-period ground motion waveforms. This represents the high-frequency spectrum corresponding to the short-period high-frequency ground motion component.

[0207] Step 53: Perform an inverse Fourier transform on the broadband spectrum combined with the phase spectrum to obtain the final ground motion time history covering the target broadband range:

[0208] ;

[0209] in, Indicates the final earthquake time history, Indicates time, This represents the inverse Fourier transform operator. Represents the phase term in complex exponential form. Represents the natural constant. Represents the imaginary unit. Represents the phase spectrum, where frequency is... The function.

[0210] Specifically, after obtaining the long-period ground motion waveform and the short-period high-frequency ground motion component, the two need to be fused in the frequency domain to obtain a broadband ground motion waveform covering the entire target frequency band. Direct splicing will produce discontinuities in the transition frequency band, resulting in artificial abrupt changes in the time-domain waveform. Therefore, this embodiment uses a cosine window weighting function for smooth transition.

[0211] In practical applications, firstly, within a preset transition frequency band, a low-frequency weighting function is defined. and high-frequency weighting function This ensures that the sum of the two is 1. In this embodiment, the transition frequency band is set to 0.8-1.2Hz, and the center frequency is... Half of the transition bandwidth .

[0212] The physical meaning of the above weighting function is: within the transition frequency band, Smoothly decrease from 1 to 0. Smoothly increasing from 0 to 1; outside the transition band. and The values ​​are either 1 or 0. This smooth transition ensures that the fused spectrum remains continuous without abrupt changes at the boundary.

[0213] Then, using the aforementioned weighting function, the low-frequency spectrum corresponding to the long-period ground motion waveform is analyzed. High spectrum corresponding to short-period high-frequency ground motion components Weighted superposition yields a broadband spectrum. In the weighted superposition formula: in the low-frequency band, , The broadband spectrum is mainly contributed by long-period waveforms; in the high-frequency band, , The broadband spectrum is mainly contributed by high-frequency components; in the transition band, the two transition smoothly according to their weights.

[0214] Finally, after obtaining the broadband spectrum Then, it needs to be converted back to the time domain using the phase spectrum to obtain the final broadband ground motion time history. Phase spectrum The phase of a long-period waveform can be used (preserving the physical phase characteristics of wave propagation) or a random phase (suitable for high-frequency components). Through inverse Fourier transform, the frequency domain information is restored to a time series, yielding the displacement (or velocity, acceleration) change curve over time at each observation point. The waveform covers the target wideband range of 0.1-10Hz, with continuous energy and smooth time-frequency characteristics without abrupt changes.

[0215] To further illustrate the simulation effect and computational efficiency of the method in this embodiment, this embodiment compares and analyzes the method with the EXSIM method and the finite difference method from two aspects: the simulation results of the seismic ground motion intensity index under different source distances and the computational cost of different methods. Figure 7 As shown, Figure 7 Tables (a) to (c) show the comparative results of ground motion indices (peak ground acceleration, peak ground velocity, and pseudospectral acceleration) under different source distances. Figure 7 In this embodiment, the method is represented by Sim, the EXSIM method by EXSIM, and the finite difference method by FDM; where the horizontal axis represents the source distance and the vertical axis represents the magnitude of the corresponding ground motion index. Figure 7 As can be seen, the method in this embodiment can better reflect the attenuation law of ground motion index with distance in different distance ranges. Its result distribution is generally closer to the finite difference method, and it has better stability and physical consistency than the EXSIM method.

[0216] Figure 8 A schematic diagram comparing the computational costs of the method in this embodiment with those of the EXSIM method and the finite difference method is shown. This diagram illustrates the computational cost comparison results for a Mw 7.3 seismic scenario within the 0.1–1 Hz frequency band. Figure 8 It is evident that the finite difference method has the highest computational cost, while the EXSIM method has the lowest. Although the computational cost of the method in this embodiment is slightly higher than that of the EXSIM method, it is far lower than that of the finite difference method, achieving approximately 1250 times higher computational efficiency. These results demonstrate that the method in this embodiment possesses a significant computational efficiency advantage while maintaining good physical plausibility and simulation performance.

[0217] In summary, the ground motion generation method based on the target spectrum and Green's function provided in this embodiment achieves a unification of physical propagation mechanism and statistical spectral constraints by fusing the environmental noise Green's function with Bayesian inversion optimization. First, the empirical Green's function extracted from environmental noise records using an unsupervised clustering algorithm can accurately reflect the wave propagation characteristics and basin amplification effect within the target area, avoiding the excessive reliance on high-precision three-dimensional velocity models in traditional methods. Second, the spatially continuous propagation tensor constructed through depth correction and spatial interpolation, combined with finite fault convolution synthesis, generates long-period ground motion waveforms with clear physical meaning, ensuring the consistency of wavefield phase. Furthermore, using the target response spectrum predicted by the empirical ground motion prediction equation as a constraint, the source parameters are globally optimized through Bayesian inversion and sequential Monte Carlo sampling, making the synthesized spectrum and the target spectrum statistically accurately match, avoiding the physical distortion caused by the forced waveform adjustment to meet spectral requirements in traditional methods. Finally, the broadband ground motion time history is generated by fusing a high-frequency radiation model with a smoothing weight function, which improves computational efficiency compared to traditional three-dimensional finite difference methods while ensuring spectral compatibility. This embodiment can quickly generate physically consistent, spectrally compatible, and regionally transferable broadband ground motions under limited observation data, providing efficient and reliable technical support for seismic design of major projects and seismic hazard assessment of urban agglomerations.

[0218] Based on the above technical solution, this embodiment also provides a seismic motion generation device based on the target spectrum and Green's function, used to implement the seismic motion generation method based on the target spectrum and Green's function described in the embodiment. Please refer to [link to relevant documentation]. Figure 9 The device includes:

[0219] The Green's function extraction module is used to preprocess the continuous environmental noise records of the target area, and to evaluate the waveform similarity of the preprocessed noise records through an unsupervised clustering algorithm to extract the empirical Green's function. The empirical Green's function is used to characterize the wave propagation characteristics and basin amplification effect within the target area.

[0220] The Green's function correction and interpolation module is used to perform depth correction on the extracted empirical Green's function based on the surface wave eigenfunction, converting it from the surface response to the Green's function at the source depth; the gradient approximation method is used to perform spatial interpolation on the depth-corrected Green's function to construct a spatially continuous propagation tensor covering the fault plane in the target area.

[0221] The long-period ground motion generation module is used to discretize the fault plane in the target area into multiple sub-source units based on the kinematic finite fault model. Using the spatial continuous propagation tensor, the source time function of each sub-source unit is convolved and superimposed with the corresponding Green's function to generate a long-period ground motion waveform that includes path and site effects.

[0222] The source parameter inversion module is used to establish a posterior distribution model of the source parameters based on the target response spectrum predicted by the empirical ground motion prediction equation as a constraint; with the goal of minimizing the residual between the synthetic spectrum of the long-period ground motion waveform and the target response spectrum, the source parameters are globally optimized and inverted using the sequential Monte Carlo annealing sampling method to obtain the optimal source parameters.

[0223] The broadband fusion module is used to generate short-period high-frequency ground motion components that match the long-period ground motion waveform based on the optimal source parameters and using a high-frequency radiation model; and to fuse the long-period ground motion waveform and the short-period high-frequency ground motion components in the frequency domain through a smoothing weighting function to generate the final ground motion time history covering the target broadband range.

[0224] It is understood that since the seismic motion generation device based on the target spectrum and Green's function described in this embodiment is a device for implementing the seismic motion generation method based on the target spectrum and Green's function described in the embodiment, the device disclosed in the embodiment is relatively simple to describe because it corresponds to the method disclosed in the embodiment. For relevant parts, please refer to the description of the method, and it will not be repeated here.

Claims

1. A method for seismic generation based on a target spectrum and a Green's function, characterized by, The method includes: Step 1: Preprocess the continuous environmental noise records of the target area, and evaluate the waveform similarity of the preprocessed noise records using an unsupervised clustering algorithm to extract the empirical Green's function. The empirical Green's function is used to characterize the wave propagation characteristics and basin amplification effect within the target area. Step 2: Based on the surface wave eigenfunctions, perform depth correction on the extracted empirical Green's function to convert it from the surface response to the Green's function at the source depth; use the gradient approximation method to perform spatial interpolation on the depth-corrected Green's function to construct a spatially continuous propagation tensor covering the fault plane within the target area; Step 3: Based on the kinematic finite fault model, the fault plane in the target area is discretized into multiple sub-source units. Using the spatial continuous propagation tensor, the source time function of each sub-source unit is convolved and superimposed with the corresponding Green's function to generate a long-period ground motion waveform that includes path and site effects. Step 4: Using the target response spectrum predicted by the empirical ground motion prediction equation as a constraint, establish a posterior distribution model of the source parameters; with the goal of minimizing the residual between the synthetic spectrum of the long-period ground motion waveform and the target response spectrum, perform global optimization inversion of the source parameters using the sequential Monte Carlo annealing sampling method to obtain the optimal source parameters. Step 5: Based on the optimal source parameters, a short-period high-frequency ground motion component matching the long-period ground motion waveform is generated using a high-frequency radiation model; the long-period ground motion waveform and the short-period high-frequency ground motion component are fused in the frequency domain using a smoothing weighting function to generate the final ground motion time history covering the target wideband range.

2. The target spectrum and Green's function based seismic generation method of claim 1, wherein, Step 1 specifically includes: Step 11: Perform segmentation, bandpass filtering, spectral whitening, and amplitude normalization on the continuous environmental noise record to obtain the preprocessed noise record. The formula for amplitude normalization is: or ; in, This indicates the original noise waveform at time [time]. amplitude, This represents the normalized amplitude. Indicates taking The absolute value, express The root mean square value; Step 12: For the preprocessed noise records, calculate the cross-correlation function between each pair of stations to obtain the cross-correlation function set. The formula for calculating the cross-correlation function is: ； in, Indicates the first The station and the first Time shift between individual stations The cross-correlation function value at the location, Indicates the first Each station is at time Preprocessed noise recordings Indicates the first Each station is at time Preprocessed noise recording; Step 13: For each pair of cross-correlation functions in the set of cross-correlation functions, calculate the dynamic time warped distance. The calculation formula is as follows: ； in, Indicates the first Article and No. The dynamic time-warped distance between the cross-correlation functions. This indicates taking the minimum value. This represents the optimal curved path. This represents the index pair of corresponding points in two cross-correlation functions. Represents corresponding points in two noise waveforms and Local distance metric Indicates the first Cross-correlation function at point amplitude, Indicates the first Cross-correlation function at point The amplitude; Step 14: Construct the kernel function matrix based on the dynamic time warping distance, the expression of which is: ； in, This represents the kernel function value after mapping by kernel principal component analysis, used to measure the... Article and No. The similarity between cross-correlation functions, This represents the kernel function bandwidth, used to control the similarity decay scale. Represents the natural exponential function; Step 15: After performing kernel principal component analysis to reduce the dimensionality of the kernel function matrix, clustering is performed using a Gaussian mixture model. The cluster number corresponding to the minimum Bayesian information criterion value is selected, and the cluster with the smallest variance is chosen. The center of this cluster is used as the empirical Green's function. The expression for the Bayesian information criterion is: ； in, Represents the Bayesian information criterion value. This represents the number of clusters in the Gaussian mixture model. This indicates the number of cross-correlation functions participating in the clustering. This represents the maximum likelihood function value of the Gaussian mixture model. It represents the natural logarithm.

3. The seismic motion generation method based on target spectrum and Green's function according to claim 1, characterized in that, In step 2, depth correction is performed on the extracted empirical Green's function based on the surface wave eigenfunctions, specifically including: Based on the eigenfunctions of Love and Rayleigh waves, a depth correction factor is calculated to correct the empirical Green's function from the surface response to the focal depth. At this point, the depth-corrected Green's function is obtained, and its approximate relationship is: ； ； in, and Representing depth Green's functions for the Love wave and Rayleigh wave at that location. and These represent the Green's functions at the Earth's surface corresponding to the Love wave and the Rayleigh wave, respectively. and These represent the eigenfunction values ​​of the Love wave and Rayleigh wave at the corresponding depths, respectively. Spatial interpolation of the depth-corrected Green's function is performed using the gradient approximation method, specifically including: For locations within the fault plane that are not directly observed, a first-order gradient approximation is used to spatially interpolate the depth-corrected Green's function, constructing a spatially continuous propagation tensor, the expression of which is: ； in, Indicates the interpolation point to be determined. Green's function value at that point, Represents a known point The depth-corrected Green's function at that point, Represents the gradient operator. Represent the Green's function at a given point The gradient vector at that point, This represents the displacement vector from the known point to the interpolation point to be determined; The interpolation scaling factor is derived from the first-order gradient approximation formula and used to evaluate the accuracy of the interpolation result. Its formula is as follows: ； in, This represents the interpolation scaling factor. The closer the value is to 1, the higher the interpolation accuracy. and This represents the spatial position vector of two known points. Indicates starting from a known point to a known point Green's function value, and Indicates the component index of the Green's function. This represents its gradient vector.

4. The seismic motion generation method based on target spectrum and Green's function according to claim 1, characterized in that, Step 3 specifically includes: Step 31: Discretize the kinematic finite tomographic model into multiple sub-source units, each sub-source unit corresponding to a spatial location; Step 32: In the frequency domain, using the spatially continuous propagation tensor, convolve the source moment tensor of each sub-source element with the Green's function corresponding to the location of that sub-source element and sum them to obtain the seismic motion spectrum at the observation point: ； in, Indicates the first observation point The seismic motion spectrum of the components, Represents angular frequency. Indicates the total number of sub-source units. Indicates the first Individual source units The spectrum of the Green's function to the observation point Indicates the first The source moment tensor spectrum of each sub-source element; Step 33: Perform an inverse Fourier transform on the seismic motion spectrum to obtain the long-period seismic motion waveform. The transformation formula is as follows: ； in, This represents the long-period ground motion waveform in the time domain. Indicates time, This represents the inverse Fourier transform operator.

5. The seismic motion generation method based on target spectrum and Green's function according to claim 1, characterized in that, In step 4, the posterior distribution model is: ； in, The posterior probability density function representing the source parameters. This represents the source parameters to be inverted. Empirical Seismic Ground Motion Prediction Equation The predicted target response spectrum, Indicates direct proportion. This represents the natural exponential function. Indicates based on current source parameters The spectral function of the synthesized long-period seismic motion waveform. Indicates matrix transpose. Represents the covariance matrix of the observed data The inverse matrix, Represents the regularization parameter. Represents the source slip vector. Represents the covariance matrix of the sliding distribution The inverse matrix.

6. The seismic motion generation method based on target spectrum and Green's function according to claim 5, characterized in that, The source slip vector The prior distribution is constrained by the Von Kármán autocorrelation model, and its spatial covariance function is: ； in, Indicates distance The covariance of the slip between two points Indicates the spatial distance between two points. Represents the sliding variance. Indicates the sliding standard deviation. This represents the power exponent, used to control the roughness of the sliding distribution. Represents the gamma function operator. This indicates the relevant scale, used to control the attenuation distance of spatial correlation. Indicates the order is The modified Bessel function.

7. The seismic motion generation method based on target spectrum and Green's function according to claim 5, characterized in that, In step 4, the source parameters are globally optimized and inverted using the sequential Monte Carlo annealing sampling method, specifically including: The posterior distribution samples of the source parameters were obtained using the sequential Monte Carlo annealing sampling method. An intermediate posterior distribution was constructed by introducing a temperature coefficient. ； in, Indicates the first The intermediate posterior distribution at each temperature step Indicates the first Temperature coefficient at each temperature step , This represents the prior distribution of the source parameters. Represents the likelihood function; With temperature coefficient The sampling process increments from 0 to 1, starting from the prior distribution. Smooth transition to posterior distribution The final converged sample distribution corresponds to the region that minimizes the residual between the synthetic spectrum and the target response spectrum of the long-period ground motion waveform. The maximum a posteriori solution or the sample mean is selected from the converged sample set to obtain the optimal source parameters. .

8. The seismic motion generation method based on target spectrum and Green's function according to claim 1, characterized in that, In step 5, the high-frequency radiation model is: The spectral model, whose spectral expression is: ； in, This represents the high-frequency spectrum corresponding to short-period high-frequency ground motion components. Indicates frequency, This represents the scalar seismic moment determined based on the optimal source parameters. Represents angular frequency. Indicates the angular frequency; The turning frequency The calculation formula is: ； in, Represents a constant. Used to adjust dimensions. Indicates the shear wave velocity of the medium. This indicates stress drop.

9. The seismic motion generation method based on target spectrum and Green's function according to claim 1, characterized in that, In step 5, the long-period ground motion waveform and the short-period high-frequency ground motion component are fused in the frequency domain using a smoothing weighting function, specifically including: Step 51: Within the preset transition frequency band, define a low-frequency weighting function and a high-frequency weighting function such that their sum is 1. The expression is as follows: ； ； in, This represents the low-frequency weighting function. Represents a high-frequency weighting function. Indicates frequency, Represents the cosine function. Pi is a constant. Indicates the center frequency of the transition band. This represents half of the transition bandwidth; Step 52: Using the low-frequency weighting function and the high-frequency weighting function, the low-frequency spectrum corresponding to the long-period ground motion waveform and the high-frequency spectrum corresponding to the short-period high-frequency ground motion component are weighted and superimposed to obtain a broadband spectrum. The superposition formula is as follows: ； in, Indicates broadband spectrum. This represents the low-frequency spectrum corresponding to long-period ground motion waveforms. This represents the high-frequency spectrum corresponding to the short-period high-frequency ground motion component. Step 53: Perform an inverse Fourier transform on the broadband spectrum combined with the phase spectrum to obtain the final ground motion time history covering the target broadband range: ； in, Indicates the final earthquake time history, Indicates time, This represents the inverse Fourier transform operator. Represents the phase term in complex exponential form. Represents the natural constant. Represents the imaginary unit. Represents the phase spectrum, where frequency is... The function.

10. A seismic motion generation device based on target spectrum and Green's function, characterized in that, The apparatus for implementing the seismic motion generation method based on target spectrum and Green's function as described in any one of claims 1 to 9, the apparatus comprising: The Green's function extraction module is used to preprocess the continuous environmental noise records of the target area, and to evaluate the waveform similarity of the preprocessed noise records through an unsupervised clustering algorithm to extract the empirical Green's function. The empirical Green's function is used to characterize the wave propagation characteristics and basin amplification effect within the target area. The Green's function correction and interpolation module is used to perform depth correction on the extracted empirical Green's function based on the surface wave eigenfunction, converting it from the surface response to the Green's function at the source depth; the gradient approximation method is used to perform spatial interpolation on the depth-corrected Green's function to construct a spatially continuous propagation tensor covering the fault plane in the target area. The long-period ground motion generation module is used to discretize the fault plane in the target area into multiple sub-source units based on the kinematic finite fault model. Using the spatial continuous propagation tensor, the source time function of each sub-source unit is convolved and superimposed with the corresponding Green's function to generate a long-period ground motion waveform that includes path and site effects. The source parameter inversion module is used to establish a posterior distribution model of the source parameters based on the target response spectrum predicted by the empirical ground motion prediction equation as a constraint; with the goal of minimizing the residual between the synthetic spectrum of the long-period ground motion waveform and the target response spectrum, the source parameters are globally optimized and inverted using the sequential Monte Carlo annealing sampling method to obtain the optimal source parameters. The broadband fusion module is used to generate short-period high-frequency ground motion components that match the long-period ground motion waveform based on the optimal source parameters and using a high-frequency radiation model; and to fuse the long-period ground motion waveform and the short-period high-frequency ground motion components in the frequency domain through a smoothing weighting function to generate the final ground motion time history covering the target broadband range.

Images