Near-fault site multi-seismic-phase free field inversion method based on sphere / cylindrical surface diffusion characteristics

Abstract

The invention provides a near-fault site multi-seismic-phase free field inversion method based on sphere / cylindrical diffusion characteristics, and a research system in the field of near-fault inversion is perfected. On the basis of spherical diffusion characteristics of body waves and cylindrical diffusion characteristics of surface waves, a near-fault site multi-seismic-phase free field inversion framework is constructed by adopting a dynamic stiffness matrix method, and an inversion problem is converted into a multi-target optimization problem for solving an epicentral distance and a seismic source depth; on the basis, the space-time difference of the free field displacement waveform is quantified by adopting a dynamic regular path distance, and a target function is constructed by taking the minimization of the three-dimensional displacement waveform difference of two similar points as a criterion; meanwhile, constraint conditions are set in combination with propagation attenuation characteristics of seismic waves, so that a near-fault inversion multi-objective optimization model is established; and then, a rapid non-dominated multi-objective genetic algorithm is introduced to solve a Pareto optimal solution of an optimization variable, and the optimal solution is substituted into an inversion framework to complete multi-seismic-phase free field inversion of the near-fault site.

Description

Technical Field

[0001] The present invention relates to the field of seismic wave field inversion, and in particular to a multi-phase free-field inversion method for near-fault sites based on spherical / cylindrical diffusion characteristics. Background Art

[0002] Under the action of earthquakes, these structures experience non-negligible dynamic interactions with the surrounding soil, forming a typical soil-structure interaction (SSI) system. Seismic response analysis that considers the SSI effect is a prerequisite for accurately assessing the seismic risk of infrastructure, and this process relies on a reasonable seismic input method and accurate seismic excitation. Under the current boundary-equivalent load seismic input method framework, the unknown underground free field forms the seismic excitation source of the SSI system, which needs to be inverted using known surface ground motions and site conditions. The accuracy of the inverted underground free field will directly affect the rationality of subsequent structural seismic input and the reliability of response assessment, thereby affecting the safety of structural design.

[0003] Existing free-field inversion methods, whether for obliquely propagating body waves or horizontally propagating surface waves, assume a plane wave assumption: the waveline direction remains consistent, and the wavefront perpendicular to the waveline is a plane. However, this plane wave assumption only applies to far-fault regions far from the epicenter. In near-fault regions, body waves and surface waves exhibit true spherical and cylindrical diffusion characteristics, respectively. Specifically, body waves generated by fault dislocations will diffuse from the hypocenter in the form of spherical waves, with wavelines in varying directions, a spherical wavefront, and energy attenuation inversely proportional to the square of the propagation distance. Surface waves, on the other hand, will diffuse in the form of cylindrical waves, with wavelines in varying directions, a cylindrical wavefront, and energy attenuation inversely proportional to the square root of the propagation distance. Specifically, spherical and cylindrical waves attenuate in all directions, with their iso-energy lines being circles formed by the same spherical and cylindrical radii. In contrast, the energy of plane waves attenuates only along their propagation direction, with iso-energy lines being lines perpendicular to the propagation direction. The aforementioned differences in attenuation characteristics lead to different free-field responses between near-fault and far-fault regions, making it difficult for far-fault free-field inversion methods based on the plane wave assumption to reflect the free-field characteristics of near-fault regions. More importantly, near-fault regions will produce special pulse-type earthquakes due to directional effects and slip effects. Pulse-type earthquakes have significant long-period components, and the energy they carry is mainly concentrated in the low-frequency band. This low-frequency characteristic makes long-period structures in near-fault regions susceptible to "resonance effects" with pulse-type earthquakes, thereby exacerbating structural damage and destruction. SSI systems typically have long periods, and existing far-fault inversion methods have difficulty accurately reflecting the long-period pulse characteristics of near-fault pulse-type ground motions. This limitation will directly affect the seismic response results of the structure. Therefore, it is urgent to propose a free-field inversion method suitable for near-fault sites to provide accurate excitation for the seismic analysis of SSI systems located in near-fault regions, thereby providing reliable guidance for structural response assessment and seismic design. Summary of the Invention

[0004] The present invention provides a multi-phase free-field inversion method for near-fault sites based on spherical / cylindrical diffusion characteristics, which improves the research system in the field of near-fault inversion.

[0005] The technical solutions of the present invention are as follows:

[0006] Based on the spherical diffusion characteristics of body waves and the cylindrical diffusion characteristics of surface waves, a dynamic stiffness matrix method is used to construct a multi-phase free-field inversion framework for near-fault sites, and the inversion problem is transformed into a multi-objective optimization problem to solve the "epicenter distance" and "focal depth". On this basis, the dynamic regularized path distance is used to quantify the spatiotemporal differences of the free-field displacement waveforms, and the objective function is constructed based on the criterion of minimizing the difference in the three-dimensional displacement waveforms between two nearby points. At the same time, constraints are set based on the propagation attenuation characteristics of seismic waves, thereby establishing a multi-objective optimization model for near-fault inversion. Subsequently, a fast non-dominated multi-objective genetic algorithm (NSGA-II) is introduced to solve the Pareto optimal solution of the optimization variables in the model, and it is substituted into the inversion framework to complete the multi-phase free-field inversion for near-fault sites.

