Continuous interpolation method of site characteristic period in seismic design of buildings

Abstract

The application relates to the field of building seismic design, and discloses a continuous interpolation method for site characteristic periods in building seismic design, which comprises the following steps: acquiring site terrain parameters, covering layer layering parameters and site bedrock parameters, calculating soil-rock impedance ratios, and building a site response model to solve and analyze a response spectrum; comparing and analyzing the response spectrum and a standard design response spectrum to extract discrete characteristic period anchor points, establishing an observation error matrix; fusing the terrain and bedrock parameters to establish a bedrock elevation grid, introducing the soil-rock impedance ratios to calculate spatial impedance adjustment characteristics; combining local terrain characteristics and the spatial impedance adjustment characteristics to calculate a directional penalty ratio, and converting to obtain a grid edge conduction weight; extracting a prior covariance matrix according to a grid topological structure, the conduction weight and the anchor points, combining the observation error matrix to solve a characteristic period prediction matrix and a prediction variance matrix; obtaining continuous interpolation results by relying on the matrix, and outputting a distribution diagram and an interpolation uncertainty evaluation diagram. The application improves the spatial prediction accuracy of a complex site.

Description

Technical Field

[0001] This invention relates to the field of seismic design of buildings, specifically a continuous interpolation method for site characteristic periods in seismic design of buildings. Background Technology

[0002] In seismic design of buildings, obtaining continuous periodic distributions of site characteristics is fundamental for structural dynamic response analysis. Since geological exploration points are typically sparsely distributed in actual engineering projects, spatial interpolation algorithms are commonly used to predict data for unexplored areas. However, existing conventional spatial interpolation methods largely rely on purely mathematical statistical models, failing to effectively incorporate physical constraints such as differences in deep soil and rock dynamic impedance and local slope and gradient direction of the bedrock surface during the calculation process. This results in insufficient prediction accuracy of the algorithms under complex terrain and geological conditions.

[0003] Meanwhile, existing interpolation methods typically employ isotropic spatial correlations that are only related to the geometric distance between nodes. They fail to differentiate the information transmission intensity along the contour lines of the bedrock surface versus those traversing the bedrock slope, based on the bedrock surface slope direction and soil-rock impedance differences. This leads to oversmoothing in areas with significant bedrock surface undulations or large differences in soil-rock impedance. Most algorithms can only output a single predicted value in continuous space, lacking a mechanism for calculating the predicted variance of each grid node under fixed spatial model parameters. This makes it difficult to assess the spatial uncertainty and distribution range of the input data in subsequent seismic design of the structure.

[0004] Furthermore, the continuous spatial prediction results obtained by existing algorithms are usually only used as theoretical numerical references and fail to be directly integrated with the benchmark system of current building seismic design codes. They lack a code benchmark comparison mechanism based on equivalent shear wave velocity and total overburden thickness to determine site categories. Faced with interpolated data for the entire site, engineers find it difficult to intuitively judge the degree of deviation between continuous predicted values ​​and code reference values ​​in local areas. There is a lack of a deviation verification and identification data layer linked to geographic spatial coordinates, and a lack of an interface mechanism to continuously extract characteristic periods from structural control point coordinates and generate response spectrum parameter files. This reduces the practicality of the predicted data in actual engineering parameter verification, control point selection, and structural analysis data retrieval. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a continuous interpolation method for site characteristic periods in seismic design of buildings. This method solves the problem that, in seismic design of buildings, due to the sparse distribution of exploration points, conventional spatial interpolation algorithms fail to fully incorporate physical mechanisms such as site topographic relief, overburden structure, and differences in deep soil and rock impedance, resulting in insufficient accuracy in predicting the distribution of site characteristic periods in continuous space and a lack of uncertainty quantification.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a continuous interpolation method for site characteristic periods in seismic design of buildings, comprising the following steps: The topographic parameters of the site are obtained, and the soil-rock impedance ratio is calculated by obtaining and using the overburden layer parameters and the bedrock parameters of the site. Then, the site response model is built based on the overburden layer parameters to solve and analyze the response spectrum. By comparing and analyzing the response spectrum with the obtained standard design response spectrum, discrete characteristic periodic anchor points are extracted, and an observation error matrix is ​​established. A bedrock elevation grid is established by integrating the site topographic parameters and the site bedrock parameters, and the soil-rock impedance ratio is introduced to calculate the spatial impedance adjustment characteristics; The grid topology of the bedrock elevation grid is established, the direction penalty ratio is calculated by combining the local topographic features and spatial impedance adjustment features of the bedrock elevation grid, and the grid edge conduction weight is obtained by inverse monotonic mapping transformation. The prior covariance matrix is ​​extracted based on the grid topology, grid edge conduction weights and discrete feature period anchor points, and then the observation error matrix is ​​introduced to jointly solve the feature period prediction matrix and prediction variance matrix. Based on the characteristic period prediction matrix and the prediction variance matrix, the interpolation results of the characteristic period of continuous sites are obtained, thereby outputting the characteristic period distribution map of continuous sites and the interpolation uncertainty assessment map.

[0007] Furthermore, the process of obtaining and utilizing the overburden layer stratification parameters and bedrock parameters to calculate the soil-rock impedance ratio, and then using the overburden layer stratification parameters to build a site response model to solve and analyze the response spectrum, includes: accumulating the thickness of each soil layer based on the overburden layer stratification parameters to obtain the total overburden thickness, dividing the total overburden thickness by the sum of the thickness of each soil layer and the shear wave velocity ratio corresponding to each soil layer to obtain the equivalent shear wave velocity of the corresponding exploration point; and simultaneously dividing the sum of the areal densities of each soil layer by the total overburden thickness to obtain the equivalent density of the overburden at the corresponding exploration point. Furthermore, the product of bedrock density and bedrock shear wave velocity in the bedrock parameters of the site is calculated and divided by the product of the equivalent density and equivalent shear wave velocity of the overburden to obtain the soil-rock impedance ratio. Based on the overburden layer parameters, a site response model for each exploration point is constructed. Under a unified bottom boundary condition and a unified ground motion amplitude modulation rule, multiple bedrock acceleration time histories are input for calculation to obtain the corresponding analytical response spectrum. This completes the equivalent transformation of physical and mechanical parameters, thereby obtaining the impedance index reflecting the dynamic characteristics of the deep site profile, allowing one-dimensional borehole data to be directly used as boundary conditions for evaluating the dynamic response characteristics of the site.

[0008] Furthermore, the process of comparing and analyzing the response spectrum with the obtained standard design response spectrum to extract discrete characteristic period anchor points and establish an observation error matrix includes: combining a preset weighting function with respect to the natural vibration period of the structure, performing a weighted square integral on the difference between the analyzed response spectrum and the standard design response spectrum to construct a spectral error function, solving for the minimum value of the spectral error function within a preset search interval to find the corresponding undetermined characteristic period value, and calculating the arithmetic mean of multiple bedrock ground motion conditions as discrete characteristic period anchor points of the exploration points; The sample variance of the numerical values ​​of the undetermined characteristic period is calculated to obtain the observation error variance corresponding to the discrete characteristic period anchor point. The sample variance is numerically compared with the preset positive minimum error variance, and the maximum value of the two is selected as the anchor point observation error variance of the corresponding exploration point. An observation error matrix is ​​constructed with the anchor point observation error variance of each exploration point as the main diagonal element, and the remaining elements are set to zero. Based on this, a numerical optimization mechanism of spectral error is introduced to determine the interpolation anchor point. At the same time, the observation error matrix is ​​directly constructed by combining the sample fluctuation variance, thereby quantifying the prior uncertainty of local exploration point data in the early stage of calculation.

[0009] Furthermore, the process of establishing a bedrock elevation grid by integrating site topographic parameters and bedrock parameters includes: constructing a bedrock elevation grid by integrating digital elevation data of the site surface and bedrock elevation data contained in the site bedrock parameters; calculating the two-dimensional spatial gradient of the elevation values ​​of the bedrock elevation grid nodes and extracting the modulus of the two-dimensional spatial gradient as the local bedrock slope; when the modulus is greater than or equal to the lower limit of the effective slope gradient, dividing the negative two-dimensional spatial gradient by the modulus of the two-dimensional spatial gradient to extract the unit direction vector as the slope direction vector; otherwise, setting the slope direction vector and the amplitude-limited normalized slope to zero or 0 simultaneously. Subsequently, the maximum value of the preset statistical quantile of non-zero slope values ​​and the lower limit of the benchmark slope is taken as the benchmark slope reference threshold. The local bedrock slope is then divided by the benchmark slope reference threshold to generate a dimensionless ratio. This dimensionless ratio is then numerically compared with a preset upper limit of 1, and the smaller value is selected as the amplitude-limited normalized slope. The slope direction vector and the amplitude-limited normalized slope are used as the local topographic features of the bedrock elevation grid. This process integrates discrete surface and bedrock elevation data, extracts the local spatial gradient and slope direction of the bedrock interface, and establishes the geometric boundary for assessing the scattering and anisotropic effects of topographic undulations on seismic wave propagation.

[0010] Furthermore, the process of calculating the spatial impedance adjustment characteristics of the soil-rock impedance ratio includes: calculating an exponentially decaying spatial reliability coefficient based on the spatial distance between the bedrock elevation grid node and the nearest soil-rock impedance ratio acquisition point; extending the natural logarithm of the soil-rock impedance ratio to the bedrock elevation grid node through spatial grid interpolation, and obtaining the soil-rock impedance ratio of the grid node through inverse exponential transformation; when the soil-rock impedance ratio of the grid node is less than or equal to 1, the impedance modulation factor is set to zero; when the soil-rock impedance ratio of the grid node is greater than 1, the natural logarithm of the soil-rock impedance ratio of the grid node is taken to extract the amplitude-limited logarithmic impedance variable; and calculating the impedance smoothing adjustment constant using a preset reference impedance ratio and a preset reference modulation level. The impedance modulation factor generated at the corresponding bedrock elevation grid node is obtained by dividing the amplitude-limited logarithmic impedance variable by the sum of the impedance smoothing adjustment constant and the amplitude-limited logarithmic impedance variable; and the spatial reliability coefficient and the impedance modulation factor are used as spatial impedance adjustment characteristics. By using the distance attenuation function and logarithmic transformation, the discrete soil-rock impedance ratio is transformed into a spatially continuous adjustment factor, thereby forming a numerical constraint on the seismic amplification effect of deep impedance differences in the calculation.

[0011] Furthermore, the grid topology of the bedrock elevation grid is established, and the direction penalty ratio is calculated by combining the local terrain features and spatial impedance adjustment features of the bedrock elevation grid. The process of obtaining the grid edge conduction weight through reverse monotonic mapping transformation includes: for each grid edge in the grid topology, at the midpoint of the grid edge, the amplitude-limited normalized slope and slope direction vector contained in the local terrain features are extracted by interpolation calculation, as well as the impedance modulation factor and spatial confidence coefficient contained in the spatial impedance adjustment features, and the square of the inner product of the grid edge direction vector and the slope direction vector at the midpoint of the corresponding grid edge is calculated. Subsequently, the dynamic penalty increment is obtained by multiplying the penalty scaling factor, spatial reliability factor, amplitude-limiting normalized slope, impedance modulation factor, and the square of the inner product. This, combined with the fundamental penalty constant, yields the dimensionless directional penalty ratio. Then, the length of the common boundary of the control area corresponding to the grid edge is divided by the product of the node center spacing and the directional penalty ratio to convert it into grid edge connectivity weights. By using the inner product of the grid edge direction and the actual slope direction, the anisotropic propagation characteristics of seismic waves caused by topographic undulations can be directly transformed into connectivity weights within the grid structure.

