Physical constraint earthquake risk analysis method for super-deep overburden site

Abstract

The application discloses a kind of super-deep overburden site physical constraint seismic risk analysis method, it is related to earthquake monitoring technical field, and the method comprises: according to potential seismogenic structure determination seismic activity parameter, constructs strong earthquake scenario space and extracts high-fidelity training scenario and to be extended scenario;Perform regional local coupling physical simulation to obtain high-fidelity training sample library;Utilize simulation sample and real earthquake event record training based on physical constraint depth learning agent model, by double source constraint acceptance agent generated extended time history, to untrusted sample trigger supplementary high-fidelity simulation;Only allow time history to enter probability seismic risk integration after acceptance, adopt multi-proposal mechanism mixed density to calculate probability correction weight, obtain annual exceedance probability, site risk curve and design ground motion parameter, and prohibit to place exceedance probability zero when support domain is missing.The application improves the integrity and accuracy of regional seismic risk assessment.

Description

Technical Field

[0001] This application relates to the field of earthquake monitoring technology, specifically to a method for analyzing the seismic hazard of sites with physical constraints in ultra-deep overburden layers. Background Technology

[0002] Traditional probabilistic seismic hazard analysis (PSHA) relies heavily on historical earthquake observation data and empirical attenuation relationships of ground motion parameters. In areas with complex topographic and geological conditions (such as those with deep overburden), empirical formulas cannot accurately reflect the local ground motion amplification effect. Meanwhile, while physical-mechanism-based numerical simulations of ground motion can reflect site effects, their computational cost is extremely high, and the limited simulation scenarios cannot meet the data requirements of probabilistic analysis for massive earthquake scenarios, thus limiting the completeness and accuracy of regional seismic hazard assessments. Summary of the Invention

[0003] To address the technical problem that traditional probabilistic seismic hazard analysis (PSHA) relies heavily on historical earthquake observation data and empirical ground motion parameter attenuation relationships, which limits the completeness and accuracy of regional seismic hazard assessment, this application provides a physical constraint seismic hazard analysis method for ultra-deep overburden sites.

[0004] The specific technical solution adopted is as follows:

[0005] S10: Determine seismic activity parameters based on potential seismogenic structures, construct a strong earthquake scenario space and divide it into scenario hierarchical units, extract a first number of high-fidelity training scenarios for physical simulation and a second number of scenarios to be expanded according to a proposal mechanism with known density, establish a unique probability sample identifier and solution state field for each strong earthquake scenario-field point, obtain real earthquake events of the target site and surrounding stations, and associate source parameters, site parameters and event identifiers.

[0006] S20, for representative high-fidelity training scenarios selected from the spatial stratification of strong earthquake scenarios, performs coupled physical simulation from the regional scale to the local site, obtains the bedrock-surface paired three-component time history, regional propagation state and site nonlinear state, and forms a high-fidelity training sample; performs unified processing on the high-fidelity spectral training sample and the three-component time history of real earthquake events respectively, and extracts the simulated earthquake field features and real event features.

[0007] S30 uses a high-fidelity training sample library and real earthquake event records to train a dual-source physical constraint deep learning surrogate model to replace the complete physical simulation of other strong earthquake scenarios; the surrogate model is configured to learn the mapping from source-path-site conditions to bedrock wavefield and surface response, and introduces a real event bias calibration mechanism.

[0008] S40, during the training process of the dual-source physical constraint deep learning surrogate model, simulated seismic field constraints and real seismic event constraints are established, which are respectively used as loss terms or acceptance measures, and the calibration energy and frequency band change are limited to not exceeding the conditional quantile threshold determined by the real events.

[0009] S50: In the scenario to be expanded, the accepted dual-source physical constraint deep learning proxy model is invoked to generate candidate samples of expanded ground motion time history. The sample credibility is determined based on the physical sensitivity distance of the scenario to be expanded relative to the high-fidelity training scenario and the real earthquake event condition domain, the cognitive discrepancy between multiple proxy models, and the deviation between the sample-by-sample simulated earthquake field residual and the real event calibration. Only when the credibility is not lower than the preset threshold and the hard physical constraints are passed, is it allowed to enter the hazard integration; otherwise, a supplementary high-fidelity physical simulation is triggered, and the process returns to S20 to update the proxy model.

[0010] S60 reads the solution status of each strong earthquake scenario-field point, and only allows time histories in the proxy effective or high-fidelity effective state to enter the probabilistic seismic hazard integral; determines the probability correction weight by the ratio of the target condition density in the layer to the mixed proposed density in the layer, and calculates the annual exceedance rate, site hazard curve and design ground motion parameters.

[0011] This application has, but is not limited to, the following technical effects:

[0012] This application takes the site hazard of ultra-deep overburden layers as the final output, placing the annual occurrence rate of potential seismic sources, scenario probability quality, time histories of earthquakes of magnitude 10,000, and exceedance events of seismic ground motion intensity indices into the same probability chain. This directly generates the annual exceedance rate of acceleration response spectra with different periods, site hazard curves, and exceedance probabilities for the design reference period. This improves the completeness and accuracy of regional seismic hazard assessment. Attached Figure Description

[0013] Figure 1 This is a schematic diagram of the closed-loop process of a small number of complete spectral element working conditions in this application—dual-source physical constraint deep learning—ten-thousand-level time history—site hazard integral;

[0014] Figure 2 A schematic diagram illustrating the work-equivalent coupling between the regional-scale bedrock propagation domain and the local fine-scale site domain;

[0015] Figure 3 A schematic diagram of the process for managing the cyclic elastoplastic update and trial calculation of local spectral elements at Gauss-Lobatto-Legendre integration points and the submission status.

[0016] Figure 4 A schematic diagram of a dual-source physical projection conditional adversarial deep learning process consisting of a condition generator, a hard physics projection operator, a simulated seismic field condition discriminator, and a real seismic event condition discriminator.

[0017] Figure 5 A schematic diagram of the process for directly generating the predicted bedrock-surface three-component time history and high-fidelity backoff without pre-performing the calculation of the spectral elements of the scenario to be expanded;

[0018] Figure 6 A schematic diagram illustrating the simulated seismic field constraints, real seismic event constraints, dual-source conditional domain, and sample credibility acceptance process;

[0019] Figure 7 A schematic diagram of the active high-fidelity compensation process driven by the expected reduction of the risk variance of unit complete spectral element node kernels;

[0020] Figure 8 A schematic diagram of the risk integration process for intra-layer hybridization proposals with probability mass conservation;

[0021] Figure 9 A schematic diagram of the process for checking the effective sample size at the tail and the generalized Pareto extension of physical adjacency;

[0022] Figure 10 A schematic diagram of the two-layer propagation process involving outer-layer cognitive uncertainty and inner-layer hierarchical weighted sampling;

[0023] Figure 11 A schematic diagram illustrating the separation of uses, version dependencies, and hierarchical acceptance processes for real event records, spectral element training records, and hazardous sample records;

[0024] Figure 12 This is a schematic diagram of the module composition of the dual-source physical constraint deep learning seismic hazard analysis system of this application;

[0025] Figure 13 This diagram illustrates the state transitions in scenario-based problem-solving, including proxy reasoning, high-fidelity supplementary calculation, atomic substitution, and risk re-integration, under the unique probability sample identification. Detailed Implementation

[0026] It should be understood that the specific embodiments described herein are merely illustrative of this application and are not intended to limit this application.

[0027] This application provides a method for analyzing the seismic hazard of ultra-deep overburden sites under physical constraints. In the first embodiment of this method, referring to... Figure 1 The methods include:

[0028] Step S10: Determine seismic activity parameters based on potential seismogenic structures, construct a strong earthquake scenario space and divide it into scenario hierarchical units, extract a first number of high-fidelity training scenarios for physical simulation and a second number of scenarios to be expanded according to a proposal mechanism with known density, establish a unique probability sample identifier and solution state field for each strong earthquake scenario-field point, obtain real earthquake events of the target site and surrounding stations, and associate source parameters, site parameters and event identifiers.

[0029] This application constructs a data stream with three separate uses.

[0030] The first category is simulated seismic field data streams: dozens to hundreds of representative working conditions are selected from the spatial stratification of strong earthquake scenarios, and complete regional-local spectral element calculations are performed to provide regional bedrock wavefields, surface three-component time histories, propagation travel times, energy envelopes, bedrock-surface transmission relationships, and overburden nonlinear states.

[0031] The second category is real earthquake event data streams: instrument response, baseline, coordinates and effective frequency band processing are performed on the measured records of the target site and surrounding stations to provide measured time histories, response spectra, Arias intensity, significant duration and multi-station spatial coherence constraints; real event records are only used for training and calibration and are not used as hazard probability samples.

[0032] The third category is the hazard probability sample data stream: based on the annual occurrence rate of potential seismic sources and the scenario probability model, tens of thousands of scenarios to be expanded are extracted. The bedrock and surface three-component time histories are directly generated by the dual-source physical constraint deep learning proxy model, and the probability mass conservation hazard integral is entered according to the actual proposed density.