[0007] The steps of the multi-phase free-field inversion method for near-fault sites based on spherical / cylindrical diffusion characteristics are as follows:

[0008] Step 1: Based on the spherical diffusion characteristics of body waves and the cylindrical diffusion characteristics of surface waves, a multi-phase free-field inversion framework for near-fault-layered media sites is established in the cylindrical coordinate system.

[0009] Step 2: The unknown “epicenter distance r” in the multi-phase free-field inversion framework of the near-fault-stratified medium site HO " and "focal depth H" are used as optimization variables, and the inversion problem is transformed into a multi-objective optimization problem. A multi-objective function is constructed based on the criterion of minimizing the difference in three-dimensional displacement waveforms of two nearby points in the site. Constraints are set based on the propagation attenuation characteristics of seismic waves, and a multi-objective optimization model for near-fault inversion is established.

[0010] Step 3: Use a fast non-dominated multi-objective genetic algorithm to solve the Pareto optimal solution of the optimization variables and substitute it into the inversion framework to complete the multi-phase free-field inversion of the near-fault site.

[0011] Furthermore, the process of the first step is as follows:

[0012] Step 1.1: Assume that the earthquake source is located in the underlying bedrock half-space, establish a cylindrical coordinate system r-θ-Z with the epicenter as the origin, the line connecting the epicenter and the surface station as the radial direction, and the direction from the epicenter to the earthquake source as the vertical direction;

[0013] Step 1.2: Use the dynamic stiffness matrix method to calculate the Green's function of the spherical diffusion P wave, SV wave, and SH wave at any position (r, θ, z) in space:

[0014]

[0015] Where k is the horizontal radial wave number; and are the radial and vertical displacements of the j-th soil layer medium particle induced by the unit amplitude P wave in the wave number domain; and are the radial and vertical displacements of the j-th soil layer medium particle induced by the unit amplitude SV wave in the wave number domain; is the displacement of the particle in the jth soil layer along the circumferential direction induced by the unit amplitude SH wave in the wave number domain; are the spatial domain displacement Green's functions obtained by the inverse Hankel transform of the above wavenumber domain displacement;

[0016] Step 1.3: Take j = 1 in formula (1) and calculate the three types of spherical diffuse body waves at the surface station position (r HO ,0,0) at Green's function solution:

[0017]

[0018] In the formula, the Green function solution depends on the unknown quantity "epicenter distance r HO ” and “focal depth H”;

[0019] Step 1.4: Use the normalized inner product method to extract the body wave displacement components collected at the surface station Based on formula (5), the incident amplitudes of P wave, SV wave and SH wave are inverted.

[0020]

[0021] Step 1.5: Multiply the incident wave amplitude obtained by formula (5) by the Green's function shown in formula (1) to construct the displacement free field induced by the spherical diffuser wave at any position in the spatial domain.

[0022]

[0023] Step 1.6: Calculate the time partial differential and spatial partial differential of the body wave displacement free field, construct the corresponding velocity free field and stress free field, and complete the body wave free field inversion;

[0024] Step 1.7: Extend the body wave dynamic stiffness matrix method in cylindrical coordinates to surface waves and derive the spatial free-field displacement induced by the Rayleigh surface waves propagating on the cylindrical surface. and

[0025]

[0026] Where, and are the radial and vertical displacements of the particle in the jth soil layer induced by Rayleigh waves in the cylindrical coordinate system; and are the spatial domain displacements induced by Rayleigh waves in the plane coordinate system; i is the imaginary number, satisfying i 2 =-1; J α (kr) is the Bessel function of the first kind of order α; and are the up- and down-wave amplitudes of the non-uniform P wave in the j-th soil layer; and are the upgoing and downgoing wave amplitudes of the inhomogeneous SV wave in the medium layer, respectively; and are the vertical wave numbers of non-uniform P waves and non-uniform SV waves, respectively; is the velocity ratio of the P wave to the Rayleigh surface wave in the jth layer of medium; is the velocity ratio of the S wave to the Rayleigh surface wave in the jth layer of medium;

[0027] Step 1.8: Similarly, derive the free-field displacement induced by the Love surface wave diffused by the cylinder.

[0028]

[0029] Where, is the circumferential displacement of the particle in the j-th soil layer induced by the Love wave in the cylindrical coordinate system; is the spatial domain displacement induced by the Love wave in the plane coordinate system; and are the up- and down-wave amplitudes of the non-uniform SH wave in the j-th soil layer medium; is the vertical wave number of the non-uniform SH wave;

[0030] Step 1.9: Perform modal decoupling on the dispersive surface wave. Based on the single-mode displacement mode shape, use formula (12) to derive the vertical attenuation factor of the single-mode surface wave free-field displacement:

[0031]

[0032] Where, the subscript m represents the mth mode of the dispersive surface wave;