[0012] Furthermore, the process of extracting the prior covariance matrix based on the grid topology, grid edge conduction weights, and discrete feature periodic anchor points includes: extracting the grid edge conduction weights corresponding to each grid edge based on the grid topology; for each bedrock elevation grid node in the grid topology, summing the conduction weights of each grid edge connected to the bedrock elevation grid node as the main diagonal element of the graph Laplace matrix, using the inverse of the grid edge conduction weights between connected bedrock elevation grid nodes as the off-diagonal element, and setting corresponding elements that are not directly connected to each other to zero, thereby generating the graph Laplace matrix; Construct an area diagonal matrix with the control area of ​​the bedrock elevation grid nodes as the main diagonal element, multiply it by the square of the preset spatial scale control parameter, add it to the graph Laplace matrix, and then multiply it by the preset accuracy scale parameter to synthesize a sparse accuracy matrix. Simultaneously calculate the bilinear interpolation weights of the adjacent bedrock elevation grid nodes around the projected position of the discrete characteristic periodic anchor point to the target discrete characteristic periodic anchor point; The bilinear interpolation weights are then written into the corresponding rows of the observation mapping interpolation matrix for the corresponding exploration points and the corresponding columns of the adjacent bedrock elevation grid nodes, with all other elements set to zero, to construct the observation mapping interpolation matrix. Finally, the observation mapping interpolation matrix, the inverse of the sparse precision matrix, and the transpose of the observation mapping interpolation matrix are multiplied consecutively to extract the prior covariance matrix for the corresponding location. By combining the graph Laplacian matrix and the area diagonal matrix, a sparse precision matrix reflecting the physical connectivity of the site can be constructed, thereby reducing the matrix solution dimensionality of spatial covariance calculations.

[0013] Furthermore, the process of jointly solving the characteristic period prediction matrix by introducing the observation error matrix includes: adding the prior covariance matrix and the observation error matrix to obtain the total covariance matrix; and calculating the weighted observation sum and the normalized weight sum based on the inverse matrix of the total covariance matrix, the anchor observation vector composed of discrete characteristic period anchor points, and the all-1 vector, and dividing the two to estimate the mean of the global constant; then using the residual vector of the anchor observation vector and the mean of the global constant, the inverse matrix of the total covariance matrix, and the transpose matrix of the observation mapping interpolation matrix, calculating the bias load vector projected onto the bedrock elevation grid, and constructing a sparse linear equation system corresponding to the bedrock elevation grid based on the synthesized sparse accuracy matrix and the bias load vector. Finally, by solving the sparse linear equations, the column vector of numerical fluctuations for the bedrock elevation grid nodes is obtained. This column vector is then added to the mean of the global constants and reconstructed according to the scale dimension of the bedrock elevation grid to generate a characteristic period prediction matrix. Based on this derivation, solving the constructed sparse linear equations allows for the fusion of sparse observation anchors and the prior structure of the entire grid, ultimately calculating the characteristic period distribution values ​​for continuous intervals.

[0014] Furthermore, the process of solving the prediction variance matrix includes: multiplying the transpose of the observation mapping interpolation matrix, the inverse of the observation error matrix, and the observation mapping interpolation matrix together to obtain the accuracy increment matrix; adding the accuracy increment matrix to the synthesized sparse accuracy matrix to obtain the posterior accuracy matrix; then, based on the size of the posterior accuracy matrix, using a selective inversion algorithm or random diagonal estimation to solve the posterior accuracy matrix, extracting the marginal posterior variance at the main diagonal position; and finally, reconstructing and arranging the obtained marginal posterior variance according to the two-dimensional row and column coordinates to generate the prediction variance matrix. The above matrix solving process extracts the marginal variance term, which reflects the spatial interpolation error caused by the uneven distribution of exploration points and the randomness of soil layer parameters, thereby obtaining the distribution variance of the prediction result.

[0015] Furthermore, after obtaining the continuous site feature period interpolation results data, it also includes outputting the deviation verification mark surface and extracting the control point feature period, specifically including: based on the overburden layer layering parameters, performing spatial interpolation on the equivalent shear wave velocity and the total thickness of the overburden layer at each exploration point to generate target distribution data, and classifying the bedrock elevation grid nodes corresponding to the target distribution data according to the preset site classification standard to determine the site category of each node, and obtaining the corresponding standard feature period value based on the site category to generate a standard reference matrix; Next, the predicted feature period value at a single grid node extracted from the feature period prediction matrix is ​​subtracted from the standard feature period value at the same grid node extracted from the standard reference matrix to obtain the predicted fluctuation deviation term. This term is then divided by the standard feature period value to obtain the relative deviation value. Finally, based on the relationship between the relative deviation value and the preset deviation verification threshold, as well as the relationship between the predicted feature period value and the standard feature period value, the corresponding grid node is assigned a feature period difference verification prompt, a data reliability prompt, or a compliance judgment value indicating in-band agreement. The node coordinate set containing the compliance judgment value is then bound to the three-dimensional spatial coordinates corresponding to the bedrock elevation grid, and the deviation verification marker surface is output. Simultaneously, the extraction of control point characteristic periods includes: determining the corresponding rectangular grid cells when the structural engineering control points are located within the site grid coverage area; constructing a control point observation mapping row vector based on the normalized coordinates of the structural engineering control points within the rectangular grid cells; and then performing matrix multiplication between the control point observation mapping row vector and the characteristic period prediction column vector of the bedrock elevation grid to calculate the control point characteristic period. This logic compares the spatial interpolation results with the benchmark values ​​in the seismic design code, and, in conjunction with deviation status indicators and control point coordinate mapping, outputs the actual characteristic period and compliance judgment information at the structural engineering control points.

[0016] This invention provides a continuous interpolation method for the site characteristic period in seismic design of buildings. It has the following beneficial effects: 1. This invention calculates the soil-rock impedance ratio by obtaining the overburden layer layer parameters and bedrock parameters of the site, combines this with the site response model to calculate and analyze the response spectrum, and integrates the site surface digital elevation and bedrock elevation data to extract local topographic features. This process incorporates the deep soil-rock dynamic impedance differences and local topographic slope as physical constraints into spatial interpolation, overcoming the shortcomings of conventional pure mathematical interpolation algorithms that are detached from the site's physical mechanisms, and improving the accuracy of predicting the periodic spatial distribution of site features under complex terrain and strata conditions.

[0017] 2. This invention derives the directional penalty ratio and grid edge conduction weights by calculating the local topographic features and spatial impedance adjustment characteristics of a bedrock elevation grid. It then extracts the prior covariance matrix and combines it with the observation error matrix, whose diagonal elements are predominantly anchor point observation error variance, to solve for the characteristic period prediction matrix and prediction variance matrix. This process transforms the anisotropic attenuation characteristics of seismic waves propagating in irregular terrain into network weights in numerical calculations. While analyzing the continuous spatial characteristic period prediction values, it also quantifies the variance of the interpolation error distribution of each grid node, providing a clear and reliable data basis for structural seismic design.

[0018] 3. This invention extracts discrete characteristic period anchor points by comparing and analyzing the response spectrum with the standard design response spectrum. After obtaining the interpolation results of continuous site characteristic periods, it determines the site category of the bedrock elevation grid nodes based on the spatial interpolation results of equivalent shear wave velocity and total overburden thickness to generate a standard reference matrix. Then, it compares the predicted value of the characteristic period extracted from the characteristic period prediction matrix with the standard reference matrix. Based on the deviation relationship of the comparison, it assigns corresponding compliance judgment state values ​​and binds them to three-dimensional terrain coordinates. Finally, it outputs the deviation verification mark surface and the characteristic period of control points, so that the continuous spatial prediction data generated by the algorithm is directly aligned with the current seismic design code for buildings. This makes it easier for engineering designers to intuitively judge the reliability and compliance of local area data and directly extract the design parameters of specific engineering control points, thus enhancing practicality. Attached Figure Description

[0019] Figure 1 This is a system architecture diagram of the present invention; Figure 2 This is a flowchart of the method of the present invention; Figure 3 The flowchart illustrates the workflow for obtaining field parameters, generating and analyzing response spectra, and extracting discrete characteristic periodic anchor points and observation error matrices for this invention. Figure 4 This is a flowchart illustrating the process of calculating the spatial physical geometry modulation coefficients and converting them into grid edge conduction weights according to the present invention. Figure 5 This is a flowchart illustrating the workflow of constructing a sparse precision matrix and solving for the prediction matrix in this invention. Figure 6 This is a flowchart illustrating the workflow of constructing a verification identification surface and outputting dynamic analysis interface data according to the present invention. Figure 7 This is a continuous site feature periodic distribution map according to an embodiment of the present invention; Figure 8 This is an interpolation uncertainty evaluation diagram according to an embodiment of the present invention. Detailed Implementation

[0020] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] See attached document Figure 1 This invention provides a continuous interpolation method for site characteristic period in seismic design of buildings. The system for executing this method includes: a data acquisition module, an anchor point extraction module, a spatial attribute calculation module, a matrix construction module, an interpolation solution module, and an interface output module.

[0022] The data acquisition module is configured to extract the overburden layer parameters and site parameters from the exploration points in the field area, calculate the equivalent dynamic parameters of the overburden and the soil-rock impedance ratio, and then construct the site response model for each exploration point. It also inputs the acquired bedrock acceleration time histories into the site response model to perform one-dimensional site seismic response analysis to obtain the analysis response spectrum.

[0023] The anchor point extraction module is configured to process the analyzed response spectrum and the standard design response spectrum based on the normalized spectral shape matching algorithm, extract discrete characteristic periodic anchor points, and determine the anchor point observation error variance based on the sample variance to establish the observation error matrix.

[0024] The spatial attribute calculation module is configured to reconstruct the bedrock elevation grid based on the basic data source and calculate spatial attribute parameters including the slope direction vector, amplitude-limited normalized slope, impedance modulation factor, and spatial reliability coefficient.

[0025] The matrix construction module is configured to establish a grid topology within the field area, calculate the direction penalty ratio based on the extracted spatial attribute parameters, and convert the direction penalty ratio into grid edge conduction weights.

[0026] The interpolation module is configured to generate a graph Laplacian matrix based on the grid topology and grid edge conduction weights, and synthesize a sparse precision matrix. It estimates the global constant mean by combining discrete feature period anchors, prior covariance matrix and observation error matrix, constructs and solves a sparse linear system of equations to generate a feature period prediction matrix, and extracts the prediction variance matrix to quantify the uncertainty of the entire field.

[0027] The interface output module is configured to output a continuous site characteristic period distribution map and an interpolation uncertainty assessment map based on the characteristic period prediction matrix and the prediction variance matrix. It generates a deviation verification marker surface for design parameter deviation verification and extracts the characteristic period of control points to generate a response spectrum function file and provide the target spectrum parameters required for the ground motion time history, so as to complete the continuous interpolation application of the site characteristic period of building seismic design.

[0028] See attached document Figure 2 This invention provides a continuous interpolation method for site characteristic periods in seismic design of buildings, comprising the following steps: S10. Acquire site parameters and generate analysis response spectrum. The data acquisition module extracts the overburden layer parameters corresponding to the coordinates of all exploration points in the site. The overburden layer parameters include at least the thickness of each soil layer, the shear wave velocity of each soil layer, the density of each soil layer, the dynamic shear modulus ratio curve, and the damping ratio curve. Based on this, the equivalent dynamic parameters of the overburden, including the total thickness of the overburden, the equivalent shear wave velocity, and the equivalent density of the overburden, are calculated. At the same time, the site topographic parameters, including the digital elevation of the site surface, and the site bedrock parameters, including the bedrock elevation, bedrock shear wave velocity, and bedrock density, are acquired. The data acquisition module calculates the soil-rock impedance ratio at each exploration point based on the equivalent density of the overburden, equivalent shear wave velocity, bedrock density, and bedrock shear wave velocity. It then constructs a site response model for each exploration point based on the overburden layer parameters and acquires multiple bedrock acceleration time histories. Using uniform bottom boundary conditions and uniform seismic amplitude modulation rules, the module inputs multiple bedrock acceleration time histories into the site response model for each exploration point, obtaining the analytical response spectrum for each bedrock acceleration time history at each exploration point.

[0029] S20. Extracting Discrete Feature Period Anchor Points and Observation Error Matrix. The anchor point extraction module employs a normalized spectral shape matching algorithm to perform peak normalization processing on the extracted analytical response spectrum and the standard design response spectrum containing undetermined feature periods, respectively, to obtain the normalized analytical response spectrum and the normalized standard response spectrum. Within a preset spectral shape fitting period interval, the anchor point extraction module constructs a spectral shape error function that uses the feature period as a variable and describes the difference between the normalized analytical response spectrum and the normalized standard response spectrum. By solving for the minimum value of the constructed spectral shape error function and calculating the arithmetic mean of the undetermined feature period values ​​corresponding to each bedrock ground motion, discrete feature period anchor points are extracted, and the sample variance of the undetermined feature period values ​​is calculated to obtain the anchor point observation error variance. An observation error matrix is ​​established based on the anchor point observation error variance of each exploration point.