[0033] In the high-fidelity spectral element scenario, the regional-scale bedrock propagation domain covers the main propagation path from the potential source to the target site, directly calculating the source excitation, path propagation, and the three-component wavefield of the target site's bedrock multi-field point. The local fine-grained site domain is established around the target site, explicitly including topography, the overburden-bedrock interface, ultra-deep overburden, and the bedrock buffer zone, responsible for calculating the cyclic nonlinear site response. The regional and local models are coupled through work-equivalent boundaries. Complete spectral element calculations are only performed on representative training scenarios and actively supplemented scenarios. For scenarios awaiting expansion that have passed dual-source physical acceptance, the surrogate model's bedrock propagation branch directly predicts the bedrock wavefield, and the overburden site response branch then obtains the surface time history from the predicted bedrock wavefield; regional or local spectral element calculations for this scenario are not performed beforehand.

[0034] In practice, the annual occurrence rate and magnitude-frequency relationship (such as the GR relationship) of potential seismogenic structures (e.g., the Medog Fault Zone and the Qiongduo River Fault Zone in the lower reaches of the Yarlung Tsangpo River) are obtained from historical earthquake catalogs. The strong earthquake scenario space is jointly determined by magnitude, rupture location and direction, focal mechanism, stress drop, rupture velocity, propagation path parameters, random rupture realization, target field coordinates, and ultra-deep overburden parameters. This scenario space is divided into multiple scenario stratification units (e.g., stratified by magnitude interval or rupture location), each unit possessing target probability mass (allocated by the annual occurrence rate of potential source areas). Following a proposal mechanism with known density (e.g., uniform sampling or importance sampling), a first number (tens to hundreds) of representative scenarios are extracted as high-fidelity training scenarios for subsequent physical simulations; a second number (much larger than the first, such as tens of thousands) of scenarios to be expanded are extracted for surrogate model inference. For each strong earthquake scenario-site point, a unique probability sample identifier (such as a hash ID) is established, and a solution status field is set (to be solved, surrogate candidate, surrogate valid, queued for supplementary calculation, high-fidelity valid, integrated). Simultaneously, real earthquake event records from the target site and surrounding stations (such as acceleration time histories recorded by strong-motion meters) are acquired, and associated with source parameters (magnitude, location, etc.), site parameters (Vs30, overburden thickness, etc.), and event identifiers. It is clarified that this record is only used for training and calibration data and not as a hazard probability sample (i.e., it does not have a proposed density).

[0035] Step S20: For representative high-fidelity training scenarios selected from the spatial stratification of strong earthquake scenarios, perform coupled physical simulation from the regional scale to the local site to obtain the bedrock-surface paired three-component time history, regional propagation state and site nonlinear state, and form a high-fidelity training sample; perform unified processing on the high-fidelity spectral training sample and the three-component time history of real earthquake events to extract the simulated earthquake field features and real event features.

[0036] In practical implementation, for each high-fidelity training scenario, coupled physical simulations from the regional scale to the local site are performed. The regional-scale bedrock propagation domain covers the main propagation path from the potential source to the target site. The spectral element method is used to solve the elastic wave equation, directly obtaining the three-component wave field of the multi-field point of the source-path-target site bedrock. The local fine-grained site domain is established around the target site, explicitly including the topographic surface, ultra-deep overburden, overburden-bedrock interface, and bedrock buffer zone. The regional and local models are coupled through work-equivalent boundaries: a kinematic projection matrix Pu and a surface force projection matrix Pt are constructed, such that the local boundary displacement and velocity are equal to the product of Pu and the regional boundary displacement and velocity, and the boundary surface force is equal to the product of Pt and the regional boundary surface force. Within the local fine-grained site domain, pressure-related multi-yield surface cyclic elastoplastic updates are performed at the GLL integration points. This constitutive model dynamically reconstructs the yield surface based on effective stress, cyclic history, and degradation variables, using local substep partitioning and return mapping. Joint acceptance is performed through state integration error, discrete energy balance, yield surface consistency, and non-negative plastic dissipation. This yields the bedrock three-component time histories, surface three-component time histories, regional propagation status, and constitutive status of ultra-deep overburden associated with the same strong earthquake scenario identifier, forming a high-fidelity training sample library. Subsequently, the high-fidelity samples and real earthquake event records are processed uniformly to extract features of the simulated seismic field (propagation travel time, energy envelope, bedrock-surface complex transfer matrix, etc.) and features of real events (response spectrum, Arias intensity, significant duration, etc.).

[0037] Specifically, in S20, the work-equivalent projection of the regional-local coupling boundary includes: constructing a regional-scale bedrock propagation domain covering the potential source to the target site and a local fine-grained site domain containing ultra-deep overburden, topographic surface, and bedrock interface; performing wave calculations at the regional scale to obtain the bedrock wave field, and transferring the regional boundary velocity and surface force to the local site domain through work-equivalent projection; and performing pressure-related multi-yield surface cyclic elastoplastic constitutive model updates in the local site domain to obtain high-fidelity samples associated with the same scenario identifier.

[0038] In S20, the mesh of the local fine field region is adaptively optimized according to the comprehensive discrete error index. The comprehensive discrete error index takes into account wavelength resolution, interface curvature, material-state gradient and coupling index, and determines the orientation scale and order of each element with the goal of minimizing the spatiotemporal computation cost.

[0039] In S20, the cyclic elastoplastic update adopts local substep partitioning and return mapping, and is jointly verified by state integral error, discrete energy balance, yield surface consistency and non-negative plastic dissipation, and the stable time step is determined based on the consistent tangent modulus that has passed the verification.

[0040] In S20, the simulated seismic field characteristics include a high-fidelity three-component complex transfer matrix, which is obtained through multi-time-window weighted estimation and satisfies conjugate symmetry, positive semi-definite residual power spectrum, and negative time energy constraints. The pressure-dependent multi-yield-surface cyclic elastoplastic constitutive model dynamically reconstructs the yield surface based on effective stress, cyclic history, and degradation variables, and updates the back stress and hardening modulus after local substep acceptance.

[0041] In S20, the integration point states are managed using a submission state buffer and a trial state buffer. The global time step is verified by the global discrete energy residual and plastic dissipation, and the data is written to the submission state only when the verification is successful. The regional-scale physical simulation results are reused through incident feature quantity decomposition and projection interpolation, with the boundary energy flux residual and complex spectrum propagation error used as coupling verification indicators.

[0042] Reference Figure 1 and Figure 2 The regional-scale bedrock propagation domain covers the potential seismogenic fault, main propagation path, and target site area, and may include three-dimensional velocity structure, density structure, attenuation structure, and topography. The regional propagation domain is discretized using spectral elements and the source is applied according to the source kinematics or dynamic parameters, directly advancing the displacement, velocity, and stress wave fields within the regional model. A local fine-grained site domain is set around the target site, with its bottom and sides extending a certain distance into the stable bedrock to form a bedrock buffer zone. Within the local domain, the topographic surface z_s(x,y) and the overburden-bedrock interface z_b(x,y) are explicitly reconstructed. The overburden thickness is:

[0043]

[0044] For regions where H(x,y)>0, normalized depth coordinates ξ=[z_s(x,y)−z] / H(x,y) are used to describe different thickness locations, ensuring that material parameters obtained from drilling, testing, and inversion are consistently mapped along the stratigraphic plane. The integration point state can be written as q_g=[H,ξ,V_s,V_p,ρ,p′_0,γ_r,θ_c], where θ_c is the set of cyclic elastoplastic parameters. Geometric tolerances and topological stitching are employed at the pinch-out points of the overburden layer to avoid negative thickness, isolated thin-sheet elements, or interface normal jumps.

[0045] For the regional and locally coupled boundaries, establish the kinematic projection matrix P_u and the surface force projection matrix P_t, respectively, and satisfy:

[0046]

[0047] Where W^L and W^R are the spectral element integral weight matrices of the local and regional boundaries, respectively. Under the unified method, the local boundary displacement, velocity, and surface force are projected by P_u and P_t, respectively. The instantaneous boundary energy flux residual is:

[0048]

[0049] When r_Γ exceeds the limit, priority is given to encrypting the projection points of the coupling boundary or shortening the interpolation interval of the local time axis of the region time axis, without eliminating the error by fitting the input waveform with arbitrary amplitude.

[0050] At each boundary integration point, an acoustic tensor A_n and an impedance matrix Z=(ρA_n)^(1 / 2) are constructed using the normal n pointing from the region to the local region. The characteristic components entering the local domain are extracted from the boundary velocity v^R and surface force t^R obtained directly from the region spectral elements.

[0051]

[0052]

[0053] The incident characteristic quantities at the regional time step are transferred to the local boundary through windowed sinc interpolation and spatial projection. The cumulative incident energies E_in^R and E_in^{L,0} before and after interpolation are calculated, and both are checked to be positive. If either cumulative incident energy is not greater than zero, it indicates that at least one of the unified normal, incident eigenvalue decomposition, or time-space projection has failed. In this case, amplitude correction is not performed, but the normal is checked and the time interpolation or boundary integration points are densified. Only when both are positive, a_Γ=[(E_in^R+ε_Γ) / (E_in^{L,0}+ε_Γ)]^(1 / 2) is used for single energy consistency correction, and the cumulative energy error and target band complex spectrum error are checked again.