[0033] Step 1.10: Construct the surface station position (r HO ,0,0) and the single-mode surface wave free field at any reference point (r,0,0) on the surface; further based on the consistency criterion of the vertical motion of the two points, the attenuation factor of the single-mode surface wave free field displacement along the horizontal radial direction is derived using formula (13):

[0034]

[0035] Step 1.11: Use the normalized inner product method to extract the surface wave displacement components collected at the surface station, and convert the decoupled single-order modal displacement Multiply the corresponding horizontal radial attenuation factor and vertical attenuation factor respectively to construct the single-mode surface wave displacement free field at any position (r, θ, z) in space:

[0036]

[0037] Where, the free field displacement of a single-mode surface wave depends on the unknown quantity “the epicenter distance r HO ”;

[0038] Step 1.12: Calculate the time partial differential and spatial partial differential of the single-mode surface wave displacement free field to construct the corresponding single-mode velocity free field and stress free field;

[0039] Step 1.13: Based on the principle of surface wave dispersion mode superposition, all single-mode free fields are superimposed to form a multi-mode free field, completing the surface wave free field inversion.

[0040] Step 1.14: Superimpose the body wave free field and the surface wave free field to construct the multi-phase total free field near the fault site.

[0041] Furthermore, the second step is specifically performed as follows:

[0042] Step 2.1: The waveform difference at similar locations near the fault site is used as the target variable of the objective function, and the dynamic regularized path distance is used as a quantitative indicator to measure the waveform difference with phase difference characteristics;

[0043] Step 2.2: Construct the multi-seismic phase displacement time history of two points close to the site and The cumulative distance matrix D is composed of:

[0044]

[0045] Where p and q represent the number of rows and columns of the matrix respectively, D(p,q) is the matrix element at that position; d pq For timeline The displacement and time history of the pth discrete point in The Euclidean distance of the displacement of the qth discrete point in ;

[0046] Step 2.3: Let the regular path traverse the cumulative distance matrix D, starting from the upper left corner element position (1,1) of the matrix and ending at the last element (n,n). Calculate the final cumulative distance along the regular path:

[0047]

[0048] Where n is the total number of discrete points in the two time histories; the smaller the WPD value, the shorter the traversal path, and thus the smaller the difference in displacement waveforms between two adjacent points;

[0049] Step 2.4: Based on formula (16), a multi-objective function is constructed with the criterion of minimizing the difference in three-dimensional displacement waveforms between two adjacent points on the site:

[0050]

[0051] Where, and is the radial multi-seismic phase displacement of two nearby points on the site; and is the multi-seismic phase displacement in the tangential direction; and is the multi-seismic phase displacement in the vertical direction;

[0052] Step 2.5: Considering the attenuation characteristics of seismic waves during propagation, the displacement amplitude at the propagation lag point is Smaller than the displacement amplitude at the first point As constraints:

[0053]

[0054] Step 2.6: Keep the values of the two optimization variables greater than 0 to further narrow the feasible region:

[0055] r HO ∈R + ,H∈R + (19)

[0056] Step 2.7: Combine the objective function in formula (17) with the constraints shown in formulas (18)-(19) to establish a multi-objective optimization model for near-fault inversion:

[0057]

[0058] The beneficial effects of the present invention are as follows:

[0059] (1) The present invention provides a multi-phase free-field inversion method for near-fault sites based on spherical / cylindrical diffusion characteristics. This method fully considers the spherical diffusion characteristics of body waves and the cylindrical diffusion characteristics of surface waves, overcomes the limitations of existing far-fault inversion methods in near-fault sites, and greatly improves the inversion accuracy. Comparative analysis results show that the free-field results obtained by inversion are highly consistent with the actual borehole records, thus providing more reliable seismic excitation for the seismic response analysis of SSI systems. It has important engineering applications in revealing the true seismic response behavior of structures, as well as damage mechanisms and failure modes.

[0060] (2) The method of the present invention can achieve accurate inversion of key parameters of underground seismic motion, reveal the free field characteristics of near-fault sites, and especially accurately determine the pulse characteristics and pulse period of underground seismic motion.

[0061] The long-period pulse component in the pulse-type seismic motion is the main cause of the damage to medium- and long-period structures in the near-fault area. In particular, when the pulse period is close to the natural vibration period of the structure, it will easily induce a resonance effect. Therefore, the determination of the pulse characteristics and the calculation of the pulse period are of great significance to the seismic design of medium- and long-period structures, especially the long-period SSI system. When the near-fault surface seismic motion has pulse characteristics, the corresponding underground seismic motion pulse characteristics and pulse period values are uncertain. This uncertainty brings great challenges to the seismic response assessment and seismic design of the SSI system in the near-fault area. The method of the present invention can truly reflect the pulse characteristics of the underground seismic motion based on the surface pulse-type seismic motion and the local site conditions, thereby providing accurate excitation for the seismic analysis of the SSI system, and has important engineering application value. The formula derivation results are credible and the calculation process is efficient.

[0062] (3) The method of the present invention transforms the overdetermined problem existing in the free-field inversion of near-fault sites into an optimization problem. By introducing the theoretical framework of multi-objective optimization, the problem of solving the overdetermined set of equations that originally had no exact solution is transformed into a parameter optimization problem with clear physical meaning, and the equilibrium solution of the optimization variables under multi-objective control is obtained.