[0030] S30. Calculate the spatial physical geometry modulation coefficient. After anchor point extraction, the spatial attribute calculation module integrates the site surface digital elevation data (DEM) from the site topographic parameters with the bedrock elevation data from the site bedrock parameters to construct a bedrock elevation grid in a unified projected coordinate system. The module calculates the two-dimensional spatial gradient of the node elevation values ​​in the bedrock elevation grid to define the slope direction vector and local bedrock slope, and then performs amplitude-limited normalization to obtain the amplitude-limited normalized slope. The slope direction vector and the amplitude-limited normalized slope are used as local topographic features. The spatial attribute calculation module generates an impedance modulation factor based on the soil-rock impedance ratio. It calculates an exponentially decaying spatial reliability coefficient based on the spatial distance between the bedrock elevation grid node and the nearest soil-rock impedance ratio acquisition point, and uses the spatial reliability coefficient and the impedance modulation factor as spatial impedance adjustment features.

[0031] S40. Calculate the directional penalty ratio and convert it into grid edge conduction weights. The matrix construction module establishes a grid topology structure within the site, consisting of bedrock elevation grid nodes and grid edges. At the midpoint of each grid edge, interpolation is used to extract the amplitude-limited normalized slope and descent direction vectors contained in the local terrain features, as well as the impedance modulation factor and spatial reliability coefficient contained in the spatial impedance adjustment features. The matrix construction module calculates the directional penalty ratio using the grid edge direction vector and descent direction vector corresponding to the location of the grid edge. The directional penalty ratio is converted into grid edge conduction weights through inverse monotonic mapping to suppress data smoothing across geological high impedance difference boundaries in the underlying logic.

[0032] S50. Construct a sparse precision matrix and solve for the prediction matrix. The interpolation module generates a graph Laplacian matrix based on the grid topology and grid edge conduction weights, and synthesizes a sparse precision matrix by operating on the introduced area diagonal matrix. The area diagonal matrix is ​​constructed by calculating the control area of ​​the bedrock elevation grid nodes. To achieve data interaction between discrete space and continuous grid, an observation mapping interpolation matrix is ​​constructed based on the spatial coordinate relationship between discrete characteristic periodic anchor points and bedrock elevation grid nodes to map the calculation dimension, and the prior covariance matrix of the corresponding location is extracted by combining the sparse precision matrix.

[0033] After introducing the observation error matrix, the global constant mean is estimated by combining the discrete characteristic periodic anchor points and the prior covariance matrix using the generalized least squares method. Based on the global constant mean, a sparse linear equation system is constructed and solved to directly generate the characteristic periodic prediction matrix for the seismic design site of the building area covering the field grid. Subsequently, the posterior accuracy matrix is ​​constructed and the prediction variance matrix is ​​extracted using selective inversion or random diagonal estimation to quantify the uncertainty of the entire field.

[0034] S60. Construct a verification marker surface and output dynamic analysis interface data. The interface output module determines the confidence interval based on the characteristic period prediction matrix and the prediction variance matrix, obtains the continuous site characteristic period interpolation results data, and then outputs a continuous site characteristic period distribution map and the corresponding interpolation uncertainty assessment map in combination with the geospatial coordinate system. At the same time, spatial interpolation is performed on the equivalent shear wave velocity and total overburden thickness of each exploration point to generate target distribution data (i.e., the distribution of equivalent shear wave velocity and total overburden thickness of the entire field), obtains the site classification standard, classifies the bedrock elevation grid nodes corresponding to the target distribution data according to the site classification standard to determine the site category of each node, and obtains the corresponding standard characteristic period value based on the site category to generate the standard reference matrix.

[0035] The system calculates the relative deviation between the characteristic period prediction matrix and the standard reference matrix and outputs a deviation verification label surface with compliance judgment status for design parameter deviation verification. Based on the actual structural engineering control points, the system extracts the corresponding control point characteristic periods to generate response spectrum function files and provides the target spectrum parameters required for seismic motion time history selection, matching, or amplitude modulation, thus completing the continuous interpolation application of site characteristic periods for building seismic design.

[0036] See attached document Figure 3 In step S10, the process of acquiring field parameters and generating analytical response spectra includes the following sub-steps: The data acquisition module extracts the overburden layer stratification parameters corresponding to the coordinates of all exploration points within the site. These parameters include at least the thickness of each soil layer, the shear wave velocity of each soil layer, the density of each soil layer, the dynamic shear modulus ratio curve, and the damping ratio curve. This is used to calculate the equivalent dynamic parameters of the overburden layer, which include the total thickness, equivalent shear wave velocity, and equivalent density of the overburden layer. The calculation formula is as follows: ; in, Indicates the first The total number of soil layers contained in the overburden of each exploration point; Indicates the first The first of the exploration sites The thickness of each soil layer; Indicates the first The first of the exploration sites Shear wave velocity of each soil layer; Indicates the first The first of the exploration sites The density of each soil layer; Indicates the first Total thickness of overburden at each exploration point; Indicates the first Equivalent shear wave velocity at each exploration point; Indicates the first Equivalent density of overburden at each exploration point; Indicates the first The first of the exploration sites The travel time of shear waves propagating in one direction within a soil layer; Indicates the first The first of the exploration sites The surface density of each soil layer; Indicates the first The calculation involves traversing and summing all soil layers within the overburden of each exploration point.

[0037] When the seismic design standard for an engineering project specifies a different depth for calculating the equivalent shear wave velocity, the equivalent travel time method shall be used within the specified depth range, and the calculation depth and standard version shall be recorded in the model configuration data. Simultaneously, site topographic parameters including digital elevation data of the site surface, and site bedrock parameters including bedrock elevation, bedrock shear wave velocity, and bedrock density shall be acquired and used together with the overburden layer parameters as the basic data source.

[0038] The data acquisition module calculates the soil-rock impedance ratio at each exploration point based on the equivalent density of the overburden, equivalent shear wave velocity, bedrock density, and bedrock shear wave velocity. For the first... The first exploration point will be the first The product of the bedrock density and bedrock shear wave velocity at each exploration point, divided by the first... The product of the equivalent density of the overburden and the equivalent shear wave velocity at each exploration point is used to obtain the product of the first and second exploration points. Soil-rock impedance ratio at each exploration point Subsequently, the data acquisition module constructs a site response model reflecting the vertical geological profile characteristics of the exploration point based on the overburden layer stratification parameters. According to the site seismic safety assessment report or engineering seismic fortification requirements, natural or synthetic seismic waves are selected as bedrock ground motions, and multiple bedrock acceleration time histories are established.

[0039] The data acquisition module employs a one-dimensional horizontal layered shear wave propagation equivalent linear site response model. The model uses an elastic half-space boundary at the bottom and adopts uniform bottom boundary conditions and ground motion amplitude modulation (AM) rules. All exploration points use the same bedrock ground motion set, input ground motion representation, and AM target. When the input data is a bedrock outcrop acceleration time history, it is uniformly converted and input according to the outcrop wave input rules of the site response analysis program. The input type, AM target, bottom boundary type, and equivalent linear iterative convergence conditions are recorded in the model configuration data. Multiple bedrock acceleration time histories are input into the site response model of each exploration point, and one-dimensional site seismic response analysis is performed to obtain the absolute acceleration response of surface nodes, generating corresponding acceleration response spectra as the analysis response spectra for each exploration point.

[0040] See attached document Figure 3 In step S20, to reduce the interference of differences in the absolute amplitude of input ground motion, the process of extracting discrete characteristic periodic anchor points and observation error matrices includes the following sub-steps: S21. Perform peak value normalization to reduce spectral amplitude interference. The anchor point extraction module uses a normalized spectral shape matching algorithm to obtain the standard design response spectrum containing undetermined characteristic periods as defined in the specification. It iterates through each exploration point in the field area, determining the normalized peak value within a preset normalized period interval (ranging from 0.0s to 6.0s) that is consistent across all exploration points. Specifically, it extracts the analytical response spectrum corresponding to each bedrock ground motion, as well as the maximum peak acceleration of the standard design response spectrum within the normalized period interval.

[0041] Subsequently, the spectral acceleration values ​​of each analytical response spectrum at each structural natural period were divided by its corresponding maximum peak ground acceleration (MPA). Similarly, the spectral acceleration values ​​of the standard design response spectrum at each structural natural period were divided by its corresponding MPA. This process removed the absolute acceleration dimension from the response spectra, extracting dimensionless spectral shape data reflecting the spectral characteristics. Normalized analytical response spectra and normalized standard response spectra were then obtained for each bedrock ground motion. The normalized period interval covered the subsequent spectral shape fitting period interval to avoid local peaks outside the fitting period interval altering the spectral shape matching results.

[0042] S22. Construct the spectral error function and solve for the discrete characteristic period anchor points. The preset spectral fitting period interval is jointly determined by the effective definition interval of the standard design response spectrum and the effective frequency band of the bedrock ground motion used, and remains consistent across all exploration points and all input ground motions in the same field area, without changing with the natural period of a single proposed structure.

[0043] Within the preset spectral shape fitting period interval, using the characteristic period as a variable, the anchor point extraction module constructs a spectral shape error function in the form of a square integral function to describe the difference between the normalized analytical response spectrum and the normalized canonical response spectrum. A weighting function related to the structure's natural vibration period is introduced to adjust the data weighting of different period frequency bands during the fitting process. The weighting function remains uniform globally across the field and does not change with the dominant natural vibration period of any individual proposed structure. In specific implementation, the weighting function is set... In an alternative implementation, a pre-fixed piecewise constant weighting function can be used, but the boundaries and weight values ​​of each segment should be determined before anchor point extraction and kept consistent across all exploration points in the field. For the first... The first exploration point and the first For each bedrock seismic motion, the minimum value is solved to find the corresponding undetermined characteristic period value: ; in, This indicates the lower limit period of the preset spectral fitting period interval; This indicates the upper limit of the preset spectral fitting period range; Represents the natural period of the structure Integral infinitesimal element; Represents the natural period of the structure The weighting function; Indicates the first The exploration point corresponds to the first... Normalized analysis response spectrum of bedrock ground motion in the natural period of the structure The value at that location; Indicates the period of the undetermined feature; Indicates a periodicity containing undetermined features. The normalized gauge response spectrum in the natural period of the structure The value at that location; Indicates the first The exploration point corresponds to the first... Numerical results of undetermined characteristic periods extracted from bedrock ground motion; This indicates the search for undetermined characteristic periods that minimize the subsequent objective function. The computational operators; This indicates the period within the preset spectral shape fitting interval. Operators for calculating definite integrals within; This indicates that the normalized analysis reaction spectrum and the normalized standard reaction spectrum are in the periodicity... The difference in spectral acceleration at that location; This represents the integral sum of the spectral error function after introducing the weighting function.