[0054] The local spectral mesh comprehensively considers the shortest shear wavelength, topographic and overburden-bedrock dual-interface curvature, material-state gradient, and coupling error. A comprehensive index η_e = max(η_λ,e,η_Γ,e,η_q,e,η_c,e) is defined for the spectral element e, minimizing the spatiotemporal computational cost under the conditions of η_e ≤ 1 and limited scale variations and spectral order jumps between adjacent spectral elements. The stable step size for explicit time advancement uses the minimum distance between adjacent Gauss-Lobatto-Legendre nodes in physical coordinates, rather than performing a second external multiplication for spectral order correction.

[0055]

[0056] This expression already incorporates the effects of spectral order, element distortion, and the non-uniform distribution of Gauss-Lobatto-Legendre nodes. For a uniform tangent modulus that may not possess principal symmetry, a directional acoustic tensor is first constructed and its symmetric portion is taken; the minimum eigenvalue is checked for positive within the discrete propagation direction set on a unit sphere, and then c_alg,e is determined using the maximum eigenvalue. A non-positive minimum eigenvalue indicates a strong ellipticity failure, and the corresponding local substep must not be submitted.

[0057] The scenario space is organized by potential source region z and intra-layer element h. Each scenario stores magnitude, rupture geometry, mechanism, stress drop, rupture velocity, random realization number, target condition density p_tar(c|h), actual proposed mechanism identifier, and corresponding proposed density q_k(c|h). The strong earthquake scenario identifier ID_s is obtained by hashing the above source parameters and random realization number, and remains unchanged for the same scenario in regional propagation, local site, high-low fidelity pairing, and multi-field output.

[0058] Reference Figure 3 The effective stress, deviatoric stress, back stress of each yield surface, plastic strain, active yield surface, cyclic memory variables, and cumulative plastic dissipation are stored at the Gauss-Lobatto-Legendre integration points of the overburden layer. The reference shear modulus and reference bulk modulus are determined from the in-situ shear wave velocity and compressive wave velocity.

[0059]

[0060]

[0061] The cumulative plastic shear memory variable κ, the bounded cyclic degradation variable d_κ, and the normalized effective pressure Π are respectively:

[0062]

[0063]

[0064]

[0065] The current small strain modulus, bulk modulus, and reference shear strain are dynamically updated based on Π and d_κ. For example:

[0066]

[0067]

[0068] Where p_min and p_t are normal stress parameters that ensure the stability of the effective pressure ratio domain, b_κ is a positive cyclic memory evolution coefficient, n_G, n_K, and n_γ are pressure exponents calibrated by cyclic tests under different confining pressures, and d_G, d_K, and d_γ are between 0 and 1 and calibrated by cyclic degradation tests of modulus and reference shear strain. Parameter calibration must simultaneously satisfy G_max>0, K_max>0, γ_r>0, monotonic skeleton curves, and non-negative plastic dissipation. When the pressure range exceeds the test coverage area, an out-of-domain parameter flag is triggered instead of unconstrained extrapolation.

[0069] The target shear skeleton curve can be written as:

[0070]

[0071] At the logarithmically distributed shear strain control points, the monotonically increasing target shear stress is obtained using constrained least squares, and the tangential shear modulus is constrained to be non-increasing from the inside to the outside, and the outer layer strength does not exceed c+p′tanφ. The l-th yield function is defined as F_l=√(3 / 2)||s−α_l||_F−R_l, so that R_l=√3τ_l is consistent with the von Mises equivalent stress under simple shear conditions; the piecewise hardening modulus is H_l=G_maxG_t,l / (G_max−G_t,l+ε_G), so that the elastoplastic series tangent recovers the target G_t,l. After each accepted local substep, the yield surface is reconstructed according to the new effective pressure and cyclic memory, and the back stress is mapped according to the normalized position of the current stress in the old yield surface.

[0072]

[0073] This mapping avoids unloaded plastic jumps caused by simply changing the yield surface radius. After reconstruction, the nesting order of the yield surfaces, the relative position of the current stress, and non-negative plastic dissipation also need to be checked.

[0074] In a global time step, the strain path is parameterized as ε(s) = ε_n + sΔε. The initial local substep proportion is limited by both the reference shear strain proportion and the degree to which the test stress crosses the yield surface. For candidate substeps, one full-step update z^(1) and two half-step updates z^(2) are calculated respectively, and the difference between the two forms the state error r_z; and the stress work ΔW, free energy increment Δψ, plastic dissipation ΔD_p, energy residual r_E, yield surface consistency residual r_F and negative dissipation residual r_D are calculated:

[0075]

[0076]

[0077] When χ≤1, two half-step states are accepted; when χ>1, candidate substeps are rejected and the substep is reduced without being written to historical variables. All local states are placed in the trial buffer, and the trial internal forces are assembled only after the local acceptance is passed.

[0078] At the end of the global time step, the kinetic energy increment ΔE_kin, recoverable internal energy increment ΔE_ψ, plastic dissipation ΔW_p, energy carried away by the absorbing boundary and linear decay ΔW_a, and input work ΔW_ext from the source and coupling boundary are calculated to form the global discrete energy residual:

[0079]

[0080] Only when all integration points pass local acceptance, r_G does not exceed the limit, and the total plastic dissipation is non-negative, is the trial buffer written entirely into the commit buffer; otherwise, displacement, velocity, stress, and all historical variables are rolled back to reduce the global time step recalculation. The consistent tangent modulus C_alg of the integration points after acceptance is used to construct the directional acoustic tensor Q(n), and the algorithm wave velocity and the next global stable time step are determined from its largest eigenvalue.

[0081] Each high-fidelity scenario also outputs the total plastic dissipation W_p, the volume fraction of plastic influence φ_p, the high quantile value of the cyclic memory variable κ_95, the high quantile value of the local substep number m_95, the maximum energy residual, and the volume ratio of each active yield surface, forming a constitutive state fingerprint bound to the scenario identifier.

[0082] Step S30: A dual-source physical constraint deep learning surrogate model is trained using a high-fidelity training sample library and real earthquake event records to replace the complete physical simulation of other strong earthquake scenarios. The surrogate model is configured to learn the mapping from source-path-site conditions to bedrock wavefield and surface response, and a real event bias calibration mechanism is introduced.

[0083] In S30, the dual-source physical constraint deep learning proxy model uses a travel-time anchored frequency domain neural operator to learn the bedrock wave field and outputs the site transfer operator on the Hermitian positive definite amplitude-unitary phase complex matrix manifold, and superimposes the local scattering residuals to form the surface response.

[0084] In S30, the surrogate model is implemented through a condition generator and a physical projection operator. The physical projection operator projects the generated results onto a set of physically feasible conditions that satisfy time travel, conjugate symmetry, energy, and semidefinite constraints.

[0085] Reference Figure 4 High-fidelity spectral training samples and real earthquake event records employ a usage-separate preprocessing workflow. For high-fidelity spectral training samples, complete strong earthquake scenario identifiers, regional-local grid versions, material and constitutive parameter versions, bedrock-surface multi-field three-component time histories, and regional propagation and overburden constitutive states are retained; the same time origin, coordinate basis vector, and overlapping time window are used for all field points within the same scenario. For real earthquake event records, instrument responses are first removed, baseline correction, anomalous impulse identification, coordinate direction unification, and effective frequency band marking are performed. Then, event source parameters, station coordinates, site parameters, and processing versions are bound to the time histories. Measured records without reliable source parameters, station coordinates, or effective frequency bands are not included in the training. Real event records do not have hazard proposal density and are only used as training and calibration evidence.

[0086] The dual-source physically constrained deep learning surrogate model includes a bedrock propagation branch, a cover layer site response branch, and a real event deviation calibration branch. The bedrock propagation branch employs a travel-time anchored three-component frequency-domain neural operator. The source condition encoder receives magnitude, mechanism, rupture geometry, stress drop, rupture velocity, and rupture inhomogeneity codes associated with a random seed; the path condition encoder receives regional velocity, density, path characteristics of attenuation structures, target field point coordinates, and frequency codes. The reference first wave travel time τ_m is determined by the regional medium model and follows:

[0087]

[0088] Output the predicted bedrock three-component complex spectrum at the m-th field point. Here, S_θ is the three-component source complex spectrum, and G_θ is the 3x3 propagation operator. Explicitly extracting exp(-i2πfτ_m) allows the network to primarily learn the propagation amplitude, phase deviation, and component transformation after travel time correction. For negative frequencies, a conjugate extension of the positive frequency results is used, and real-valued vectors are output for zero-frequency and Nyquist frequency points. The pre-arrival energy ratio is defined as:

[0089]

[0090] Δτ is determined by the travel time error quantile of the high-fidelity spectral calibration scenario. Predicting the bedrock wavefield also requires satisfying the source-path conditions, including the bandgap energy ratio, cumulative energy envelope, and the quantile range of the three-component complex coherence at multiple field points. These constraints allow the bedrock propagation branch to directly replace the regional spectral calculation of the scenario to be expanded, rather than requiring the spectral bedrock wavefield of that scenario to be obtained first.