[0063] In order to solve the problem of overdetermined equations without exact solutions, the least squares method is usually used in the field of mathematics to construct a constrained least squares problem, and then the least squares problem is transformed into an optimization problem, by finding an optimal solution to approximately replace the non-existent exact solution. Drawing on this solution idea, the method of the present invention transforms the overdetermined problem of near-fault free-field inversion into the problem of finding the "epicenter distance r HO The optimization problem of the optimal solution of "and" focal depth H". This idea has clear physical significance, sufficient theoretical basis, accurate and objective results, and provides a scientific analysis method for near-fault earthquake inversion research. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 Flowchart of the method of the present invention. DETAILED DESCRIPTION

[0065] The present invention will be further described below with reference to the accompanying drawings in the embodiments of the present invention.

[0066] Attachment Figure 1 In the first step, the inversion frameworks of the body wave free field of spherical diffusion and the surface wave free field of cylindrical diffusion are constructed respectively according to the dynamic stiffness matrix method in the cylindrical coordinate system, and then the two are superimposed to construct the multi-phase total free field of the near-fault site.

[0067] According to wave theory, the frequency-domain dynamic stiffness matrix of a site in a cylindrical coordinate system is consistent with that in a plane coordinate system. Therefore, the frequency-domain dynamic stiffness matrix method based on the plane wave assumption can be extended to solve the free field induced by cylindrical waves. In addition, the free field induced by spherical waves can be converted into a superposition of cylindrical wave free fields in the frequency-wavenumber domain, thereby respectively realizing the construction of spherical wave and cylindrical wave free fields. Based on the above principles, Ba Zhenning et al. proposed Green's function solutions for spherical P waves, SV waves, and SH waves in three-dimensional layered media under the assumption of a point source. The overall idea is: first, the point source is fixed to solve the intra-layer solution and the fixed end surface reaction force, then the fixed end surface reaction force is reversed and applied to the other layers of the layered medium site to solve the reaction force solution, and finally the fixed intra-layer solution and the fixed end surface reaction force solution are superimposed to obtain the total Green's function solution.

[0068] Drawing on this concept, the method presented in this paper constructs a free-field inversion framework for body and surface waves by fixing the earthquake source in the underlying bedrock half-space of a layered medium site. The spatial domain potential functions of spherical P-waves, SV-waves, and SH-waves induced by a point earthquake in the underlying bedrock half-space can be expanded using the Hankel transform as an integral superposition of cylindrical wave potential functions in the frequency-wavenumber domain:

[0069]

[0070] Where, are the potential functions of spherical P wave, SV wave and SH wave in the spatial domain respectively; are the incident amplitudes of three types of spherical body waves; k P With k S are the wave numbers of P and S waves, and the complex compression wave velocity of the half-space field. and complex shear wave velocity Related, satisfying and R is the distance of seismic wave propagation, that is, the distance between the propagation position and the earthquake source. When converted from the spatial domain to the frequency-wavenumber domain, the seismic wave propagation distance R can be expressed as the square root of the sum of the squares of the horizontal radial propagation distance r and the vertical propagation distance z; k represents the horizontal radial wave number in the frequency-wavenumber domain; J0(kr) is the zero-order Bessel function (first kind); and are the vertical wave numbers of P waves and S waves in the frequency-wavenumber domain, respectively, satisfying:

[0071]

[0072] According to formula (21), spherical P wave, SV wave, and SH wave have the same amplitude as the unit incident wave. The potential functions of the cylindrical wave components under are:

[0073]

[0074] Based on the displacement-potential function relationship in the cylindrical coordinate system, the free-field displacement induced by the unit amplitude body wave in the bedrock half-space can be derived as follows:

[0075]

[0076] Where, and They represent the free-field displacement amplitudes of the unit-amplitude body wave in the radial, vertical, and circumferential directions in the frequency-wavenumber domain. It is worth noting that the P-wave and SV-wave are coupled with each other in the P-SV plane. Therefore, when calculating the free-field displacements induced by the P-wave and SV-wave, the potential functions of each other must be set to 0. For example, when solving the free-field displacement induced by the P-wave independently, let In contrast, when solving for the free-field displacements induced independently by SV waves, let Based on the free-field displacement obtained from the above formula, the free-field stress induced by the unit amplitude body wave can be further calculated according to the stress-displacement relationship as follows:

[0077]

[0078] In the formula, λ and μ are the Lame constants of the underlying bedrock half-space. In particular, when the upper interface of the bedrock half-space is at a depth of d from the earthquake source, R When z=-d R Substituting into formulas (24)-(25), the displacement amplitude and stress amplitude of the particle in the upper interface soil layer of the half-space can be calculated. This part is decomposed into a special solution, represented by the superscript 's':

[0079]

[0080] At this time, the amplitude of the external load on the upper interface of the half space satisfies as well as Substituting the particular solution of the displacement amplitude of the upper interface of the half-space into the local dynamic equilibrium equation of the half-space, we can obtain the homogeneous solution of the external load amplitude of the upper interface, which is also the general solution, as shown in Equations (27) and (28). The homogeneous solution is denoted by the superscript 'h'.