[0044] The numerical solution process for the undetermined characteristic period is as follows: Undetermined characteristic period Within the preset search range Internal solution (where, Indicates the lower limit period of the preset search interval; This represents the upper limit of the preset search interval, which is set according to the proposed site category and design seismic group, preferably between 0.1s and 1.5s. The global minimum of the spectral error function is determined using a combination of equal-step traversal and local one-dimensional minimization, with relevant control parameters recorded in the model configuration data. Specifically, the search step size is set based on the engineering cycle accuracy (typically 0.01s to 0.05s); the termination tolerance is set based on the numerical calculation accuracy (typically ranging from 10). -6 Up to 10 -4 When the difference in spectral error corresponding to multiple candidate feature periods is less than the preset equivalent tolerance (typically 10), -6 Up to 10 -4 When selecting a candidate feature period with the smallest value, choose it as the minimum value. To ensure the uniqueness of the calculation results. When the number of input seismic ground motions is... At that time, the anchor point extraction module calculates the arithmetic mean of the values ​​of the undetermined characteristic periods corresponding to each bedrock seismic motion, and uses the calculated arithmetic mean as the first... The discrete characteristic periodic anchor points ultimately obtained from each exploration point.

[0045] S23. Calculate the variance of anchor point observation errors and construct the observation error matrix. Input the number of bedrock ground motions. Not less than 2. For the first For each exploration point, the arithmetic mean and the sample variance of the undetermined characteristic period values ​​obtained from the bedrock seismic motion calculations are calculated: ; ; ; in, This indicates the total number of bedrock ground motions input. Indicates the first The exploration point corresponds to the first... Numerical values ​​of undetermined characteristic periods for solving bedrock ground motions; Indicates the first Discrete characteristic period anchors for each exploration point (i.e., the arithmetic mean of the undetermined characteristic period values ​​of each bedrock seismic motion). Indicates the first Sample variance of undetermined characteristic period values ​​for each exploration point; This represents the preset positive minimum error variance, used to prevent the observation error matrix from becoming singular due to zero sample variance. Its value ranges from 10. -6 s 2 Up to 10 -4 s 2 ; Indicates the first Variance of anchor point observation error at each exploration point; The function operator that represents the maximum value; This indicates that the calculations for all bedrock ground motions are summed. This represents the deviation between the undetermined characteristic period value and its arithmetic mean value corresponding to a single bedrock ground motion.

[0046] The process of constructing the observation error matrix is ​​as follows: the dimension of the anchor point observation error variance is the square of time, and the minimum value of the spectral error function is output separately as the spectral matching quality index. Based on the anchor point observation error variance of each exploration point, they are arranged sequentially by number on the main diagonal, and the remaining elements are set to zero, thus constructing an invertible positive definite observation error matrix. The observation error matrix serves as the error covariance in the conditional estimation during subsequent interpolation of the surface, allowing for corresponding smoothing deviations at anchor points in continuously interpolated surfaces. This enables the transformation of local spectral matching deviations into spatial model parameters.

[0047] See attached document Figure 4 After extracting the discrete feature periodic anchor points and the observation error matrix, the process of calculating the space physical geometric modulation coefficients in step S30 specifically includes the following sub-steps: S31. Construct a bedrock elevation grid and extract spatial geometric gradients. The spatial attribute calculation module acquires the site's digital elevation data (DEM) and bedrock elevation data, unifying the data sources to the same projected coordinate system. When the bedrock elevation data represents the bedrock burial depth, the bedrock burial depth data is first interpolated to the same projected coordinates and grid position as the site's DEM. Then, the bedrock elevation is calculated by subtracting the interpolated bedrock burial depth from the site's DEM. When the bedrock elevation data already represents the bedrock elevation, grid interpolation is performed directly under the unified elevation datum. The initially formed bedrock elevation grid is processed using local polynomial smoothing with control point constraints or Gaussian smoothing with control point weights. This allows borehole exposure points and other high-confidence control points to participate in the calculation as hard constraint points or high-weight soft constraint points during the smoothing process, thereby reducing geophysical interpretation noise and avoiding the direct overwriting of local high-frequency gradients generated by control point values ​​after smoothing. Spatial grid interpolation processing uses the Kriging interpolation algorithm or the inverse distance weighted interpolation algorithm.

[0048] S32. Extract the slope direction vector and the amplitude-limited normalized slope. The spatial attribute calculation module traverses the bedrock elevation grid nodes and calculates the two-dimensional spatial gradient of the node elevation values ​​of the bedrock elevation grid nodes. The magnitude of the two-dimensional spatial gradient is used to characterize the local rate of change of bedrock elevation, and the negative two-dimensional spatial gradient direction is used to characterize the direction of the fastest descent of bedrock elevation in the horizontal projection plane, and serves as a geometric reference direction for determining whether the grid edge crosses the bedrock slope direction.

[0049] To avoid numerical divergence when calculating direction vectors in flat regions, an effective slope gradient lower bound is introduced. For regular rectangular grids, the two-dimensional spatial gradient of internal nodes is calculated using central difference, while that of boundary nodes in the computational domain is calculated using one-sided difference. The effective slope gradient lower bound is the maximum of a preset numerical calculation lower bound and the ratio of the effective resolution of the elevation data to the minimum grid spacing. The preset numerical calculation lower bound is determined based on the floating-point arithmetic precision of the computer, and its value ranges from 10. -16 Up to 10 -7 .

[0050] When the magnitude of the 2D spatial gradient is greater than or equal to the effective lower limit of the gradient, the negative 2D spatial gradient is divided by the magnitude of the 2D spatial gradient to extract the unit direction vector as the gradient direction vector. When the magnitude of the 2D spatial gradient is less than the effective lower limit of the gradient, the corresponding gradient direction vector and the amplitude-limited normalized slope are simultaneously set to zero or 0 to avoid elevation fluctuations within the measurement resolution from participating in the direction penalty calculation. The gradient direction vector is calculated using the following piecewise formula: ; in, The two-dimensional spatial projection abscissa of the bedrock elevation grid node; The two-dimensional spatial projection ordinate of the bedrock elevation grid node; Indicates at coordinates The bedrock elevation function value at the location; Represents the two-dimensional spatial gradient operator; This represents the lower limit of the effective gradient. Indicates coordinates as The slope direction vector extracted at the bedrock elevation grid nodes; The function representing the bedrock elevation in coordinates is The two-dimensional spatial gradient at the location; The function representing the bedrock elevation in coordinates is The negative two-dimensional spatial gradient at a given point is used to indicate the direction of the fastest descent in elevation; Represents the magnitude of the gradient in two-dimensional space; It represents the zero vector in the two-dimensional plane.

[0051] The local bedrock slope at each node location is acquired synchronously, and the local bedrock slope is equated to the modulus of the two-dimensional spatial gradient. The reference slope threshold is determined using the preset statistical quantiles of the slope values ​​of all valid grid nodes. In practice, the 95th quantile of non-zero slope values ​​is extracted and compared with the preset lower limit of the reference slope. The maximum value of the two is selected as the reference slope threshold. The preset lower limit of the reference slope is used to prevent division by zero errors caused by the 95th quantile of non-zero slope values ​​approaching zero in flat terrain. In practice, a value of 0.01 to 0.05 is preferred.

[0052] The modulus of the two-dimensional spatial gradient, which serves as the local bedrock slope, is divided by a reference slope threshold to generate a dimensionless ratio. This dimensionless ratio is then compared with a preset upper limit of 1, and the smaller value is selected as the amplitude-limited normalized slope. When the magnitude of the two-dimensional spatial gradient of all grid nodes is less than the lower limit of the effective slope gradient, the total amplitude-limited normalized slope is... Set it directly to 0, and use the slope direction vector and the amplitude-limited normalized slope as local terrain features.

[0053] S33. Calculate the bounded impedance modulation factor. The spatial attribute calculation module acquires the soil-rock impedance ratio at each exploration point within the field area calculated by the data acquisition module. , for the natural logarithm Spatial gridding interpolation is performed to extend to the bedrock elevation grid nodes, and the soil-rock impedance ratio of the grid nodes that are always positive is obtained through inverse exponential transform. This implementation uses unilateral high impedance ratio modulation. When the soil-rock impedance ratio of the grid nodes... When the impedance modulation factor is zero, no additional edge penalty is added due to the high impedance bedrock; when the soil-rock impedance ratio of the grid nodes is zero... The impedance modulation factor is increased based on the degree of increase in the soil-rock impedance ratio of the grid nodes relative to 1. Unilateral high impedance ratio modulation is used to characterize the additional spatial modulation when the bedrock shear wave impedance is higher than the overburden shear wave impedance. As an optional implementation, when bidirectional impedance differences need to be considered, the absolute value of the natural logarithm of the soil-rock impedance ratio of the grid nodes is used to extract the amplitude-limited logarithmic impedance variable.

[0054] To ensure the numerical stability of spatial interpolation calculations under conditions of significant differences in geological impedance and to prevent the interpolation weight calculation results from becoming invalid due to local anomalies, a logarithmic fractional function is used to perform nonlinear mapping compression on the soil-rock impedance ratio to generate a bounded impedance modulation factor. When the soil-rock impedance ratio of a grid node is greater than 1, the natural logarithm of the soil-rock impedance ratio of the grid node is taken to extract a limited logarithmic impedance variable. Subsequently, combined with the impedance smoothing adjustment constant determined according to the preset reference impedance ratio and reference modulation level, the limited logarithmic impedance variable is divided by the sum of the impedance smoothing adjustment constant and the limited logarithmic impedance variable to obtain the impedance modulation factor generated at the corresponding bedrock elevation grid node. The formula set for generating the impedance modulation factor is as follows: ; ; ; in, Indicates at coordinates Soil-rock impedance ratio at grid nodes; This indicates a preset reference impedance ratio with a value greater than 1, and its value ranges from 2.0 to 5.0. This represents a preset reference modulation level with a value between 0 and 1, preferably 0.5; represents the impedance smoothing adjustment constant used to adjust the smooth transition characteristics of the function, so that when the soil-rock impedance ratio at the grid node reaches the preset reference impedance ratio, the generated impedance modulation factor is the reference modulation level; 1 represents the lower limit constant to ensure that the logarithmic operation result is non-negative. The function operator that represents the maximum value; and Both represent the natural logarithm operator; Indicates coordinates as The limiting logarithmic impedance variable at the point; This indicates a limit on the soil-rock impedance ratio to ensure that the argument involved in the logarithmic operation remains greater than or equal to 1. Indicates coordinates as The impedance modulation factor generated at the bedrock elevation grid nodes.

[0055] S34. Extract the spatial reliability coefficient exhibiting distance decay characteristics. Since the soil-rock impedance ratio is obtained from discrete exploration boreholes, the reliability of the inferred data decreases in unknown areas far from the measured data points. To limit the data fitting effect of the interpolation algorithm in areas lacking data support, the spatial attribute calculation module introduces a decay radius to construct a Gaussian distance decay function. The spatial reliability coefficient exhibiting exponential decay is calculated based on the spatial distance from the bedrock elevation grid node to the nearest soil-rock impedance ratio acquisition point. The variogram fitting process of spatial semi-variogram analysis is a well-known technique in this field and will not be elaborated here. The formula for calculating the spatial reliability coefficient is as follows: ; in, Indicates coordinates as The spatial distance between the bedrock elevation grid node and the nearest point where the soil-rock impedance ratio is obtained; This represents the preset confidence decay radius, the value of which is determined by spatial semivariogram analysis of the distribution characteristics of geological parameters in the field area; Indicates coordinates as The spatial reliability coefficient calculated at the bedrock elevation grid nodes; The exponentiation operator represents the exponentiation operation with the natural constant as the base; This represents the ratio of spatial distance to the credibility decay radius. The term representing the square of the ratio of spatial distance to the credibility decay radius; This represents the negative exponential decay power of the spatial credibility coefficient.

[0056] The spatial reliability coefficient smoothly approaches zero as the distance increases, thus enabling the overall algorithm to smoothly degenerate into an isotropic conventional distance-weighted interpolation form in areas lacking measured geological properties. The spatial reliability coefficient and impedance modulation factor are used as spatial impedance adjustment features.

[0057] See attached document Figure 4 After completing the calculation of the space physics geometry modulation coefficients, the process of calculating the direction penalty ratio and converting it into mesh edge conduction weights in step S40 specifically includes the following sub-steps: S41. Construct the grid topology and extract grid edge attribute parameters. This embodiment uses a regular rectangular grid consistent with the bedrock elevation grid. The grid spacing in the horizontal and vertical directions is set as the first directional spacing and the second directional spacing, respectively. The matrix construction module uses a four-neighborhood connection structure that only connects adjacent grid nodes in the horizontal or vertical directions, discretizing the continuous field space into a grid topology composed of grid nodes and grid edges.

[0058] Each grid edge in the grid topology is traversed to obtain its length, and the unit direction vector from the starting point to the ending point is calculated as the grid edge direction vector. To obtain the physical and geometric properties of the grid edge locations, attribute interpolation is performed at the midpoint of each grid edge. Interpolation is used to extract the amplitude-limited normalized slope and slope direction vector contained in the local terrain features, as well as the impedance modulation factor and spatial confidence coefficient contained in the spatial impedance adjustment features. Specifically, for the amplitude-limited normalized slope, impedance modulation factor, and spatial confidence coefficient, the arithmetic mean of the corresponding attribute parameters at the two end nodes or the midpoint bilinear interpolation value is used as the corresponding parameters.

[0059] For the descent direction, midpoint interpolation is first performed on the two-dimensional gradient components of the two end nodes to obtain the midpoint gradient of the grid edge. Then, the magnitude of the midpoint gradient is calculated: when the magnitude of the midpoint gradient is greater than or equal to the effective lower bound of the descent gradient, the midpoint descent direction vector is equal to the negative vector of the corresponding gradient divided by the gradient magnitude; when the magnitude of the midpoint gradient is less than the effective lower bound of the descent gradient, the midpoint descent direction vector is zero. The effective lower bound of the descent gradient prevents overflow errors (denominator equal to zero) in numerical calculations.

[0060] S42. Calculate the direction penalty ratio. After the grid edge attribute parameters are extracted, the matrix construction module calculates the dot product between the grid edge direction vector and the slope direction vector at the corresponding midpoint of the grid edge. The square of the dot product reflects the consistency between the grid edge direction vector and the slope direction vector: when the grid edge direction is parallel to the slope direction of the bedrock terrain, the square of the dot product tends to its maximum value; when the grid edge direction extends along the contour lines of the bedrock terrain, the square of the dot product tends to 0.

[0061] By multiplying the inner product square term with the amplitude-limited normalized slope, impedance modulation factor, and spatial reliability coefficient, and adding a basic penalty constant, a dimensionless directional penalty ratio is calculated. This nonlinearly amplifies the interpolation resistance along the slope direction in areas with steep slopes and large impedance differences. The basic penalty constant is set to 1 to ensure that the calculated directional penalty ratio is greater than or equal to 1. The penalty scaling factor is used to adjust the penalty intensity, with a value ranging from 0 to 100. Its value is positively correlated with the degree of geological variation in the site area. It is specifically calibrated using a combination of grid search and cross-validation, and jointly estimated with other model parameters. The formula for calculating the directional penalty ratio is as follows: ; in, Indicates the first The direction penalty ratio calculated for each mesh edge; 1 represents the basic penalty constant; Indicates the penalty scaling factor; Indicates the first Spatial reliability coefficient at the midpoint of a grid edge; Indicates the first Limiting normalized slope at the midpoint of each grid edge; Indicates the first Impedance modulation factor at the midpoint of a grid edge; Indicates the first The direction vector of the grid edge; Indicates the first The slope direction vector at the midpoint of a grid edge; This represents the square of the dot product between the grid edge direction vector and the slope direction vector; This indicates a dynamic penalty increment.

[0062] Since the spatial reliability coefficient smoothly approaches zero with increasing distance from the measured impedance ratio points, when it approaches zero, the direction penalty ratio in the formula approaches the base value of 1. This restores the grid edge conduction weights of the corresponding region to the baseline isotropic graph conduction weights without impedance and direction modulation. The degenerate result is an isotropic graph Laplace spatial condition estimate, rather than an inverse distance weight interpolation. To calibrate the attenuation range of the spatial reliability coefficient, when the number of measured soil-rock impedance ratio points is not less than the preset minimum number, the reliability attenuation radius is taken as the range obtained by fitting the experimental semivariogram of the soil-rock impedance ratio, where the preset minimum number is set to 30 based on statistical experience. When the number of measured points is insufficient for stable fitting or the semivariogram fitting does not converge, twice the median nearest neighbor distance of the measured soil-rock impedance ratio points is taken as the reliability attenuation radius. When there are fewer than 3 measured soil-rock impedance ratio points, spatial impedance modulation is not performed, all spatial reliability coefficients are set to 0, and an insufficient impedance data indicator is output.

[0063] S43. The direction penalty ratio is converted into grid edge conduction weights through inverse monotonic mapping. Based on the constructed four-neighborhood regular rectangular grid, the matrix construction module, based on the finite volume conduction form, converts the calculated direction penalty ratio into grid edge conduction weights through inverse monotonic mapping. The specific conversion rule is set as follows: the grid edge conduction weight is equal to the common boundary length of the corresponding control region divided by the product of the node center spacing and the direction penalty ratio.

[0064] The inverse monotonic mapping operation places the directional penalty ratio in the denominator of the calculation. This results in an increased directional penalty ratio and a corresponding decrease in the conduction weight of the grid edge when the grid edge crosses a geological boundary with a large slope and high impedance difference. Conversely, when the grid edge extends parallel to the geological boundary, the directional penalty ratio approaches the basic penalty constant, and the conduction weight of the corresponding grid edge maintains its initial proportion. This inverse monotonic mapping operation suppresses data smoothing across geologically high impedance difference boundaries and preserves the numerical distribution differences of site characteristic periods in geologically volatile regions.

[0065] See attached document Figure 5 In step S50, the process of constructing the sparse precision matrix and solving the prediction matrix specifically includes the following sub-steps: S51. Generate the graph Laplacian matrix, area diagonal matrix, and observation mapping interpolation matrix. Based on the regular rectangular mesh and mesh topology constructed in step S41, the interpolation module generates the graph Laplacian matrix. The specific calculation logic is as follows: For horizontal mesh edges, the node center spacing is set as the first directional spacing, and the corresponding control region common boundary length is set as the second directional spacing. The formula for calculating the edge conduction weight is as follows: ; in, This represents the edge conduction weight of the horizontal grid edges; Indicates the spacing in the second direction; Indicates the spacing in the first direction; This represents the direction penalty ratio corresponding to the horizontal grid edge.

[0066] For vertical grid edges, the node center spacing is set as the second directional spacing, and the length of the common boundary of the corresponding control region is set as the first directional spacing. The formula for calculating the edge conduction weight is as follows: ; in, This represents the edge conduction weight of the vertical grid edges; Indicates the spacing in the first direction; Indicates the spacing in the second direction; This represents the directional penalty ratio corresponding to the vertical grid edge.

[0067] Iterate through all bedrock elevation grid nodes, summing the conduction weights of all grid edges connected to a single bedrock elevation grid node, and use this sum as the element on the main diagonal of the graph Laplacian matrix; use the negative of the conduction weight of the grid edge connecting two connected bedrock elevation grid nodes as the element on the off-diagonal; set the elements corresponding to bedrock elevation grid nodes that are not directly connected to each other to 0. The graph Laplacian matrix represents the data smoothness penalty mechanism between adjacent bedrock elevation grid nodes in the spatial model.

[0068] Simultaneously, the interpolation module arranges the control areas of all bedrock elevation grid nodes sequentially along the main diagonal according to their node numbers, and sets the remaining off-diagonal elements to 0 to construct an area diagonal matrix. The control area of ​​internal grid nodes is the product of the spacing in the first direction and the spacing in the second direction; the control areas of boundary nodes and corner nodes are determined by trimming the corresponding rectangular control units according to the computational domain boundary. This embodiment does not use diagonal connecting edges or irregular grid edges.

[0069] The interpolation module, based on the spatial coordinate relationship between discrete characteristic periodic anchor points and bedrock elevation grid nodes, uses a bilinear interpolation algorithm to calculate the bilinear interpolation weights of the four adjacent bedrock elevation grid nodes around the projected location of the discrete characteristic periodic anchor point to the target discrete characteristic periodic anchor point. Since the discrete characteristic periodic anchor point is extracted from the corresponding exploration point, its spatial coordinates are equivalent to the planar geometric coordinates of that exploration point. Therefore, for the first [anchor point] located within a regular rectangular grid cell... The exploration point (i.e., the first exploration point) The location of each discrete characteristic periodic anchor point), and the formula for calculating the normalized coordinates are defined as follows: ; in, Indicates the serial number of the exploration point; Indicates the first The first-dimensional normalized coordinates of each exploration point; Indicates the first The second-dimensional normalized coordinates of each exploration point; Indicates the first The first-dimensional planar geometric coordinates of each exploration point; Indicates the first The second-dimensional planar geometric coordinates of each exploration point; The smallest coordinates representing the first-dimensional boundary of a rectangular grid cell; This represents the maximum coordinate of the first-dimensional boundary of a rectangular grid cell; The minimum coordinates representing the second-dimensional boundary of a rectangular grid cell; This represents the maximum coordinates of the second-dimensional boundary of a rectangular grid cell; The value is equivalent to the spacing in the first direction. ; The value is equivalent to the spacing in the second direction. .