[0091] The overburden site response branch takes as input the predicted bedrock wavefield, ultra-deep overburden thickness and velocity structure, topography and overburden-bedrock interface parameters, bedrock input intensity, and frequency encoding. The 3x3 complex transfer matrix H^HF is inverted from the high-fidelity spectral sample, and the low-fidelity matrix H^LF is calculated from the co-originating one-dimensional three-component nonlinear overburden column. Polar decomposition H=QP is performed respectively, where Q is the unitary phase factor and P is the Hermitian positive definite amplitude factor, and the following definitions are provided:

[0092]

[0093]

[0094] The amplitude correction head and phase correction head output the original complex matrices Z_A and Z_Omega, respectively. Hermitian amplitude difference and oblique Hermitian phase difference are obtained using A_hat=(Z_A+Z_A^H) / 2 and Omega_hat=(Z_Omega−Z_Omega^H) / 2, respectively. The multi-fidelity field transfer matrix is ​​then reconstructed via matrix exponentiation. Therefore, the positive definiteness of the amplitude and the unitary nature of the phase are guaranteed by the network structure, without relying on post-processing removal of anomalous frequencies.

[0095] The multi-fidelity site transfer matrices for each field point are combined into a block diagonal operator H_script, and the principal physical response of the surface is Y_hat^(s,phy)=H_scriptY_hat^b. The conditional residual branch outputs the residual mean mu_eta, complex low-rank factor B_eta, non-negative diagonal residual d_eta, and constitutive state gate G based on the local scattering in the high-fidelity samples that is not explained by the transfer matrix. A positive semi-definite residual covariance is constructed according to C_eta=B_etaB_eta^H+diag(d_eta). The total energy of the candidate local scattering residuals must not exceed the upper limit of the residual energy ratio of the high-fidelity calibration samples in adjacent conditional domains.

[0096] The real-event bias calibration branch does not replace the simulated seismic field main branch. Instead, it utilizes the conditional differences between the measured surface time history and the simulated constrained surface response under the same source-station conditions to output the observation calibration transfer increment and observation residuals. The observation calibration transfer increment is still Hermitian-oblique Hermitian projection, and the observation residuals are limited by energy, effective frequency band, and pre-time energy. When there is no adjacent real-event support, the calibration gate automatically decays to unit transfer and zero observation residuals, causing the model to degenerate into an output constrained only by the simulated seismic field, prohibiting unbounded extrapolation corrections for conditions without measured coverage.

[0097] The training samples of real events are used to construct event condition vectors based on their event source parameters, station coordinates, propagation path characteristics, and site parameters. These condition vectors are then input into a condition generator to obtain simulated constrained surface responses. These responses are then compared with measured surface time histories under the same event and station identifier to calculate observational calibration quantities. This process does not require registering real event records as hazard probability samples, nor does it require pre-performing full-area-local spectral calculations for each real event. The conditional neighborhood of a real event is determined by both the physical sensitivity distance and effective frequency band on the event calibration set. If either of these conditions is not met, observational calibration is disabled; adjacency is not determined solely by the names or geometric distances of neighboring stations.

[0098] As an improved conditional adversarial implementation, the bedrock propagation branch, the overburden site response branch, and the real event deviation calibration branch are all used as the condition generator G_θ, and a simulated seismic field condition discriminator D_sim and a real seismic event condition discriminator D_obs are set up. The parameters of the two discriminators are not shared: D_sim determines whether the predicted bedrock-surface paired three-component time-frequency tensor is consistent with the high-fidelity spectral element simulated field under given source-path-field-overburden conditions; D_obs determines whether the predicted surface three-component time-frequency tensor is consistent with the conditional distribution of the real event record under given event source-station-site conditions and effective frequency band mask. This avoids the confusion of conditional meaning caused by directly mixing the internal physical state of the simulated field and the measured record into a single discriminator.

[0099] The conditional generation results are not directly used in adversarial training as the original network output. Instead, they are first processed by the physical projection operator P_phy. P_phy sequentially performs positive and negative frequency conjugate projection, reference travel time and pre-time energy projection, Hermitian positive definite amplitude-unitary phase projection of the site transfer matrix, positive semi-definite projection of the local scattering residual covariance, and upper bound projection of the residual energy. Therefore, the discriminator compares the candidate time histories that satisfy the hard physical structure with the true time histories of the corresponding data source. The adversarial gradient can only adjust the conditional distribution within the set of physically feasible outputs and cannot exchange pre-time spurious energy, non-causal transfer matrices, or non-positive semi-definite residual covariance for a discriminant score.

[0100] During training, simulated samples are grouped into D_sim,tr and D_sim,cal based on the complete strong earthquake scenario identifier, while real records are grouped into D_obs,tr and D_obs,cal based on the real event identifier. Different field points, components, and time windows of the same simulated scenario or the same real event must not cross sets. Simulation fitting terms include bedrock and surface complex spectra, response spectra, Arias intensity, and significant duration errors; simulation physics terms include travel time, pre-arrival energy, energy envelope, spatial coherence, complex matrix projection, and residual energy; real event fitting terms calculate the conditional time-frequency distribution, response spectrum, Arias intensity, and significant duration differences within the effective station-band mask; real event physics terms include multi-station spatial coherence, pre-arrival energy, observation calibration energy, and band change. The total loss is:

[0101]

[0102] For r∈{sim,obs}, the r-th conditional discriminator adopts a conditional Wasserstein objective with gradient penalty:

[0103]

[0104] Where y_bar_r = ξy_r + (1−ξ)y_hat_r, and ξ follows a uniform distribution on [0,1]. The objective of the condition generator is:

[0105]

[0106] Where g_obs∈[0,1] is the real event conditional domain gating; when there is no adjacent event-station-site condition support or the effective frequency band is empty, g_obs=0, thus closing the real event adversarial branch. Each supervised loss term is first divided by the median value of the initial training phase to make it dimensionless, and λ_1 to λ_5, α_sim, α_obs, and λ_gp are non-negative weights. During training, D_sim, D_obs, and G_θ are updated alternately, and augmented Lagrangian updates are used to update the constraint multipliers for travel time, causality, energy, semidefiniteness, and spatial coherence. The condition generator enters the inference phase only when both the simulated calibration set and the real event calibration set pass the specified error quantile threshold; when inferring scenarios of tens of thousands, the two conditional discriminators are not loaded, and only G_θ and P_phy are run. The discriminator output is only used for training and acceptance and is not included in the hazard probability weights.

[0107] In the executable parameter determination process, the number of positive frequency discrete points K is first determined based on the unified time history length, sampling interval, and common effective frequency band of the two types of data. The real and imaginary parts of the three-component complex spectrum are used as channels respectively. The source, path, field point, and overburden scalars are robustly standardized using only the median and interquartile range of the simulation training set. The real event data retains the same physical dimensions and coordinate definitions but is saved separately in the event domain standardized version. The number of frequency domain neural operator blocks, hidden channel width, residual covariance rank, and capacity of the two discriminators for the bedrock propagation branch are selected from a pre-given finite candidate set in ascending order. The smallest structure that simultaneously satisfies the simulation calibration error, real event calibration error, hard physics over-limit rate, and inference resource constraints is selected. If no candidate structure meets the conditions, the model must not enter S50. Training batches are composed of complete simulated scenario identifiers and real event identifiers. The condition generator is updated once after each of the two discriminators is updated. The number of updates, learning rate, gradient penalty coefficient, dimensionless loss baseline, augmented Lagrange multiplier, and stopping rounds are all stored as model version fields, so that the same training data version and random seed can reproduce training and inference results.

[0108] Step S40: During the training process of the dual-source physical constraint deep learning proxy model, simulated seismic field constraints and real seismic event constraints are established, which are used as loss terms or acceptance measures, respectively, and the calibration energy and frequency band change are limited to not exceeding the conditional quantile threshold determined by leaving out real events.

[0109] As an example, simulated seismic field constraints are used to calculate propagation travel time residuals, energy ratios before arrival, positive and negative frequency conjugate residuals, regional propagation energy envelope residuals, bedrock-surface transfer causal residuals, changes in physical projection of the site transfer operator, energy ratios of local scattering residuals, and three-component spatial coherence errors for spectral element samples. Real seismic event constraints are used to calculate differences in time-frequency distribution, response spectrum, Arias intensity, significant duration, and spatial coherence differences across multiple stations for measured samples. These constraints are used as loss terms or output projection constraints, respectively, and limit the calibration energy and frequency band changes of the real event deviation calibration branch to no more than the conditional quantile thresholds determined by the reserved real events. During training, the training set and calibration set are divided according to complete simulated scenario identifiers and real event identifiers; different field points and components of the same scenario do not cross sets. A surrogate model can only be used to replace spectral element calculations when the errors and conditional domain coverage on the calibration set meet the allowable thresholds.

[0110] In step S40, a dual-source physical projection condition adversarial training framework is adopted, a simulated seismic field condition discriminator and a real seismic event condition discriminator are set up, and the physical constraint limit is incorporated into the training objective through augmented Lagrange multipliers.

[0111] Step S50: Invoke the accepted dual-source physical constraint deep learning proxy model for the scenario to be expanded, generate candidate samples of expanded ground motion time history, and determine the sample credibility based on the physical sensitivity distance of the scenario to be expanded relative to the high-fidelity training scenario and the real earthquake event condition domain, the cognitive discrepancy between multiple proxy models, and the deviation between the sample-by-sample simulated earthquake field residual and the real event calibration. Only when the credibility is not lower than the preset threshold and the hard physical constraints are passed, is it allowed to enter the hazard integration; otherwise, trigger the supplementary high-fidelity physical simulation and return to S20 to update the proxy model.