[0081]

[0082] In the formula, [K R ] P-SV and are the local dynamic stiffness matrix in the P-SV plane of the half-space and the local dynamic stiffness value out of the plane, respectively. By superimposing the special solution of the external load amplitude on the interface of the half-space with the homogeneous solution, the total external load can be obtained:

[0083]

[0084] The upper interface of the underlying bedrock half-space is also the lower interface of the adjacent overburden medium. According to the interlayer continuity criterion, the external load acting on the upper interface of the half-space shown in Equation (29) will be applied in the opposite direction to the lower interface of the adjacent overburden, thereby causing a dynamic response of the particles in the overburden medium. Since there is no wave source in the overburden medium, the external load at each interface is zero. Therefore, by integrating the external load amplitudes at all interfaces of the layered medium, the external load amplitude vector in the overall cylindrical coordinate system is obtained:

[0085]

[0086] At this point, the displacement amplitude of each interface can be calculated according to the overall dynamic balance equation as shown below:

[0087]

[0088] Where, [K] P-SV with [K] SH are the overall dynamic stiffness matrices for in-plane and out-of-plane motion in the cylindrical coordinate system, respectively. They can be obtained by integrating the local dynamic stiffness matrices of the overlying N layers of soil medium and the underlying bedrock half-space, and they are exactly the same as the expressions in the plane coordinate system. The overall interface displacement obtained above is converted into the local displacement of the upper and lower interfaces of each soil layer medium, and then the amplitude of the upward and downward waves in the layer medium is calculated using the local dynamic equilibrium equation of each soil layer. Then, by back-substituting the amplitude and changing the coordinates, the free field displacement of any spatial position in the frequency-wavenumber domain can be obtained. Since this part of the free field displacement is caused by the reaction force shown in formula (29), it is called the reaction force solution. Among them, the reaction force solutions of the radial and vertical free field displacements in the j-th layer of soil medium induced by the unit amplitude P wave are respectively expressed as and The free-field displacement reaction solutions in the radial and vertical directions within the j-th soil layer induced by the unit amplitude SV wave are expressed as and The free field displacement reaction solution along the circumferential direction in the j-th soil layer induced by the unit amplitude SH wave is expressed as Based on the relationship between displacement-velocity and displacement-stress, the corresponding free-field velocity reaction force solution and free-field stress reaction force solution can be further derived. Finally, the free field in the frequency-wavenumber domain obtained above is converted back to the spatial domain using the inverse Hankel transform, and the body wave free field of the spherical diffusion body wave in the spatial domain is obtained. For the underlying bedrock half-space, the free field consists of three parts: the special solution, the homogeneous solution, and the reaction force solution. For the overlying layered soil layer, the free field consists only of the reaction force solution, and its displacement Green's function is shown in Equation (1).

[0089] Based on the above displacement Green's function, by multiplying it with the amplitude in the frequency domain, the body wave free field displacement at any spatial position near the fault site is obtained, as shown in formula (6). For the unknown amplitude, it can be obtained by inverting the Green's function solution at the surface station location shown in formulas (2)-(4) and the displacement record collected at that point, as shown in formula (5). Finally, the body wave free field displacement is calculated by performing time partial differential and spatial partial differential calculations respectively, constructing the body wave velocity free field and stress free field, and completing the body wave free field inversion near the fault site.

[0090] By extending the body wave dynamic stiffness matrix method in cylindrical coordinates to surface waves, the free-field displacement induced by the cylindrically diffused Rayleigh / Love surface waves can be derived, as shown in formulas (7)-(11). Based on the independence of horizontal radial and vertical motion in the free-field displacement of surface waves, inversion can be performed in two directions. In the vertical direction, surface waves have multimodal dispersion characteristics, which are decoupled by modal participation coefficients, and the displacement mode shape corresponding to each single mode is calculated. The displacement mode shape characterizes the ratio of the displacement amplitude at any vertical position to the surface, so the ratio of other elements to the first element can be defined as the vertical attenuation factor, as shown in formula (12). In the horizontal radial direction, based on the consistency of vertical motion, the horizontal radial displacement relationship expression between the surface station and any reference point can be established, which can be defined as the horizontal radial attenuation factor, as shown in formula (13). Furthermore, the single-mode surface wave displacement components collected at the surface station are multiplied by the horizontal radial attenuation factor and the vertical attenuation factor, respectively, to establish a single-mode surface wave displacement free field at any spatial position, as shown in formula (14). The time partial differential and spatial partial differential calculations of the single-mode surface wave free field displacement are performed, respectively, to construct the corresponding single-mode velocity free field and stress free field. Finally, according to the modal superposition principle, all single-mode surface wave free fields are superimposed to form a multi-mode surface wave free field, completing the surface wave free field inversion.