[0070] The bilinear interpolation weights for the four corresponding grid nodes are calculated sequentially as follows: subtract the difference between 1 and the first-dimensional normalized coordinate, then multiply by the difference between 1 and the second-dimensional normalized coordinate; multiply the first-dimensional normalized coordinate by the difference between 1 and the second-dimensional normalized coordinate; multiply the first-dimensional normalized coordinate by the difference between 1 and the second-dimensional normalized coordinate; and multiply the first-dimensional normalized coordinate by the second-dimensional normalized coordinate.

[0071] The interpolation solution module writes the four bilinear interpolation weights obtained from the calculation into the observation map interpolation matrix. The The matrix rows corresponding to each exploration point and the column positions corresponding to the four adjacent bedrock elevation grid nodes are set to 0, with all other elements set to 0. The bilinear interpolation weights are determined by the normalized coordinates within the grid cells, and inverse distance weights are not used.

[0072] S52. Model Parameter Optimization and Prior Covariance Matrix Extraction Based on Maximum Likelihood Estimation. After constructing and fixing the observation mapping interpolation matrix according to the established bilinear interpolation algorithm logic, in order to achieve statistical modeling independent of local grid resolution, the interpolation solution module obtains spatial scale control parameters and accuracy scale parameters based on the grid search results for penalized scaling coefficients within a preset candidate interval, combined with negative log-likelihood estimation. The accuracy scale parameter is used to adjust the scale of the sparse accuracy matrix as a whole; when other parameters remain unchanged, increasing the accuracy scale parameter corresponds to an overall decrease in the prior covariance. The preset candidate interval for the penalized scaling coefficient is empirically set based on the degree of geological variation in the field area, usually ranging from 0.1 to 10.0, and is divided into discrete candidate parameter points by equal logarithmic intervals and traversed.

[0073] In the traversal loop of this parameter optimization, the interpolation module dynamically pre-calculates the corresponding total covariance matrix and the estimated global constant mean for each set of candidate parameters. The formula for synthesizing the sparse precision matrix is ​​as follows: ; in, Represents the sparse precision matrix of the synthesis; Indicates the precision scale parameter; Represents the graph Laplace matrix; Indicates spatial scale control parameters; Represents an area diagonal matrix; This represents the area diagonal matrix after scaling. This represents the basic accuracy matrix. Therefore, the Graph Laplace matrix, and the matrix obtained by multiplying the area diagonal matrix by the squared terms of the spatial scale control parameters, both employ a consistent finite volume discrete scale.

[0074] Subsequently, the interpolation module combines the sparse precision matrix to extract the prior covariance matrix at the corresponding location, reflecting the statistical correlation of the spatial locations of the discrete feature periodic anchor points. The prior covariance matrix is ​​theoretically equivalent to the continuous product of the observation mapping interpolation matrix, the inverse of the sparse precision matrix, and the transpose of the observation mapping interpolation matrix.

[0075] The inverse of the sparse precision matrix is ​​not explicitly extracted during the calculation. Instead, the columns of the prior covariance matrix are obtained by solving the sparse linear equations in a single-column or block-based manner. In the traversal loop for parameter optimization, the prior covariance matrix and the total covariance matrix are calculated using the following formulas: ; in, Represent the prior covariance matrix; The inverse matrix of a sparse precision matrix; Represents the observation mapping interpolation matrix; This represents the transpose of the observation mapping interpolation matrix; Represent the total covariance matrix; This represents the observation error matrix.

[0076] After generating the total covariance matrix based on the sum of the prior covariance matrix and the observation error matrix, minimize the following formula for the negative log-likelihood function: ; in, Represents the negative log-likelihood function; Represents a constant proportionality coefficient; Represents the anchor point observation vector; This represents the estimated global constant mean; 1 represents a column vector with all elements equal to 1. Indicates the number of exploration points; The matrix representing the inverse of the total covariance matrix; The determinant of the total covariance matrix; The natural logarithm of the determinant of the total covariance matrix; Represents the global constant mean vector; Represents the transpose of the residual vector; The squared Mahalanobis distance of the residual vector is represented by the matrix. It represents twice the value of pi; The natural logarithm of pi (twice the value of pi); This represents the constant normalization term.

[0077] Subsequently, the root mean square error of the characteristic period prediction was calculated through cross-validation with one exploration point left over. The model parameter combination with the smallest error was selected to synthesize the final sparse accuracy matrix. The spatial scale control parameter, with the dimension of the reciprocal of the length, is used to control the range of data influence.