[0112] Reference Figure 5 For the scenarios to be expanded required for hazard analysis, the system inputs the strong earthquake scenario conditions, target field conditions, and random realization numbers. The bedrock propagation branch directly outputs the predicted three-component bedrock wavefield, the overburden site response branch outputs the simulated constrained surface response, and the real event deviation calibration branch outputs bounded calibration results when supported by measured conditions. The entire inference process does not pre-execute regional-scale or local fine-grained spectral element calculations for the scenario to be expanded. Therefore, after training and calibration with dozens to hundreds of complete spectral element cases, three-component earthquake time histories can be generated in batches for combinations of strong earthquake scenarios and field points in the tens of thousands range.

[0113] Each predicted sample stores the strong earthquake scenario identifier, field identifier, model version, random seed, predicted bedrock and surface three-component time histories, simulated seismic field acceptance parameters, real event calibration deviation, and cognitive discrepancy. Random implementation coding only characterizes rupture heterogeneity or local random scattering already defined in S10 and with known proposal mechanisms; it does not generate additional hazardous samples without probabilistic identities. The same strong earthquake scenario identifier and random seed should yield reproducible output.

[0114] A sample-by-sample hard acceptance test must include at least the following: the real-valuedness and baseline stability of the predicted bedrock and surface time histories; reference travel time and energy before arrival; band energy and cumulative energy envelope; three-component complex coherence between field points; conjugate symmetry, causality, singular values, and energy range of the site transfer matrix; positive semidefiniteness of the local scattering covariance and upper limit of residual energy; and energy, band, and conditional domain support for the real event bias calibration. If any hard acceptance test fails, the time history must not be included in the hazard integral.

[0115] In one executable implementation, representative complete spectral cases are selected in layers based on potential source region, magnitude range, rupture direction, propagation path, overburden thickness, bedrock input intensity, and degree of plastic mobilization, with the number of cases ranging from dozens to hundreds. Each case outputs three-component time histories for multiple field points. The surrogate model inference generates no fewer than ten thousand strong earthquake scenario-field-point surface three-component time histories. The quantity is used to illustrate the computational scale relationship and does not constitute a fixed lower limit; the actual quantity is jointly determined by the target year exceedance rate accuracy, the tail effective sample size, and the computational budget.

[0116] The verification process recorded the bedrock and surface time-frequency errors, response spectrum errors, travel time errors, energy, and spatial coherence errors on the simulated calibration set, as well as the response spectrum, Arias intensity, significant duration, and multi-station spatial coherence errors on the real-event calibration set. The reduction in computational resources was calculated based on the estimated node core time if all scenarios to be expanded performed full-region-local spectral element calculations, minus the initial high-fidelity spectral element calculations, active supplementary calculations, and actual node core times for training and inference. This resource indicator is only used to record technical effects and does not replace physical acceptance.

[0117] In S50, the dual-source applicability domain and sample credibility of the scenario to be expanded are determined according to the following method: the source parameters, propagation path parameters, target field coordinates, overburden thickness and velocity structure, predicted bedrock input intensity, singular values ​​and group delays of the predicted site transfer matrix, and predicted ground motion intensity indices are used to form conditional embeddings; the simulation conditional domain and the real event conditional domain are established by using the conditional embeddings of each high-fidelity strong earthquake scenario in the simulation training set and the conditional embeddings of each real earthquake event-station record in the event training set, respectively; and the physical sensitivity metric matrix is ​​constructed based on the Jacobi norm of the predicted bedrock wavefield, site transfer matrix, and ground motion intensity indices according to each input direction.

[0118] Reference Figure 6 The conditional vectors are robustly standardized and a conditional embedding phi(c) is formed, which includes the source, path, field point, overburden, predicted bedrock strength, singular values ​​and group delay of the predicted transfer matrix, and predicted ground motion intensity index. Simulated conditional domains and real-event conditional domains are established using simulated training scenarios and real-event station records, respectively. A physical sensitivity metric M(c) is constructed based on the Jacobian norm of the predicted bedrock wavefield, field transfer matrix, and intensity index for each input direction. The nearest physical sensitivity distances d_sim(c) and d_obs(c) from the scenario to be extended to the two data domains are calculated, respectively.

[0119] Multiple surrogate models, resampled and trained according to simulated scenario labels and real event labels, are used to calculate the inter-model deviations for predicting bedrock complex spectrum, site transfer matrix, and ground motion intensity indices, yielding the cognitive discrepancy u_epi. The sample-by-sample simulated seismic field acceptance quantity and real event calibration exceedance quantity are then denoted as a_sim and a_obs, respectively. Using the finite sample ranking of the two calibration sets, d_sim, d_obs, u_epi, a_sim, and a_obs are converted into rank scores between 0 and 1, and the minimum value is taken to obtain the sample confidence level.

[0120]

[0121] Sample confidence is used only to determine whether the proxy time history can be used and whether to trigger recalculation; it does not represent the probability of a strong earthquake scenario occurring and is not included in the hazard weighting. If the real event domain does not cover the conditions to be expanded, the real event deviation calibration branch is turned off and a conservative event domain score is assigned; if the overall confidence is below the threshold, the time history cannot be used.

[0122] In S50, supplementary high-fidelity scenario batches that return to S20 for full region-local spectral element calculation are selected according to the following method: scenarios to be expanded are formed from scenario hierarchical units whose sample confidence is below the threshold, whose simulated conditional domain distance or real event conditional domain distance exceeds the corresponding threshold, whose hard physical constraints are not passed, and whose tail-weighted effective sample size at the target intensity threshold is insufficient; for the current high-fidelity dataset, the current variance of the previous year's exceedance rate of the design intensity threshold set is calculated using multiple cognitive uncertainty branches, and virtual high-fidelity updates are performed on each candidate scenario in the candidate batch using the posterior prediction results of the surrogate model.

[0123] Reference Figure 7The candidate scenarios to be expanded, including those with low confidence, failure in hard physics acceptance, out-of-domain conditions due to dual sources, and insufficient effective sample size in the tail of hazards, are grouped into a supplementary calculation candidate set. Full spectral element calculations have not yet been performed on these candidate scenarios. For candidate batch B, the expected reduction in annual exceedance rate variance U(B) on the design intensity threshold set is estimated using a surrogate model with posterior virtual updates. The full spectral element cost simultaneously accounts for regional propagation, local fluctuations, and integral point constitutive updates.

[0124]

[0125] Each non-negative cost coefficient is calibrated using the actual node core time of historical high-performance computing tasks. The batch selection objective is:

[0126]

[0127] It also satisfies single-round node kernel time budget, physical sensitive distance diversity, and critical tail layered cover constraints. Selecting the scenario and returning to S20 executes full-region—local spectral calculation, adding results to update bedrock propagation branches, overburden site response branches, and simulation condition domains; real event records and their event identifiers remain unchanged. High-fidelity spectral time histories replace the surrogate time histories of the same probability sample, without adding new probability samples.

[0128] High-fidelity recalculation is stopped only when the changes in the hazard curve, the annual exceedance rate coefficient of variation, the tail-weighted effective sample size, the simulation calibration error, the real event calibration error, and the decrease in the unit node kernel time variance of the remaining candidates all meet the stopping criteria. This stopping rule ensures that the number of spectral element scenarios is jointly controlled by hazard results and dual-source physical evidence, rather than assuming that a fixed training size can cover all scenarios.

[0129] Step S60: Read the solution status of each strong earthquake scenario-field point, and only allow time histories in the proxy effective or high-fidelity effective state to enter the probabilistic seismic hazard integral; determine the probability correction weight by the ratio of the target condition density in the layer to the mixed proposed density in the layer, and calculate the annual exceedance rate, site hazard curve and design ground motion parameters.

[0130] In S60, the probability mass conservation multi-proposal hazard integration is performed as follows: For the h-th scenario stratification unit in the z-th potential source area, the samples are divided into initial stratification proposal samples, scenario proposal samples to be expanded, and support domain supplementary proposal samples according to the actual mechanism for generating strong earthquake scenario condition vectors, and the k-th actual proposal mechanism is denoted as k; For the support domain supplementary proposal mechanism established based on previous samples and acceptance results, the historical information, condition proposal density, sample allocation number, and random seed on which the proposal mechanism is based are frozen and versioned before sampling. Subsequent proxy model updates, high-fidelity supplementary calculation results, or hazard results shall not revert the already executed condition proposal density.

[0131] As an example, refer to Figure 8 Within each scenario stratification unit, hazard samples are divided according to the proposal mechanism of the actual generated strong earthquake scenario condition vector c, including the initial stratification proposal mechanism, the scenario proposal mechanism to be expanded, and the support domain supplementary proposal mechanism. Real earthquake event records do not have hazard proposal density and are not included in this sample set. S50 only selects complete spectral element supplementation objects from strong earthquake scenarios that have already been generated and have actual proposal density. Active supplementation only changes the time history solution source and does not constitute a new proposal mechanism. When the same strong earthquake scenario is regressed from proxy inference to complete high-fidelity spectral element calculation, the expanded ground motion time history is replaced with the high-fidelity spectral element time history, and the replaced time history is invalidated. The scenario condition vector is counted as only one probability sample. The proposal mechanism mixing coefficient is determined by the actual number of selected samples for each mechanism. The mixing proposal density and probability correction weights are:

[0132]

[0133]

[0134]

[0135]

[0136] The conditional exceedance probability, annual exceedance rate, and design reference period exceedance probability are as follows:

[0137]

[0138]

[0139]

[0140] The exceedance probability for the design reference period in the above formula adopts the implementation method of the earthquake occurrence process in each potential source area following the Poisson model. When using the time-dependent recurrence model, the update process model, or other occurrence processes calibrated by seismic activity data, the annual exceedance rate ν(y) is still obtained by integrating the exceedance probability under the above scenario conditions, while the exceedance probability for the design reference period is determined by the survival function of the corresponding occurrence process, and the exponential form is not mandatory.

[0141] If q_mix ≤ q_min in the subdomain where p_tar > 0, then the scenario's hierarchical unit has a missing support domain or the proposed density value is indistinguishable. It is necessary to return to S10 to establish a supplementary proposal mechanism with known density in the missing subdomain and extract the supplementary strong earthquake scenario, first through S50 surrogate inference. If the dual-source applicable domain or sample-by-sample physical acceptance fails, then S20 is used to complete the full region-local high-fidelity spectral element calculation. Regularization constants must not be added to q_mix to mask the missing data. Passing the S50 acceptance means that the time history has sufficient dual-source physical consistency, but it does not change the target probability quality it represents. When regressing from surrogate inference to high-fidelity spectral element calculation, the strong earthquake scenario condition vector is not re-extracted, thus retaining the original proposal mechanism and its density.

[0142] Support domain supplementation proposal mechanisms can be adaptively established based on the previous round's missing support domains, tail effective sample size, and physical acceptance results. However, before each round of sampling, the historical information, conditional proposal density, sample allocation number, and random seed upon which the proposal mechanism is based must be frozen and versioned. Intra-layer hybrid proposal density uses only the versioned conditional proposal density that has been actually implemented and can be evaluated across the entire target support domain. Subsequent surrogate model updates, high-fidelity rollbacks, or hazardous results must not retroactively modify the proposal mechanisms and densities of already sampled scenarios. If the density of a supplementation mechanism cannot be determined during sample generation, the generated scenario must not be included in the integration as a hazardous probability sample.

[0143] Reference Figure 9 The tail-weighted effective sample size and the maximum normalized weight at the target threshold y are respectively:

[0144]

[0145]

[0146] The weighted empirical exceedance probability is only effective when the support domain is covered, N_eff meets the target, and κ_w does not exceed the limit. If there are no samples above the target threshold, the event is not interpreted as "probability equal to zero". When the lower candidate starting threshold u_h has enough effective samples and the shape parameter is stable with the threshold, a weighted generalized Pareto distribution can be jointly fitted among the physical adjacency scenario hierarchical units. The likelihood term uses the normalized tail weights within the hierarchical unit ω̃_i,h=ω_i,hI(IM_i>u_h) / [Σ_jω_j,hI(IM_j>u_h)+ε_ω] to avoid the number of samples in different hierarchical units changing the graph regularization intensity only through the likelihood scale; the adjacency weights are given by the Gaussian kernel of the physical sensitive distance, and graph smoothing is applied to the shape and logarithmic scale parameters, but each hierarchical unit still needs to satisfy the parameter support and calibration interval. The extension when y>u_h is:

[0147]

[0148] When ξ_h approaches 0, the exponential limit is used, and continuity at u_h is guaranteed. Extension is prohibited if there is no valid initial threshold, and the sample is returned.

[0149] In each round, after the newly added complete regional-local high-fidelity spectral results are added to the high-fidelity spectral training sample library, the extended ground motion time history with the same strong earthquake scenario identifier and field point identifier is replaced by the high-fidelity spectral time history. The replaced extended ground motion time history and the high-fidelity spectral time history are prohibited from being simultaneously included in the hazard integral as two probability samples. The newly added results are used to sequentially update the bedrock propagation branch, the overburden site response branch, the physical constraint allowable threshold of the simulated seismic field, and the dual-source applicable domain. The real earthquake event training set and its event identifiers remain unchanged. The real event deviation calibration branch is only used again if it passes the event calibration set acceptance after retraining. The dual-source physical constraint deep learning surrogate model is used to replace the regional-scale bedrock propagation domain and the local fine-scale site domain complete spectral calculation for the scenario to be extended, rather than just replacing the local overburden solution. The complete simulated scenario identifier and the real event identifier are then grouped and cross-reserved again for recalculation. Determine the simulated seismic field error, real event error, sample credibility, and target condition domain coverage status; update the mixed proposal density, probability correction weight, conditional exceedance probability, annual exceedance rate, and its uncertainty within the layer based on the selected sample records corresponding to the actual proposal mechanism; stop supplementing high-fidelity spectral element calculations only when the following conditions are met simultaneously for multiple consecutive rounds: the relative change in annual exceedance rate within the design intensity threshold set is not greater than the hazard convergence threshold; the annual exceedance rate coefficient of variation is not greater than the uncertainty threshold; the minimum tail-weighted effective sample size of all scenario stratification units with non-zero target probability quality within the design intensity threshold set is not less than; the specified error quantile values ​​of both the simulated calibration set and the real event calibration set are not greater than the corresponding allowable values; and the expected reduction in the annual exceedance rate variance per unit node core time within the remaining node core time budget for any candidate batch is not greater than the preset stop threshold; if any condition is not met, continue executing S20 to S60.

[0150] Establish interconnected and separately used real event records, high-fidelity spectral element training records, and hazard sample records. Real event records should include at least the real event identifier, source parameters, station identifiers and coordinates, site parameters, instrument response processing information, effective frequency band, measured three-component time histories and their response spectra, Arias intensity, significant duration, and spatial coherence characteristics, and should be marked as being used only for training and calibration and not participating in probability weight calculations. High-fidelity spectral element training records should include at least the strong earthquake scenario identifier, regional and local spectral element mesh versions, material and constitutive parameter versions, random seed, bedrock-surface three-component time histories, regional propagation status, overburden constitutive status, numerical acceptance parameters, high-performance computing task identifiers, and actual... When the node core is in operation, the hazardous sample record should include at least the potential source area and scenario layer unit identifier, strong earthquake scenario identifier, field point identifier, target condition density, actual proposal mechanism, actual proposal density, time history solution source, dual-source physical constraint surrogate model version, sample confidence, effective three-component time history and ground motion intensity index; when the same strong earthquake scenario-field point is regressed from surrogate inference to high-fidelity spectral element calculation, its actual proposal mechanism and actual proposal density should be retained, the extended ground motion time history should be replaced with high-fidelity spectral element time history and the replaced time history should be marked as invalid; when any simulation data, real event data, network parameters, grid, constitutive or calibration version changes, all dependent fields should be invalidated and the corresponding steps should be re-executed.

[0151] A tiered acceptance process is implemented for unified records: If the instrument response, baseline, coordinate orientation, or effective frequency band of the real event record is unqualified, return to S10 to reprocess the real event data; if the grid resolution, region-local coupling, cyclic elastoplastic integral, negative plastic dissipation, discrete energy, or stable time step of the high-fidelity spectral element training record is unqualified, or if the simulated seismic field characteristics or real event feature extraction is unqualified, return to S20; if the travel time, causality, conjugate symmetry, energy, spatial coherence, complex matrix physical projection, response spectrum, Arias intensity, or significant duration error of the bedrock propagation branch, overburden site response branch, or real event deviation calibration branch is unqualified... Return to S30 or S40; if the dual-source applicability domain, cognitive discrepancy, or sample credibility of the scenario to be expanded is not up to standard, return to S50 and trigger S20 to supplement the complete region—local high-fidelity spectral element scenario; if the intra-layer mixed proposal support domain, probability correction weight, tail-weighted effective sample size, or tail stability is not up to standard, return to S60 and trigger S10 to supplement the strong earthquake scenario with known proposal density, or trigger S50 and S20 to supplement the tail key high-fidelity scenario, depending on the reason for the missing information; each failure invalidates the result fields and version identifiers after the failed step, prohibits the use of invalid results to continue calculation, and prohibits replacing missing conditional exceedance probability or danger results with zero values.

[0152] The annual exceedance rate, site hazard curve, and design ground motion parameters are used as the results of the seismic hazard analysis for the target site. Outputs include peak ground acceleration, peak velocity, Arias intensity, the mean annual exceedance rate, annual exceedance rate quantile curve, site hazard curve, and exceedance probability for the design reference period corresponding to acceleration response spectra with 5% to 95% significant duration and multiple periods. Design ground motion parameters, rare ground motion parameters, and extremely rare ground motion parameters are inversely solved based on the combination of the given design reference period and exceedance probability. For multiple sites, the same strong earthquake scenario identifier, outer cognitive uncertainty branch, and inner weighted sampling index are used to simultaneously calculate the hazard of each site. The correlation matrix of intensity indices between sites, joint exceedance probability, and spatial hazard distribution of peak ground acceleration or multiple-period acceleration response spectra are output. The annual exceedance rate and its variance are decomposed according to the potential source area, magnitude interval, rupture location and direction, nonlinear state of ultra-deep overburden, time history solution source, and uncertainty source. Before output, all hazard curves are checked for non-negative monotonicity, and quantile curves are checked for non-crossover.