[0091] The above-mentioned body wave free field and surface wave free field are superimposed to form the multi-phase total free field near the fault site. However, in the process of body wave free field inversion, there are two unknown quantities in the Green function at the surface station location, namely, the epicenter distance r and the HO ” and “focal depth H”. The set of motion equilibrium equations shown in formula (5) is used to solve these two unknown quantities. However, the number of equations (N f =3) more than the number of unknowns (N x =2), which constitutes an overdetermined set of equations in mathematics. In addition, in the process of surface wave free field inversion, it is also necessary to solve the unknown “epicenter distance r HO In the system of motion equilibrium equations shown in formula (14), the number of equations (N f =3) is also greater than the number of unknowns (Nx =1), which also constitutes an overdetermined set of equations in mathematics. This means that the multi-phase free-field inversion near the fault site is essentially an overdetermined problem. At this time, the unknown quantity "epicenter distance r HO There is no exact analytical solution for "and" focal depth H".

[0092] Attachment Figure 1 The second step is to transform the overdetermined problem in the inversion framework of the first step into a multi-objective optimization problem and establish a multi-objective optimization model by constructing reasonable objective functions and constraints. The specific process is as follows:

[0093] In order to solve the problem of overdetermined equations without exact solutions, the least squares method is usually used in the field of mathematics to construct a constrained least squares problem, and then the least squares problem is transformed into an optimization problem, by finding an optimal solution to approximately replace the non-existent exact solution. Drawing on this solution idea, the method of the present invention transforms the overdetermined problem of near-fault free-field inversion into a problem of finding the optimal variable ("epicenter distance r HO In the optimization problem, the optimization goal is characterized by the objective function, and the convergence of the iteration is controlled by imposing constraints.

[0094] In the near-fault free-field inversion framework, the objective function is defined as: the multi-phase seismic records at any two points on the site obtained by forward modeling based on the optimized variables have the minimum error with respect to the target variable. Therefore, the key to constructing the objective function lies in determining the target variable that characterizes the seismic motion error between the two points. This target variable should be selected as a seismic motion index with similar characteristics at the two points. In addition, the selected target variable should also have a high sensitivity to the two optimized variables. Based on the above criteria, the method of the present invention selects the spectrum as the target variable from the three elements that characterize seismic motion characteristics (amplitude, spectrum, and duration). The spectrum is reflected as a waveform in the time domain, and thus the objective function can be constructed by minimizing the difference between the waveforms at two points. In order to solve the problem of phase difference between the waveforms at any two points, traditional difference measurement indicators, such as the root mean square error, cannot quantitatively calculate the difference between two time series with such phase difference. Based on this, the method of the present invention adopts the dynamic time warping algorithm (DTW) to define the warp path distance (WPD) as an indicator to quantify the waveform difference, which is used to characterize the difference between two time series with phase difference characteristics.

[0095] The dynamic time warping algorithm is an effective method to measure the difference between two time series, and is particularly suitable for situations where time series data have different time quantities and time offsets. The algorithm overcomes the limitations of traditional difference measurement methods when processing time series by allowing the time series to be warped on the time axis to find the optimal match between the two series. The core of the algorithm is to match the corresponding points in the two sequences by warping the time axis to find the best correspondence between them. Then, based on the best correspondence, a "warped" path is found so that the sum of the distances between all points on this path is minimized. This warped path is composed of non-overlapping elements in the cumulative distance matrix D and is quantified by the warped path distance. In the calculation process of WPD, the key step is to construct a two-dimensional cumulative distance matrix D. Taking the minimization of the displacement waveform difference between two close points in the site as the measurement criterion, the time series The displacements of the discrete points are arranged from left to right along the horizontal dimension, and the time course The displacements of each discrete point are arranged vertically from top to bottom. When the number of discrete points in the two displacement time histories is n, a two-dimensional matrix of n×n is formed. Each element in the matrix represents the cumulative distance between corresponding points in the two displacement time histories. The DTW regularized path begins at the top-left corner element (1,1) of the matrix and then traverses along its adjacent elements and diagonal elements, adhering to the continuity criterion. The continuity criterion requires that any two adjacent points (p1,q1) and (p2,q2) on the path must satisfy the conditions 0≤|p1–p2|≤1 or 0≤|q1–q2|≤1. Furthermore, the traversal process adheres to the monotonicity criterion, which requires that the two points (p1,q1) and (p2,q2) on the path must satisfy p2–p1≥0 or q2–q1≥0. Ultimately, the regularized path terminates at the bottom-right corner element (n,n).

[0096] In the cumulative distance matrix D, the cumulative distance corresponding to the first element value D(1,1) is equal to the time course The displacement and time history of the first discrete point in The Euclidean distance d of the first discrete point displacement 11 :

[0097]

[0098] Starting from the element position (1,1), there are three possible paths for the first traversal step: the horizontally adjacent element position (2,1), the vertically adjacent element position (1,2), and the diagonal adjacent element position (2,2). The final traversal direction will be along the path with the minimum value among the three element values D(2,1), D(1,2), and D(2,2). The element value D(2,1) represents the cumulative distance of the traversal path from the first element position to the second element position, which is given by and The base distance d between two points 21 And the cumulative distance D(1,1) of the previous traversal step:

[0099]

[0100] Similarly, D(1,2) is calculated as shown in formula (36):

[0101]

[0102] Without loss of generality, when the path traverses along the bottom or left end of the matrix, the cumulative distance to any element position (p,1) and (1,q) is:

[0103]

[0104] in,

[0105]

[0106] When the path traverses along the internal element positions of the matrix, the cumulative distance to any element position (p,q) is shown in formula (15). Next, the traversal path will continue to move toward one of the matrix elements (p+1,q), (p,q+1), (p+1,q+1) until it reaches the end position (n,n) of the path. Accordingly, the element value D(n,n) is the final cumulative distance along the regularized path, that is, the distance WPD of the two displacement time regularized paths. As shown in formula (16), the smaller the value of WPD, the shorter the traversal path, and thus the smaller the difference in displacement waveforms at two nearby points. Based on this, a multi-objective function is established based on minimizing the differences in displacement waveforms of two nearby points in three dimensions, as shown in formula (17).

[0107] Constraints are restrictions on the solution domain of optimization variables in an optimization problem and are usually expressed in the form of equations or inequalities. Considering the attenuation characteristics of body and surface wave propagation, the method of the present invention uses the displacement peak at the propagation lag point to be smaller than that at the first arrival point as a constraint, thus forming the inequality constraint shown in formula (18). In addition, the values of both optimization variables should be greater than 0, as shown in formula (19). Finally, the constraints are combined with the objective function to form the multi-objective optimization model shown in formula (20).

[0108] Attachment Figure 1 In the third step, the fast non-dominated multi-objective genetic algorithm (NSGA-II) was used to solve the optimization variables in the model (epicenter distance r HO ” and “focal depth H”) and substituted it into the inversion framework to complete the multi-phase free-field inversion of the near-fault site.

[0109] Here’s how it works:

[0110] In multi-objective optimization problems, each objective is usually constrained by each other, and there may be a large negative correlation between any two objectives. This means that in the optimization process, the improvement of one objective performance is often at the expense of the loss of other objective performances, and it is usually impossible to have a solution that makes all objective performances reach the optimal level. Therefore, for multi-objective optimization problems, the optimal solution is usually a balanced solution under multi-objective control, and the solution is not improveable, that is, it is impossible to change a certain target value of the solution without making other objectives worse. Based on this, the method of the present invention adopts a non-dominated sorting genetic algorithm (NSGA-II) based on simulated biological evolution to solve. This algorithm is a classic algorithm for solving multi-objective optimization problems. Its core is to coordinate the relationship between each optimization objective and find a balanced solution that satisfies each objective as much as possible. The resulting balanced solution is called a Pareto optimal solution. After all Pareto optimal solutions are integrated, a set of non-inferior solutions is formed - the Pareto solution set.

Claims

1. A multi-phase free-field inversion method for near-fault sites based on spherical / cylindrical diffusion characteristics, characterized by: Here are the steps: Step 1: Based on the spherical diffusion characteristics of body waves and the cylindrical diffusion characteristics of surface waves, a multi-phase free-field inversion framework for near-fault-layered media sites is established in the cylindrical coordinate system. Step 2: The unknown "epicenter distance r" in the multi-phase free-field inversion framework of the near-fault-stratified medium site HO " and "focal depth H" are used as optimization variables, and the inversion problem is transformed into a multi-objective optimization problem. A multi-objective function is constructed based on the criterion of minimizing the difference in three-dimensional displacement waveforms of two nearby points on the site. Constraints are set based on the propagation attenuation characteristics of seismic waves, and a multi-objective optimization model for near-fault inversion is established. Step 3: Use a fast non-dominated multi-objective genetic algorithm to solve the Pareto optimal solution of the optimization variables and substitute it into the inversion framework to complete the multi-phase free-field inversion of the near-fault site.

2. The multi-phase free-field inversion method for near-fault sites based on spherical / cylindrical diffusion characteristics according to claim 1 is characterized in that: The process of the first step is as follows: Step 1.1: Assume that the earthquake source is located in the underlying bedrock half-space, establish a cylindrical coordinate system r-θ-Z with the epicenter as the origin, the line connecting the epicenter and the surface station as the radial direction, and the direction from the epicenter to the earthquake source as the vertical direction; Step 1.2: Use the dynamic stiffness matrix method to calculate the Green's function of the spherical diffusion P wave, SV wave, and SH wave at any position (r, θ, z) in space: Where k is the horizontal radial wave number; and are the radial and vertical displacements of the j-th soil layer medium particle induced by the unit amplitude P wave in the wave number domain; and are the radial and vertical displacements of the j-th soil layer medium particle induced by the unit amplitude SV wave in the wave number domain; is the displacement of the particle in the jth soil layer along the circumferential direction induced by the unit amplitude SH wave in the wave number domain; are the spatial domain displacement Green's functions obtained by the inverse Hankel transform of the above-mentioned wavenumber domain displacement; Step 1.3: Take j = 1 in formula (1) and calculate the three types of spherical diffuse body waves at the surface station position (r HO ,0,0) at Green's function solution: In the formula, the Green function solution depends on the unknown quantity "epicenter distance r HO ” and “focal depth H”; Step 1.4: Use the normalized inner product method to extract the body wave displacement components collected at the surface station Based on formula (5), the incident amplitudes of P wave, SV wave and SH wave are inverted. Step 1.5: Multiply the incident wave amplitude obtained by formula (5) by the Green's function shown in formula (1) to construct the displacement free field induced by the spherical diffuser wave at any position in the spatial domain. Step 1.6: Calculate the time partial differential and spatial partial differential of the body wave displacement free field, construct the corresponding velocity free field and stress free field, and complete the body wave free field inversion; Step 1.7: Extend the body wave dynamic stiffness matrix method in cylindrical coordinates to surface waves and derive the spatial free-field displacement induced by the Rayleigh surface waves propagating on the cylindrical surface. and Where, and are the radial and vertical displacements of the particle in the jth soil layer induced by Rayleigh waves in the cylindrical coordinate system; and are the spatial domain displacements induced by Rayleigh waves in the plane coordinate system; i is the imaginary number, satisfying i 2 =-1; J α (kr) is the first kind α-order Bessel function; and are the up- and down-wave amplitudes of the non-uniform P wave in the j-th soil layer; and are the upgoing and downgoing wave amplitudes of the inhomogeneous SV wave in the medium layer, respectively; and are the vertical wave numbers of non-uniform P waves and non-uniform SV waves, respectively; is the velocity ratio of the P wave to the Rayleigh surface wave in the jth layer of medium; is the velocity ratio of the S wave to the Rayleigh surface wave in the jth layer of medium; Step 1.8: Similarly, derive the free-field displacement induced by the Love surface wave diffused by the cylinder. Where, is the circumferential displacement of the particle in the j-th soil layer induced by the Love wave in the cylindrical coordinate system; is the spatial domain displacement induced by the Love wave in the plane coordinate system; and are the up- and down-wave amplitudes of the non-uniform SH wave in the j-th soil layer medium; is the vertical wave number of the non-uniform SH wave; Step 1.9: Perform modal decoupling on the dispersive surface wave. Based on the single-mode displacement mode shape, use formula (12) to derive the vertical attenuation factor of the single-mode surface wave free-field displacement: Where, the subscript m represents the mth mode of the dispersive surface wave; Step 1.10: Construct the surface station position (r HO ,0,0) and the single-mode surface wave free field at any reference point (r,0,0) on the surface; further based on the consistency criterion of the vertical motion of the two points, the attenuation factor of the single-mode surface wave free field displacement along the horizontal radial direction is derived using formula (13): Step 1.11: Use the normalized inner product method to extract the surface wave displacement components collected at the surface station, and convert the decoupled single-order modal displacement Multiply the corresponding horizontal radial attenuation factor and vertical attenuation factor respectively to construct the single-mode surface wave displacement free field at any position (r, θ, z) in space: Where, the free field displacement of a single-mode surface wave depends on the unknown quantity "epicenter distance r HO ”; Step 1.12: Calculate the time partial differential and spatial partial differential of the single-mode surface wave displacement free field to construct the corresponding single-mode velocity free field and stress free field; Step 1.13: Based on the principle of surface wave dispersion mode superposition, all single-mode free fields are superimposed to form a multi-mode free field, completing the surface wave free field inversion. Step 1.14: Superimpose the body wave free field and the surface wave free field to construct the multi-phase total free field near the fault site.

3. The multi-phase free-field inversion method for near-fault sites based on spherical / cylindrical diffusion characteristics according to claim 1, characterized in that: The second step is specifically performed as follows: Step 2.1: The waveform difference at similar locations near the fault site is used as the target variable of the objective function, and the dynamic regularized path distance is used as a quantitative indicator to measure the waveform difference with phase difference characteristics; Step 2.2: Construct the multi-seismic phase displacement time history of two points close to the site and The cumulative distance matrix D is composed of: Where p and q represent the number of rows and columns of the matrix respectively, and D(p,q) is the matrix element at that position; d pq For timeline The displacement and time history of the pth discrete point in The Euclidean distance of the displacement of the qth discrete point in ; Step 2.3: Let the regular path traverse the cumulative distance matrix D, starting from the upper left corner element position (1,1) of the matrix and ending at the last element (n,n). Calculate the final cumulative distance along the regular path: Where n is the total number of discrete points in the two time histories; the smaller the WPD value, the shorter the traversal path, and thus the smaller the difference in displacement waveforms between two adjacent points; Step 2.4: Based on formula (16), a multi-objective function is constructed with the criterion of minimizing the difference in three-dimensional displacement waveforms between two adjacent points on the site: Where, and is the radial multi-seismic phase displacement of two nearby points on the site; and is the multi-seismic phase displacement in the tangential direction; and is the multi-seismic phase displacement in the vertical direction; Step 2.5: Considering the attenuation characteristics of seismic waves during propagation, the displacement amplitude at the propagation lag point is Smaller than the displacement amplitude at the first point As constraints: Step 2.6: Keep the values of the two optimization variables greater than 0 to further narrow the feasible region: r HO ∈R + ,H∈R + (19) Step 2.7: Combine the objective function in formula (17) with the constraints shown in formulas (18)-(19) to establish a multi-objective optimization model for near-fault inversion:

Images