[0078] S53. Perform condition estimation and estimate the mean of global constants. The interpolation module extracts the final total covariance matrix obtained by adding the prior covariance matrix and the observation error matrix determined in the parameter optimization process. Combined with the discrete characteristic periodic anchor points, the generalized least squares method is used to estimate the mean of global constants that takes into account spatial correlation and observation uncertainty, providing an unbiased and minimum variance benchmark average state estimate for the discrete characteristic periodic anchor points.

[0079] S54. Solve the sparse linear equation system and generate the characteristic period prediction matrix. The interpolation module calculates the intermediate bias vector along the observation anchor point dimension based on the estimated global constant mean, and constructs the sparse linear equation system corresponding to the bedrock elevation grid based on the synthesized sparse accuracy matrix and the bias load vector. Then, the global constant mean is substituted into the sparse linear equation system for solving to obtain the column vector of numerical fluctuations of the bedrock elevation grid nodes relative to the global constant mean. The estimation formulas for the global constant mean, the intermediate bias vector formula, and the calculation formulas for the sparse linear equation system are as follows: ; in, This represents the estimated mean of the global constants; Represents a row vector whose elements are all 1s; The matrix representing the inverse of the total covariance matrix; 1 represents the anchor observation vector; 1 represents a column vector with all elements being 1. Indicates the intermediate bias vector; Represents the sparse precision matrix of the synthesis; The column vector representing the numerical fluctuation of the bedrock elevation grid nodes is to be determined. This represents the transpose of the observation mapping interpolation matrix; Represents the global constant mean vector; This represents the residual vector between the anchor point observation vector and the global constant mean. This represents the weighted sum of observations for the mean estimate; This represents the normalized weighted sum of the mean estimates; This represents the precision-weighted term for numerical fluctuations; This represents the bias load vector projected onto the bedrock elevation grid.

[0080] The interpolation module uses relative residuals as the convergence criterion. This applies to the offset load vector. The algorithm calculates the L2 norm of the difference between the accuracy weighted term of the numerical fluctuation and the bias load vector, and divides it by the L2 norm of the bias load vector (with an additional minimum constant lower bound to prevent division by zero). Iteration stops when the resulting ratio is less than or equal to the convergence tolerance. The convergence tolerance of the relative residual is determined based on the allowable numerical error of the prediction results and convergence sensitivity experiments, and is typically taken as 10. -6 The maximum number of iterations is set based on the total number of grid nodes and preset computation time constraints, typically ranging from 1000 to 5000. The system outputs a solution failure flag when convergence is not achieved, without directly altering the linear equation solution tolerance based on geological complexity. Those skilled in the art can use the preconditioned conjugate gradient method to solve the system of linear equations containing a symmetric positive definite sparse coefficient matrix.

[0081] After obtaining the specific values ​​of the numerical fluctuations, the interpolation module adds each value in the column vector to the global constant mean to directly generate the initial characteristic period vector of the bedrock elevation grid nodes. The initial characteristic period vector is then reconstructed and arranged according to the two-dimensional row and column coordinates of the bedrock elevation grid to generate the initial characteristic period prediction matrix. In this embodiment, anchor point error modeling, conditional mean prediction, and conditional variance calculation are uniformly performed in the linear domain of the characteristic period.

[0082] After generating the initial characteristic period prediction matrix, the system performs a physical effective range check on all grid nodes. The lower limit of the characteristic period is set based on the physical limit of the ideal rigid bedrock in the site area, typically taking a value of 0.05s; the upper limit of the effective period is set based on the maximum possible soft soil geological response limit in the region, typically taking a value of 2.0s or 3.0s. When the predicted value of a node is outside the effective range, the corresponding node is set to a model anomaly verification state, its value is forcibly limited to the corresponding effective range boundary value, and an out-of-bounds warning message is output for verifying model parameters, anchor point data, and observation errors. Finally, a characteristic period prediction matrix covering the site grid for seismic design of buildings is generated.

[0083] See attached document Figure 5 After obtaining the characteristic periodic prediction matrix of the covered field grid, step S50, which involves constructing the sparse accuracy matrix and solving for the prediction matrix, also includes full-field condition estimation based on the sparse equation solution. This includes the following sub-steps: S55. Construct the posterior accuracy matrix to characterize the overall uncertainty. The interpolation module extracts the prior sparse accuracy matrix and linearly superimposes it with the observation mapping information to synthesize the posterior accuracy matrix. When the main diagonal elements of the observation error matrix are greater than zero, the posterior accuracy matrix is ​​constructed according to the following extraction formula: ; in, Represents the posterior accuracy matrix; Represents the sparse precision matrix of the synthesis; This represents the transpose of the observation mapping interpolation matrix; The matrix representing the inverse of the observation error matrix; This represents the error-weighted observation mapping matrix; This represents the precision increment matrix provided by the observation data.

[0084] S56. Perform selective inversion and extract the prediction variance matrix. For grid sizes capable of sparse Kolesky decomposition, a selective inversion algorithm is used. Based on the Kolesky decomposition factor of the posterior precision matrix (i.e., decomposing the symmetric positive definite posterior precision matrix into the product of a lower triangular matrix and its transpose), the elements in the posterior covariance matrix corresponding to the non-zero pattern positions of the posterior precision matrix are recursively calculated from bottom to top, thereby extracting the marginal posterior variance at the main diagonal position.

[0085] For grids exceeding the preset memory size, random diagonal estimation is used. Preset minimum sample size (typically 30 to 50), maximum sample size (typically 300 to 500), and convergence threshold (typically 1% to 5%) are defined. The solution to the equations and update formulas are as follows: ; in, Represents the posterior accuracy matrix; Indicates the corresponding iteration number is The solution vector at that time; Indicates the corresponding iteration number is A random vector whose elements take values ​​of 1 or -1. Represents the total number of random samples; This represents the Hadamard element-wise multiplication operation; The matrix representing the inverse of the posterior precision matrix; This represents the estimated vector of the diagonal of the posterior variance; This represents the Hadamard dot product of a random vector and the solution vector. This represents the sum of the Hadamard product.

[0086] Once the minimum number of samples is reached, the calculation stops if the relative change between two consecutive estimates is below a preset convergence threshold; otherwise, it iterates to the maximum number of samples. Variance estimates less than zero are forcibly set to zero. The interpolation module reconstructs and arranges all acquired marginal posterior variances according to two-dimensional row and column coordinates to generate a prediction variance matrix. The prediction variance matrix reflects the statistical dispersion of the spatial interpolation results in different regions.

[0087] See attached document Figure 6 After completing the full-field condition estimation based on the sparse equation solution, the process of constructing the verification marker surface and outputting the dynamic analysis interface data in step S60 specifically includes the following sub-steps: S61. Determine the confidence interval and generate a set of characteristic period spatial distribution results. The interface output module determines the confidence interval based on the characteristic period prediction matrix and the prediction variance matrix, and extracts the prediction standard deviation by taking the square root of the corresponding elements in the prediction variance matrix. Based on the Gaussian distribution assumption, the preset confidence level quantiles are multiplied by the prediction standard deviation, and then the product is added to and subtracted from the node prediction values ​​respectively. This is used to calculate the upper and lower confidence limits for each bedrock elevation grid node. When the lower confidence limit is less than the preset lower limit of the physical effective range, the result value is forcibly set to the corresponding physical lower limit. The preset lower limit of the physical effective range is determined based on the theoretical lower bound or theoretical zero value of the characteristic period of the extremely hard rock and soil in the field.

[0088] The confidence levels are preset by the engineering design task (selectable as 90%, 95%, or 99%, corresponding to quantiles of 1.645, 1.960, and 2.576 respectively), and the relevant values ​​must be recorded synchronously in the deliverables. After completing the extreme value calculation, the upper confidence limit matrix and the lower confidence limit matrix are constructed using all the upper and lower confidence limit values.

[0089] The interface output module binds the feature periodicity prediction matrix, prediction variance matrix, upper confidence limit matrix, and lower confidence limit matrix to the site's geospatial coordinate system, converting them into a geographic raster image format. Based on the numerical range, it establishes continuous color gradient bands for pseudo-color mapping, outputting a continuous site feature periodicity distribution map covering the work area, along with the corresponding interpolation uncertainty assessment map. The output continuous site feature periodicity distribution map and interpolation uncertainty assessment map provide a basis for parameter evaluation of the site area, and the automatic association mechanism for building seismic fortification categories is isolated in the output results.

[0090] S62. Interpolation Distribution of Dynamic Parameters and Classification of Site Nodes. The interface output module obtains the equivalent dynamic parameters of the overburden at the exploration points in the site area. These parameters include the total overburden thickness, equivalent shear wave velocity, and equivalent density. A spatial distance-inverse weighted interpolation algorithm is used to perform continuous spatial mapping on a two-dimensional grid plane of the site area based on discrete exploration point data. The algorithm sets the values ​​of the grid nodes in the site area to the weighted average of the values ​​of surrounding known discrete exploration points. The weight of each known discrete exploration point is inversely proportional to a specific power of the distance between the discrete exploration point and the grid node in the site area. The distance attenuation power exponent ranges from 1.5 to 2.5, with a preferred value of 2.

[0091] When the planar Euclidean distance between a field grid node to be calculated and a discrete exploration point is less than the preset coordinate tolerance (the preset coordinate tolerance is determined based on the spatial grid resolution and the exploration point positioning accuracy), the system directly uses the measured value of the corresponding discrete exploration point without performing distance reciprocal calculation. Other field grid nodes are calculated using the inverse distance weighting formula. For field grid nodes located outside the planar convex hull of all discrete exploration points, the corresponding field grid node is marked as an extrapolation node, and its extrapolation state attribute is retained in the site classification matrix and standard reference matrix to indicate that the numerical acquisition environment of the corresponding grid node lacks surrounding boundary constraints. This allows for the acquisition of the full-field equivalent shear wave velocity and total overburden thickness distribution based on the interpolation of the equivalent dynamic parameters of the overburden layer.

[0092] Subsequently, the site classification criteria are obtained, which are determined based on the cross-delimitation table of equivalent shear wave velocity threshold and overburden thickness range in a pre-specified version of the engineering seismic design standard. The interface output module traverses the bedrock elevation grid, and based on the extracted values ​​of the full-field equivalent shear wave velocity and overburden thickness distribution at each node, it performs node classification operations on the bedrock elevation grid according to the site classification criteria and the set physical limits to determine the site category of each node, and converts the classification results into a site classification matrix with attribute labels.

[0093] S63. Generation of the Standard Reference Matrix. The interface output module obtains the design seismic group attributes of the area where the site is located. The design seismic group attributes are determined by the National Engineering Site Seismic Ground Motion Parameter Zoning Map or the Regional Seismic Safety Evaluation Report. The interface output module writes the pre-specified engineering seismic design standard name, version number, and adoption date into the model configuration data. Based on the correspondence table between site categories and design seismic groups given in the specified standard version, the classification label of each node in the site classification matrix is ​​directly mapped to a fixed standard characteristic period value. After completing the numerical mapping of all nodes, the standard characteristic period values ​​are arranged according to the spatial topology row and column dimensions of the grid to generate a standard reference matrix with a stepped spatial distribution.

[0094] S64. Relative Deviation Calculation and Deviation Verification Marker Output. The interface output module calculates the relative deviation between the characteristic period prediction matrix and the standard reference matrix. The calculation logic for the relative deviation is as follows: the predicted characteristic period value at a single grid node extracted from the characteristic period prediction matrix is ​​subtracted from the standard characteristic period value at the corresponding grid node extracted from the standard reference matrix to obtain the predicted fluctuation deviation term. Then, the predicted fluctuation deviation term is divided by the standard characteristic period value to obtain the relative deviation value.