[0153] As an example, refer to Figure 10 The outer cognitive uncertainty branch covers seismic activity, rupture, regional medium, cyclic elastoplastic constitutive model, bedrock propagation pathway, site transfer matrix, local scattering residuals, and real event bias calibration parameters. The inner layer maintains scenario stratification and performs multiplier resampling within each outer branch, adhering to the actual proposal mechanism structure. For samples within the same stratified unit and proposal mechanism, commutative nonnegative multipliers a_i are generated and normalized so that the sum of multipliers within a group equals the number of samples in that group. Each resampling uses a_iω_i as a one-time combined weight to calculate the conditional exceedance probability, avoiding repeated weighting by multiplying by ω_i after sampling. The sum of sampling variance and cognitive variance is used as the total annual exceedance rate variance. A nonnegative monotonic projection is performed on the annual exceedance rate over discrete intensity thresholds, and non-crossing constraints are applied simultaneously to multiple quantile curves.

[0154] Reference Figure 11 The versioned data storage area establishes separate records for real events, high-fidelity spectral element training, and hazard sample records. Real event records store event identifiers, station information, instrument processing versions, measured three-component time histories and their intensity and spatial coherence characteristics, and are marked as being used only for training and calibration. High-fidelity spectral element training records store strong earthquake scenario identifiers, regional and local grids, constitutive and material versions, bedrock-surface time histories, regional propagation and overburden status, numerical acceptance parameters, supercomputing task identifiers, and node kernel times. Hazard sample records store scenario layering units, target condition density, actual proposal mechanisms and densities, time history solution sources, surrogate model versions, credibility, and effectively expressed time histories and intensity indices. These three types of records cannot be directly merged due to filenames or similar site locations; they can only be associated through explicit event identifiers, strong earthquake scenario identifiers, site identifiers, and version dependencies.

[0155] The three types of records employ a data access relationship that is segregated by purpose. The real event record area does not have target condition density, proposal mechanism, or proposal density fields. Real earthquake event processing nodes have write permissions, while the training and inference engines and acceptance controllers have read-only permissions. The site hazard analysis and result output engine does not have an interface to read real event records as probability samples. The high-fidelity spectral element training record area only opens records that have passed regional-local coupling, cyclic elastoplasticity, and discrete energy acceptance as training labels. The hazard sample record area must simultaneously possess a unique probability sample identifier, target condition density, actual proposal density, solution source, and solution status. Only records in a proxy-valid or high-fidelity-valid state can be read by the hazard analysis and result output engine.

[0156] Acceptance testing is performed in the order of S10 to S60. Issues related to instrument response to real events, baseline, coordinates, or effective frequency bands should be returned to S10; issues related to spectral grids, coupling, constitutive model, negative dissipation, global energy, or time step should be returned to S20; issues related to dual-source characteristics should be returned to S25; issues related to bedrock propagation, site response, real-event calibration, travel time, causality, energy, spatial coherence, or response spectrum should be returned to S30 or S40; issues related to dual-source applicability domains, cognitive discrepancies, or credibility should be returned to S50 and trigger a supplementary calculation in S20; issues related to proposed support domains should be returned to S10 to establish a supplementary proposal mechanism; issues related to tail evidence should be returned to S50 and S20 to supplement key operating conditions. All result fields after a failed step are invalid and must not be replaced with zero values ​​to represent missing hazard evidence.

[0157] Reference Figure 12 This system consists of a probabilistic scenario and sample manager, real earthquake event processing nodes, full-region-local spectral element calculation nodes, a dual-source physical projection condition adversarial training and inference engine, a physical acceptance and active supplementary calculation controller, a site hazard analysis and result output engine, a scenario solution status register, and a versioned data storage area. The computing system has at least one or more processors, memory, data communication interfaces, deep learning training and inference computing interfaces, and high-performance computing task scheduling and node time measurement interfaces. Each node can be deployed on the same high-performance computing and artificial intelligence computing platform or connected via a high-speed network.

[0158] The method and system are mapped item by item according to the same data object: S10 is executed by the probabilistic scenario and sample manager and the real earthquake event processing node; S20 is executed by the full region-local spectral element calculation node; S25, S30, and S40 are executed collaboratively by the versioned data storage area, the training and inference engine, and the acceptance controller; S50 is executed by the proxy inference task queue, the training and inference engine, the scenario solution state register, and the active supplementary calculation controller; S60 is executed by the site hazard analysis and result output engine. A unique probabilistic sample identifier runs through all the above nodes. The target condition density and the actual proposal density are passed along with this identifier and do not change with the source of the solution. Therefore, the system does not rename the method steps one by one, but implements the method through controlled interfaces, state registers, task queues, and result failure relationships.

[0159] The probabilistic scenario and sample manager receives annual occurrence rates of potential seismic source areas, magnitude and rupture probability models, regional medium models, target field points, and ultra-deep overburden parameters. It establishes scenario hierarchical units and extracts high-fidelity training scenarios and scenarios to be expanded, generating a unique probability sample identifier for each strong earthquake scenario-field point. Representative high-fidelity training scenarios enter the complete spectral element task queue, while a larger number of scenarios to be expanded enter the proxy inference task queue. Both queues store the same format of probability sample identifiers and proposed density fields, but different task types and solution source fields. The real earthquake event processing node receives station measurement records and event-station metadata, completes instrument response, baseline, coordinate, and effective frequency band processing, and provides it to the dual-source physical projection condition adversarial training and inference engine through a read-only training interface.

[0160] The full-region-local spectral element computation node only receives representative training scenarios and active supplementary calculation scenarios. It performs local cyclic nonlinear calculations of regional bedrock propagation and ultra-deep overburden, writing the three-component time histories, physical states, numerical acceptance quantities, and actual node kernel times of bedrock-surface multi-field points into a versioned data storage area. The dual-source physical projection conditional adversarial training and inference engine reads two types of constraint data with different purposes from high-fidelity spectral element training records and real event records. During training, the hard physics projection operator projects the output of the condition generator to the physically feasible set, and then the simulated seismic field condition discriminator and the real seismic event condition discriminator perform dual-condition Wasserstein adversarial training respectively. When inferring scenarios to be expanded, the two condition discriminators are turned off, and no corresponding tasks are sent to the full-region spectral element computation node. Only the condition generator and hard physics projection operator are run, directly outputting the predicted bedrock wavefield, predicted surface time histories, dual-source physical acceptance quantities, and cognitive divergence quantities.

[0161] The physical acceptance and active recalculation controller only transmits a valid time history table that has passed the simulated seismic field constraints, real event constraints, dual-source condition domain, and credibility acceptance to the site hazard analysis and result output engine. If the acceptance fails, or if the expected reduction in the annual exceedance rate variance on the design intensity threshold set reaches the recalculation trigger threshold, a recalculation task is generated based on the expected reduction in the annual exceedance rate variance per unit full spectral node. This task is sent to the full region-local spectral calculation node via the high-performance computing task scheduling interface. After the recalculation is completed, the proxy time history for the same scenario-site is replaced with a high-fidelity spectral time history.

[0162] The system is configured with separate high-performance computing (HPC) resource pools and artificial intelligence (AI) computing resource pools. Full spectral element task queues are scheduled only to the HPC resource pool, while proxy inference task queues are scheduled only to the AI ​​computing resource pool. The HPC task scheduling and node core time measurement interface records task identifier, number of nodes, number of processor cores, runtime, mesh version, constitutive version, and actual node core time; the deep learning training and inference computing interface records model version, input version, inference batch, and per-sample physical acceptance quantity. The active compensation controller determines task priority based on the ratio of the expected reduction in the variance of the previous year's exceedance rate of the candidate batch against the design strength threshold set to the estimated full spectral element node core time. This ensures that the system prioritizes limited HPC resources for scenarios outside the dual-source condition domain, those failing hard physical constraints, or those with critical tail hazards, rather than initiating spectral element tasks for every scenario in the tens of thousands.

[0163] The site hazard analysis and results output engine establishes an intra-layer mixed proposal density based on the actual proposal mechanism rather than the time history solution source. It completes probabilistic mass conservation integrals, tailing, and uncertainty propagation, and outputs peak ground acceleration, peak velocity, Arias intensity, significant duration, annual exceedance rate of acceleration response spectra for different periods, site hazard curves, quantile curves, and exceedance probabilities for different design reference periods. It also outputs design ground motion parameters, rare and extremely rare ground motion parameters, spatial hazard distribution at multiple sites, and hazard contributions decomposed by source, path, overburden state, time history solution source, and uncertainty source. All output curves must be non-negative and monotonically non-increasing with the intensity threshold, and quantile curves must not cross.