[0095] After obtaining the relative deviation values ​​of all grid nodes, a deviation verification threshold is set to classify the deviation verification levels. This threshold is determined by the engineering-specific technical plan, model verification results, or design unit verification rules. Based on the comparison between the relative deviation values ​​and the deviation verification threshold, the interface output module constructs a deviation verification identifier surface with a compliance determination status for the bedrock elevation grid nodes. The specific status assignment logic is implemented through the following piecewise functions: ; in, Indicates the compliance determination status value; This represents the calculated relative deviation value; This indicates the preset deviation verification threshold and ; This indicates that the predicted characteristic period value is higher than the standard characteristic period value, and the impact of the characteristic period difference on the design response spectrum within the corresponding structural period range needs to be reviewed (i.e., characteristic period difference review prompt). This indicates the matching status value where the relative deviation value is within the preset verification zone (i.e., matching within the verification zone). This indicates that the predicted value of the characteristic period is lower than the standard characteristic period value, and that the data reliability and design parameter adoption conditions need to be reviewed (i.e., data reliability warning). Represents the absolute value of the relative deviation; This represents the negative of the preset deviation verification threshold.

[0096] The node coordinate set containing three state values ​​(1, 0, and -1) is bound to three-dimensional terrain coordinates to output a deviation verification marker surface with compliance judgment status. This surface is used for parameter deviation verification in site seismic design. The deviation verification marker surface highlights the warning area and simultaneously marks the node areas with extrapolated state attributes to indicate that the data in that area lacks surrounding boundary constraints. The deviation verification marker surface provides spatial location prompts for the differences between the predicted results and the standard lookup table results, but it is not the sole basis for confirming geological genetic attributes.

[0097] S65. Control Point Feature Period Extraction and Subsequent Interface File Generation. The interface output module obtains the planar geometric coordinates of the actual structural engineering control points imported from the outside. When the structural engineering control point is located within the site grid coverage area, the rectangular grid cell is determined. Based on the normalized coordinates of the structural engineering control point within the grid cell, a control point observation mapping row vector is constructed. This row vector is then multiplied by the bedrock elevation grid feature period prediction column vector to calculate the control point feature period. The control point condition prediction variance is not obtained by direct bilinear interpolation of the grid node variance, but rather by solving the following equation: ; And the calculations show that: ; in, Represents the posterior accuracy matrix; This represents the intermediate auxiliary column vector obtained through the solution process; Represents the transpose of the row vectors of the control point observation mapping; This represents the control point condition prediction variance, which includes the influence of the posterior covariance among the four nodes of the grid cell containing the control point. For cases with a large number of imported structural engineering control points, a block-based solution method with multiple right-hand terms is used to calculate the control point condition prediction variance for each structural engineering control point.

[0098] When a structural engineering control point is located on a grid node, the node value is read directly; when it is located on a grid edge, one-dimensional linear interpolation is used between the corresponding two end nodes. When the actual structural engineering control point is located outside the site grid boundary, the interface output module stops extraction and outputs an out-of-bounds flag. This embodiment does not perform characteristic period extrapolation on structural engineering control points outside the site grid coverage area.

[0099] After successfully extracting the characteristic periods of the control points, the interface output module calculates the dynamic amplification factor corresponding to the natural periods of different structures using the seismic influence coefficient attenuation function given by the seismic design code. Based on the discretization mapping result of the continuous prediction matrix, the continuous natural period interval is discretized into a data sequence with a fixed step size. In addition to the characteristic periods of the control points, the interface output module obtains the maximum seismic influence coefficient, design seismic group, damping ratio, descent attenuation parameter, and effective period range of the response spectrum given by the engineering fortification parameter file, seismic ground motion parameter zoning data, or seismic safety evaluation results. Based on the obtained maximum seismic influence coefficient, design seismic group, damping ratio, descent attenuation parameter, effective period range of the response spectrum, and characteristic periods of the control points, a two-dimensional data table consisting of the natural period and the seismic influence coefficient is generated. The two-dimensional data table is packaged and output to generate a response spectrum function file.

[0100] The interface output module outputs the control point characteristic period, the standard characteristic period value, the compliance judgment status, and the control point condition prediction variance, respectively. The response spectrum function file should include at least the control point number, control point coordinates, characteristic period used, damping ratio, period sequence, corresponding seismic influence coefficient, unit, and applicable standard version, and be output in a text file, a comma-separated file, or a data format supported by the target structure analysis software.

[0101] The characteristic period used in the response spectrum function file is determined by preset engineering parameter adoption rules, and the file records the source of the adopted characteristic period as the control point characteristic period, the standard characteristic period value, or the design adoption value determined by special review. The system maintains the independence of the control point characteristic period and the standard characteristic period value. At the same time, relevant parameters are extracted to construct a dynamic time history target spectrum parameter set. The dynamic time history target spectrum parameters are only used as the target response spectrum conditions for the selection, matching, or amplitude modulation of seismic time histories of control points or support areas of various structural engineering projects. When multi-support time history analysis is required, spatiotemporal related parameters are further determined in the external seismic motion generation or matching module. The response spectrum function file and the dynamic time history target spectrum parameters are directly transferred to the external finite element structural calculation system to realize cross-software data interaction and complete the continuous interpolation application of site characteristic periods for seismic design of buildings.

[0102] The following is in conjunction with the appendix Figure 7 With appendix Figure 8 This paper presents a specific application example of continuous interpolation of site characteristic periods in the seismic design of a building site with a design area of ​​1000m × 800m. The site includes a distinct bedrock concealed fault zone and a geological high impedance difference boundary caused by ancient riverbed erosion.

[0103] First, the site parameters were discretized. Within a defined 1000m × 800m rectangular area, 15 discrete engineering geological exploration points were established. The overburden thickness, initial shear wave velocity, and density of each soil layer were extracted for each exploration point, based on… Calculate the equivalent shear wave velocity of the overburden layer based on Calculate the equivalent density of the overburden. Bedrock in the site area is defined as hard rock mass with a shear wave velocity greater than 500 m / s. Calculate the ratio of bedrock impedance to the equivalent impedance of the overburden at each exploration point to obtain the soil-rock impedance ratio.

[0104] After obtaining the basic parameters, a one-dimensional equivalent linearized frequency domain analytical algorithm was used as the lower-level feature implementation of the site response model. The acceleration time histories of seven bedrock outcrops were then analyzed. The frequency domain was converted using Fast Fourier Transform, and the propagation of wave motion in the horizontally layered medium was calculated using the complex impedance transfer function. 65% of the maximum shear strain calculated at the midpoint of each soil layer was taken as the effective shear strain. The shear modulus and damping ratio of each soil layer were iteratively updated using pre-acquired dynamic shear modulus ratio and damping ratio curves. The iteration was terminated when the relative rate of change of the shear modulus and damping ratio of all soil layers in two consecutive iterations was less than a preset 5% tolerance limit. The absolute surface acceleration response was obtained using inverse Fourier transform, and an analytical response spectrum with a damping ratio of 5% was generated.

[0105] Subsequently, anchor point characteristic period extraction and observation error quantification were performed. The standard design response spectrum containing undetermined characteristic periods, as defined in the specification, was obtained. The normalized period interval was set to 0.0 s to 6.0 s. The maximum peak acceleration was extracted from both the analytical response spectrum and the standard design response spectrum, and a division operation was performed to obtain dimensionless normalized spectral shape data. The spectral shape fitting period interval was set to 0.1 s to 1.5 s. Within this interval, a spectral shape error function in the form of a square integral with the characteristic period as the variable was constructed. (Setting the weight function) ).

[0106] Within the aforementioned spectral fitting period, a grid traversal method with a step size of 0.01s is used to find the undetermined characteristic period value that minimizes the integral within a preset search interval of 0.15s to 0.90s. For each exploration point, the arithmetic mean of the undetermined characteristic period values ​​corresponding to seven bedrock ground motions is calculated as the discrete characteristic period anchor point, and the corresponding sample variance is calculated. To prevent numerical singularities, a preset positive minimum error variance of 10 is set. -4 s 2 The maximum value of this value is compared with the sample variance to form the anchor observation error variance, and a diagonal observation error matrix is ​​constructed based on this. .

[0107] In the spatial geological feature extraction stage, the digital elevation data of the site surface and the geophysical bedrock burial depth data are spatially aligned and subtracted to obtain discrete bedrock elevation data. A Kriging interpolation algorithm with hard constraints on exploration holes is used to reconstruct the bedrock elevation grid on a regular 20m×20m rectangular grid. At the grid nodes, the central difference method is used to determine the modulus and direction of the two-dimensional spatial gradient. The effective gradient lower limit is set to 10. -6 Extract the slope direction vector of the grid nodes. Extract the 95th percentile of the non-zero slope in the field area, compare it numerically with the preset lower limit of the benchmark slope, and select the maximum value of the two as the benchmark slope reference threshold. Divide the magnitude of the two-dimensional spatial gradient by the benchmark slope reference threshold to generate a dimensionless ratio, and compare the generated dimensionless ratio with the preset upper limit of the dimensionless ratio, 1, and select the smaller value of the two as the amplitude-limited normalized slope.

[0108] Furthermore, a reference impedance ratio is set. Reference modulation level The soil-rock impedance ratios of 15 exploration points were mapped to a bedrock elevation grid using natural logarithmic spatial interpolation and inverse exponential transform, based on the logarithmic fractional function. An impedance modulation factor is generated and maintained between 0 and 1. Based on the field region semivariogram fitting results, a confidence attenuation radius is set. Calculate the spatial reliability coefficient of distance decay in Gaussian form. This causes impedance modulation smoothing to fail in areas far from the exploration point.

[0109] Based on this, the conduction weights of anisotropic grid edges are calculated. A four-neighborhood grid topology connecting the bedrock elevation grid nodes is established, with a spacing of 20m in both the first and second directions. At the midpoint of each grid edge, bilinear interpolation is used to extract the amplitude-limited normalized slope, impedance modulation factor, and spatial reliability coefficient, and the square of the inner product of the grid edge direction vector and the slope direction vector is calculated. A penalty scaling factor is set. According to the formula Calculate the directional penalty ratio. Divide the length of the common boundary of the control area by the product of the node center spacing and the directional penalty ratio, and generate the grid edge conduction weights through inverse monotonic mapping. This makes the grid edge conduction weights across steep slopes of concealed fault zones approach 0.04, while the grid edge conduction weights parallel to the geological strike remain at 1.0.

[0110] Entering the full-field equation solution stage, the Laplacian matrix of the assembled graph based on the grid edge conduction weights is used. Constructing an area diagonal matrix based on grid node control area Based on the local geometric coordinates of the exploration points within the grid cells, a sparse observation mapping interpolation matrix composed of bilinear interpolation weights is constructed. For the accuracy scale parameter in the range of 0.1 to 10.0 Spatial scale control parameters Perform maximum likelihood estimation to determine the optimal parameter combination and synthesize a sparse accuracy matrix. .

[0111] Estimating the mean of global constants based on generalized least squares method Calculate the intermediate bias vector The preconditional conjugate gradient method is employed, using relative residuals as the convergence criterion, with a convergence tolerance of 10. -6 Solving large sparse linear equation systems under the given conditions Obtain the vector of numerical fluctuations across the entire field. The fluctuation is added to the mean of the global constant to generate a characteristic period prediction matrix, and a physical effective range verification is performed. The Hadamard dot product accumulation form of Radmach-distributed random vectors is used as the lower-level feature for random diagonal estimation, and this is achieved by solving... Estimate the diagonal elements of the posterior accuracy matrix to generate the prediction variance matrix.

[0112] Combined with appendix Figure 7 With appendix Figure 8 It can be seen that, Figure 7The horizontal axis represents the abscissa of the field area, and the vertical axis represents the ordinate of the field area. This figure represents the spatial continuous variation characteristics of the predicted periodic values ​​of the entire field area. The scattered points in the figure are the locations of the measured exploration points. Figure 8 The horizontal and vertical axes are also the horizontal and vertical axes of the field area. This figure represents the distribution of the prediction standard deviation at the corresponding spatial location, reflecting the evolution of uncertainty spreading to areas without observation data.

[0113] Finally, confidence intervals are generated and interface data is output. A 95% confidence level (corresponding to a quantile of 1.960) is extracted. The predicted standard deviation is extracted by taking the square root of the corresponding element in the predicted variance matrix. The quantile is then multiplied by the predicted standard deviation, and the upper and lower confidence limits are calculated by adding and subtracting the predicted values ​​from the characteristic period prediction matrix. The equivalent shear wave velocity and total overburden thickness distribution across the entire field are generated using an inverse spatial distance weighted interpolation algorithm. Site categories are classified according to seismic design codes, and a code reference matrix for the characteristic period values ​​is generated. The relative deviation between the characteristic period prediction matrix and the code reference matrix is ​​calculated for engineering verification. The characteristic periods at the coordinates of the control points of the proposed building structure are extracted, and a response spectrum function file containing the corresponding natural periods and seismic influence coefficient sequences is generated and imported into the structural finite element analysis software.

[0114] To verify the spatial interpolation performance of the sparse accuracy matrix generated by this scheme under complex geological conditions, a leave-one-out cross-validation experiment was conducted using data from the aforementioned 1000m×800m engineering site containing concealed fault zones. Three different algorithms were used to process the same data from 15 exploration points: the traditional inverse distance weighted interpolation method, the ordinary kriging interpolation method, and the sparse accuracy continuous interpolation method combining direction penalty and impedance modulation provided in this invention.

[0115] According to the error statistics of cross-validation, the root mean square error (RMSE) of the traditional inverse distance weighted interpolation method is 0.082 s, and the mean absolute error (MAE) is 0.065 s; the RMSSE of the ordinary kriging method is 0.061 s, and the MAE is 0.048 s; the RMSSE of the method of this invention is reduced to 0.034 s, and the MAE is reduced to 0.027 s. (See attached figure.) Figure 7 As can be seen from the data characteristics, traditional methods, due to their use of isotropic distance attenuation rules, produce a "data halo" smoothing effect at the high impedance difference boundaries crossing geological faults, which is inconsistent with geological origins, resulting in the flattening of predicted values ​​on both sides of the fault. This invention, however, introduces a direction penalty ratio, causing the conduction weight of the grid edges crossing steep high impedance boundaries to nonlinearly decay to near zero. This effectively cuts off the excessive smoothing of data on both sides of the fault plane in the Laplacian operator at the bottom layer of the sparse precision matrix. The final spatial surface obtained after solving the full-field linear equations at the bottom layer successfully preserves the step characteristics of the characteristic periods on both sides of the concealed fault zone.

[0116] In summary, this invention innovatively constructs a sparse precision matrix based on physical geometry modulation by integrating multi-source heterogeneous geological information from surface topography, overburden parameters, and bedrock spatial interfaces. (Appendix) Figure 7 The continuous distribution results fully demonstrate the effectiveness of the grid edge conduction weight transformation mechanism in preserving the geological boundary effect, overcoming the mathematical distortion of traditional interpolation methods under conditions of drastic local impedance changes. Simultaneously, the prediction variance extracted by jointly deriving the graph Laplacian matrix and conditional estimation provides a basis for further analysis. Figure 8 This approach provides a foundation for probabilistic quantification, clearly identifying high-risk blind spots in the seismic design of the site that urgently require supplementary geological investigation. Without altering the existing structural software interface parameter format, this technical solution utilizes a generated full-field data matrix with probability intervals and a dynamic analysis response spectrum to transform independently distributed discrete geological borehole data into spatially continuous medium mechanical load parameters suitable for structural finite element calculations. This significantly improves the reliability and precision of spatial selection of seismic fortification parameters for the engineering site.

Claims

1. A continuous interpolation method for site characteristic periods in seismic design of buildings, characterized in that, Includes the following steps: The overburden layer parameters, bedrock parameters, and topographic parameters of each exploration point in the site are obtained. Based on the overburden layer parameters, a site response model is constructed, and multiple bedrock acceleration time histories are input to calculate and analyze the response spectrum. Discrete characteristic periodic anchor points are extracted by comparing the analyzed response spectrum with the obtained standard design response spectrum, and an observation error matrix is ​​established; A bedrock elevation grid is established by integrating the site topographic parameters and the site bedrock parameters. Establish the grid topology of the bedrock elevation grid and extract the corresponding grid edge conduction weights; Based on the grid topology, the grid edge conduction weights and the discrete feature period anchor points, the prior covariance matrix is ​​extracted, and then the observation error matrix is ​​introduced to jointly solve the feature period prediction matrix and the prediction variance matrix. Based on the characteristic period prediction matrix and the prediction variance matrix, the interpolation results of the characteristic period of the continuous site are obtained, thereby outputting the characteristic period distribution map of the continuous site and the interpolation uncertainty evaluation map.

2. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 1, characterized in that, The process of comparing the analyzed response spectrum with the acquired canonical design response spectrum to extract the discrete characteristic periodic anchor points and establishing the observation error matrix includes: Peak normalization was performed on the analytical reaction spectrum and the normalized design reaction spectrum to obtain the normalized analytical reaction spectrum and the normalized normalized reaction spectrum. A spectral error function is constructed from the normalized analysis response spectrum and the normalized standard response spectrum to extract the undetermined characteristic period values. The arithmetic mean of the undetermined characteristic period values ​​corresponding to the multiple bedrock acceleration time histories is calculated and used as the discrete characteristic period anchor point for each of the exploration points. Calculate the sample variance of the undetermined characteristic period values ​​to obtain the anchor point observation error variance corresponding to the discrete characteristic period anchor points, and construct the observation error matrix with the anchor point observation error variance corresponding to each exploration point as the main diagonal element.

3. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 1, characterized in that, The process of extracting the grid edge conduction weights includes the step of calculating the soil-rock impedance ratio using the overburden layer layering parameters and the bedrock parameters of the site area, specifically including: Based on the overburden layer stratification parameters, the equivalent shear wave velocity and the equivalent density of the overburden layer at each exploration point are obtained. Then, the soil-rock impedance ratio is calculated by combining the bedrock density and bedrock shear wave velocity in the bedrock parameters of the field area.

4. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 3, characterized in that, The process of integrating the site topographic parameters and the site bedrock parameters to establish the bedrock elevation grid, and extracting the grid edge conduction weights, also includes a step of extracting local topographic features based on the bedrock elevation grid, specifically including: Obtain the site surface digital elevation contained in the site topography parameters, and merge the site surface digital elevation with the bedrock elevation data contained in the site bedrock parameters to construct the bedrock elevation grid; Calculate the two-dimensional spatial gradient of the node elevation values ​​of the bedrock elevation grid, and extract the modulus of the two-dimensional spatial gradient as the local bedrock slope; When the magnitude of the two-dimensional spatial gradient is greater than or equal to the preset effective slope gradient lower limit, the negative two-dimensional spatial gradient is divided by the magnitude of the two-dimensional spatial gradient to extract the unit direction vector as the slope direction vector, and a limit-normalized slope is generated based on the preset benchmark slope reference threshold and the local bedrock slope. When the magnitude of the two-dimensional spatial gradient is less than the preset effective gradient lower limit, the gradient direction vector is set to zero and the amplitude-limited normalized gradient is set to 0. The slope direction vector and the amplitude-limited normalized slope are used as the local terrain features of the bedrock elevation grid.

5. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 4, characterized in that, The process of extracting the grid edge conduction weights also includes the step of introducing the soil-rock impedance ratio to calculate the spatial impedance adjustment characteristics, specifically including: The spatial reliability coefficient, which decays exponentially, is calculated based on the spatial distance between the nodes of the bedrock elevation grid and the nearest point where the soil-rock impedance ratio is obtained. Spatial grid interpolation is performed on the soil-rock impedance ratio to obtain the soil-rock impedance ratio and the amplitude-limited logarithmic impedance variable of the grid nodes; The impedance smoothing adjustment constant is calculated using a preset reference impedance ratio and a preset reference modulation level. The corresponding impedance modulation factor is obtained by combining the impedance smoothing adjustment constant with the limiting logarithmic impedance variable. The spatial reliability coefficient and the impedance modulation factor are used as the spatial impedance adjustment characteristics.

6. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 5, characterized in that, The process of extracting the grid edge conduction weights also includes the step of calculating the direction penalty ratio by combining the local terrain features and the spatial impedance adjustment features, and obtaining the grid edge conduction weights through an inverse monotonic mapping transformation, specifically including: For each grid edge in the grid topology, at the midpoint of the grid edge, the local terrain features include the amplitude-limited normalized slope and the slope direction vector, as well as the impedance modulation factor and the spatial confidence coefficient included in the spatial impedance adjustment features, and the square of the inner product of the grid edge direction vector and the slope direction vector at the corresponding midpoint of the grid edge is calculated. The dynamic penalty increment is obtained by multiplying the preset penalty scaling factor, the spatial reliability factor, the amplitude limiting normalization slope, the impedance modulation factor, and the square of the inner product. The dimensionless directional penalty ratio is obtained by adding the preset basic penalty constant to the dynamic penalty increment. Then, based on the common boundary length of the control area corresponding to the grid edge, the node center spacing and the direction penalty ratio, the grid edge conduction weight is generated.

7. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 1, characterized in that, The process of extracting the prior covariance matrix based on the grid topology, the grid edge conduction weights, and the discrete feature periodic anchor points includes: Based on the grid topology, the grid edge conduction weights corresponding to each grid edge are extracted, and a graph Laplacian matrix is ​​generated for each bedrock elevation grid node using the grid edge conduction weights. Extract the control area of ​​the bedrock elevation grid nodes to construct an area diagonal matrix, and use the area diagonal matrix and the graph Laplacian matrix to synthesize a sparse precision matrix; Calculate the bilinear interpolation weights of the four adjacent bedrock elevation grid nodes around the projected position of the discrete characteristic periodic anchor point to the target discrete characteristic periodic anchor point, and write the bilinear interpolation weights into the matrix rows of each exploration point and the corresponding column positions of the four adjacent bedrock elevation grid nodes in the observation mapping interpolation matrix to construct the observation mapping interpolation matrix; Finally, the prior covariance matrix at the corresponding position is generated by the continuous product of the observation mapping interpolation matrix, the inverse matrix of the sparse precision matrix, and the transpose of the observation mapping interpolation matrix.

8. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 7, characterized in that, In the process of jointly solving the characteristic period prediction matrix and the prediction variance matrix using the observation error matrix, the process of solving the characteristic period prediction matrix includes: The total covariance matrix is ​​obtained by adding the prior covariance matrix and the observation error matrix, and the global constant mean is estimated by combining the anchor observation vector composed of the discrete feature periodic anchor points. Then, using the residual vector between the anchor point observation vector and the global constant mean, the inverse matrix of the total covariance matrix, and the transpose matrix of the observation mapping interpolation matrix, the bias load vector projected onto the bedrock elevation grid is calculated, and a sparse linear equation system corresponding to the bedrock elevation grid is constructed based on the sparse accuracy matrix and the bias load vector. Finally, the sparse linear equation system is solved to obtain the column vector of numerical fluctuation of the bedrock elevation grid nodes. The column vector of numerical fluctuation is added to the mean of the global constant to generate the characteristic period prediction matrix.

9. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 7, characterized in that, In the process of jointly solving the characteristic periodicity prediction matrix and the prediction variance matrix using the observation error matrix, the steps for solving the prediction variance matrix include: The accuracy increment matrix is ​​calculated using the observation mapping interpolation matrix, the transpose of the observation mapping interpolation matrix, and the inverse of the observation error matrix. The posterior accuracy matrix is ​​obtained by adding the accuracy increment matrix to the sparse accuracy matrix. Solve for the posterior accuracy matrix and extract the marginal posterior variances at the main diagonal positions. Reconstruct and arrange these variances according to the two-dimensional row and column coordinates to generate the prediction variance matrix.

10. The continuous interpolation method for site characteristic periods in seismic design of buildings according to claim 3, characterized in that, After obtaining the continuous site feature period interpolation result data, the method further includes the steps of outputting the deviation verification marker surface and extracting the control point feature period, specifically including: Based on the overburden layer layering parameters, the total thickness of the overburden at each exploration point is obtained. Spatial interpolation is performed between the equivalent shear wave velocity at each exploration point and the total thickness of the overburden layer to determine the site category of each bedrock elevation grid node and obtain the corresponding standard characteristic period value, thereby generating a standard reference matrix. The relative deviation value is obtained by comparing the predicted value of the feature period extracted from the feature period prediction matrix with the value of the standard feature period. Based on the relationship between the relative deviation value and the preset deviation review threshold, the corresponding node is assigned a compliance judgment status value. The compliance judgment status value is then bound to the three-dimensional terrain coordinates of the corresponding node, and a deviation review identifier surface with a compliance judgment status is output. Construct a row vector of control point observation mapping, and multiply the row vector of control point observation mapping with the column vector of feature period prediction extracted based on the feature period prediction matrix to calculate and obtain the feature period of the control points.

Images