[0164] Reference Figure 13The scenario-solving state register stores a current state for each unique probability sample identifier. New samples are in a pending state; surrogate inference output is written to the surrogate candidate state; after passing dual-source physical acceptance, it transitions to a surrogate valid state; if acceptance fails, it is outside the condition domain, or it is selected as a critical supplementary calculation object in the tail of the hazard, it transitions to a supplementary calculation queue state; after the complete region-local spectral result passes numerical acceptance, it transitions to a high-fidelity valid state; after the hazard engine completes one round of integration, it transitions to an integrated state. If an integrated surrogate sample is subsequently selected as a supplementary calculation object, the hazard result dependent on that surrogate timeline is first marked as invalid, and then enters the supplementary calculation queue state. When a high-fidelity result is returned, the original proposal mechanism and proposal density are retained in the same state update transaction, the original surrogate valid identifier is revoked, the high-fidelity valid identifier is written, and hazard re-integration at the affected intensity threshold is triggered. At any given time, the same unique probability sample identifier cannot have two valid timelines that can be read by the hazard engine simultaneously.

[0165] In this embodiment, the system includes a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of any of the methods described above.

[0166] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or system that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or system. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or system that includes that element.

[0167] The sequence numbers of the embodiments in this application are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments.

[0168] The above are merely preferred embodiments of this application and do not limit the scope of this application. Any equivalent structural or procedural transformations made based on the description and drawings of this application, or direct or indirect applications in other related technical fields, are similarly included within the scope of protection of this application.

Claims

1. A method for analyzing the seismic hazard of sites with physical constraints in ultra-deep overburden layers, characterized in that, The method includes: S10: Determine seismic activity parameters based on potential seismogenic structures, construct a strong earthquake scenario space and divide it into scenario hierarchical units, extract a first number of high-fidelity training scenarios for physical simulation and a second number of scenarios to be expanded according to a proposal mechanism with known density, establish a unique probability sample identifier and solution state field for each strong earthquake scenario-field point, obtain real earthquake events of the target site and surrounding stations, and associate source parameters, site parameters and event identifiers. S20, for representative high-fidelity training scenarios selected from the spatial stratification of the strong earthquake scenario, perform coupled physical simulation from the regional scale to the local site to obtain the bedrock-surface paired three-component time history, regional propagation state and site nonlinear state, and form a high-fidelity training sample; perform unified processing on the high-fidelity spectral training sample and the three-component time history of the real earthquake event respectively to extract the simulated earthquake field features and real event features. S30, a dual-source physical constraint deep learning surrogate model is trained using the high-fidelity training sample library and real earthquake event records to replace the complete physical simulation of other strong earthquake scenarios; the surrogate model is configured to learn the mapping from source-path-site conditions to bedrock wavefield and surface response, and a real event bias calibration mechanism is introduced. S40, during the training process of the dual-source physical constraint deep learning proxy model, simulated seismic field constraints and real seismic event constraints are established, which are respectively used as loss terms or acceptance measures, and the calibration energy and frequency band change are limited to not exceeding the conditional quantile threshold determined by leaving out real events. S50: For the scenario to be expanded, call the accepted dual-source physical constraint deep learning proxy model to generate extended ground motion time history candidate samples. The sample credibility is determined based on the physical sensitivity distance of the scenario to be expanded relative to the high-fidelity training scenario and the real earthquake event condition domain, the cognitive discrepancy between multiple proxy models, and the deviation between the sample-by-sample simulated earthquake field residual and the real event calibration. Only when the credibility is not lower than the preset threshold and the hard physical constraints are passed, is it allowed to enter the hazard integration; otherwise, a supplementary high-fidelity physical simulation is triggered, and the process returns to S20 to update the proxy model. S60 reads the solution status of each strong earthquake scenario-field point, and only allows time histories in the proxy effective or high-fidelity effective state to enter the probabilistic seismic hazard integral; determines the probability correction weight by the ratio of the target condition density in the layer to the mixed proposed density in the layer, and calculates the annual exceedance rate, site hazard curve and design ground motion parameters.

2. The method for analyzing the seismic hazard of ultra-deep overburden sites with physical constraints as described in claim 1, characterized in that, In step S20, the work-equivalent projection of the regional-local coupling boundary specifically includes: constructing a regional-scale bedrock propagation domain covering the potential seismic source to the target site and a local fine-grained site domain containing ultra-deep overburden, topographic surface, and bedrock interface; performing wave calculations at the regional scale to obtain the bedrock wave field, and transferring the regional boundary velocity and surface force to the local site domain through work-equivalent projection; and performing pressure-related multi-yield surface cyclic elastoplastic constitutive model updates in the local site domain to obtain high-fidelity samples associated with the same scenario identifier.

3. The method for analyzing the seismic hazard of ultra-deep overburden sites based on physical constraints according to claim 2, characterized in that, In step S20, the mesh of the local fine field region is adaptively optimized according to a comprehensive discrete error index. The comprehensive discrete error index takes into account wavelength resolution, interface curvature, material-state gradient, and coupling index, and determines the directional scale and order of each element with the goal of minimizing spatiotemporal computational cost. The cyclic elastoplastic update adopts local substep partitioning and return mapping, and is jointly verified by state integral error, discrete energy balance, yield surface consistency, and non-negative plastic dissipation. The stable time step is determined based on the consistent tangent modulus that has passed the verification. The simulated seismic field characteristics include a high-fidelity three-component complex transfer matrix, which is obtained by multi-time window weighted estimation and satisfies conjugate symmetry, positive semi-definite residual power spectrum, and negative time energy constraints.

4. The method for analyzing the seismic hazard of ultra-deep overburden sites with physical constraints according to claim 3, characterized in that, The pressure-dependent multi-yield-surface cyclic elastoplastic constitutive model dynamically reconstructs the yield surface based on the effective stress, cyclic history, and degradation variables, and updates the back stress and hardening modulus after local substep acceptance.

5. The method for analyzing the seismic hazard of ultra-deep overburden sites with physical constraints according to claim 4, characterized in that, In S20, the status of integration points is managed by a submission status buffer and a trial status buffer. The global time step is verified by global discrete energy residual and plastic dissipation, and the submission status is written only when the verification is successful.

6. The method for analyzing the seismic hazard of ultra-deep overburden sites with physical constraints according to claim 1, characterized in that, In S30, the dual-source physical constraint deep learning proxy model uses a travel-time anchored frequency domain neural operator to learn the bedrock wave field, and outputs a site transfer operator on a Hermitian positive definite amplitude-unitary phase complex matrix manifold, and superimposes local scattering residuals to form the surface response.

7. The method for analyzing the seismic hazard of ultra-deep overburden sites with physical constraints according to claim 6, characterized in that, In S30, the surrogate model is implemented through a condition generator and a physical projection operator. The physical projection operator projects the generated result onto a set of physically feasible conditions that satisfy time travel, conjugate symmetry, energy, and positive semidefinite constraints.

8. The method for analyzing the seismic hazard of ultra-deep overburden sites based on physical constraints according to claim 7, characterized in that, In S40, a dual-source physical projection condition adversarial training framework is adopted, a simulated seismic field condition discriminator and a real seismic event condition discriminator are set, and the physical constraint limit is incorporated into the training objective through augmented Lagrange multipliers.

9. The method for analyzing the seismic hazard of ultra-deep overburden sites based on physical constraints according to claim 8, characterized in that, In S50, the dual-source applicability domain and sample credibility of the scenario to be expanded are determined according to the following method: the source parameters, propagation path parameters, target field coordinates, overburden thickness and velocity structure, predicted bedrock input intensity, singular values ​​and group delay of the predicted site transfer matrix, and predicted ground motion intensity indices are embedded as conditions. The simulation condition domain and the real event condition domain are established by using the conditional embedding of each high-fidelity strong earthquake scenario in the simulation training set and the conditional embedding of each real earthquake event-station record in the event training set, respectively. The physical sensitivity metric matrix is ​​constructed based on the Jacobi norm of the predicted bedrock wave field, site transfer matrix and ground motion intensity index according to each input direction. In S50, the supplementary high-fidelity scenario batches that return to S20 to perform full region-local spectral element calculation are selected according to the following method: the scenarios to be expanded in the scenario hierarchical unit that have sample confidence below the threshold, simulated conditional domain distance or real event conditional domain distance exceeding the corresponding threshold, hard physical constraints not passed, and tail-weighted effective sample size at the target intensity threshold are insufficient are formed into a candidate set; for the current high-fidelity dataset, the current variance of the previous year's exceedance rate of the design intensity threshold set is calculated using multiple cognitive uncertainty branches, and virtual high-fidelity updates are performed on each candidate scenario in the candidate batch using the posterior prediction results of the surrogate model.

10. The method for analyzing the seismic hazard of ultra-deep overburden sites based on physical constraints according to claim 9, characterized in that, In S60, the probability mass conservation multi-proposal hazard integral is performed as follows: for the h-th scenario stratification unit in the z-th potential source area, the sample is divided into initial stratification proposal sample, scenario proposal sample to be expanded and support domain supplementary proposal sample according to the actual mechanism of generating strong earthquake scenario condition vector, and the k-th actual proposal mechanism is denoted as k. For the support domain supplementary proposal mechanism established based on previous samples and acceptance results, the historical information, conditional proposal density, sample allocation number and random seed on which the proposal mechanism is based shall be frozen and versioned before sampling. Subsequent proxy model updates, high-fidelity supplementary calculation results or dangerous results shall not revert the already executed conditional proposal density.

Images