Free-field inversion method of near-fault site multiple seismic phase based on spherical / cylindrical diffraction characteristics

Abstract

The application provides a near-fault site multi-seismic-phase free-field inversion method based on spherical / cylindrical diffusion characteristics, and improves the research system in the near-fault inversion field. The application is based on the spherical diffusion characteristics of body waves and the cylindrical diffusion characteristics of surface waves, adopts a dynamic stiffness matrix method to construct a near-fault site multi-seismic-phase free-field inversion framework, and converts the inversion problem into a multi-objective optimization problem of solving the "epicentral distance" and "focal depth"; on this basis, the dynamic regular path distance is used to quantify the space-time difference of the free-field displacement waveform, and a target function is constructed by taking the minimization of the three-dimensional displacement waveform difference of two points close to each other as the criterion; meanwhile, a constraint condition is set in combination with the propagation and attenuation characteristics of seismic waves, so as to establish a near-fault inversion multi-objective optimization model; subsequently, a fast non-dominated multi-objective genetic algorithm is introduced to solve the Pareto optimal solution of the optimization variable, and the optimal solution is substituted into the inversion framework to complete the near-fault site multi-seismic-phase free-field inversion.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of wave field inversion of seismic waves, and relates to a near-fault site multi-seismic phase free field inversion method based on spherical / cylindrical diffusion characteristics. BACKGROUND

[0002] Under the action of earthquakes, these structures and the surrounding soil have non-negligible dynamic interaction, forming a typical soil-structure interaction (SSI) system. Seismic response analysis considering SSI effect is the premise of accurately assessing the seismic risk of infrastructure, and this process relies on a reasonable seismic input method and accurate seismic excitation. Under the framework of the existing boundary-equivalent load seismic input method, the unknown underground free field forms the seismic excitation source of the SSI system, which needs to be obtained by inversion through the known ground motion and site conditions. The accuracy of the inverted underground free field will directly affect the rationality of the subsequent structure seismic input and the reliability of the response evaluation, and further affect the safety of the structure design.

[0003] In the existing free-field inversion method, whether the body wave propagating in the tilt direction or the surface wave propagating in the horizontal direction, the assumption of plane wave is adopted, that is, the wave line direction remains consistent, and the wave front perpendicular to the wave line is a plane. However, the assumption of plane wave is only applicable to the far-fault region far away from the epicenter. In the near-fault region, the body wave and the surface wave respectively exhibit the characteristics of real spherical diffusion and cylindrical diffusion. Specifically, the body wave generated by the fault dislocation will spread in the form of spherical wave from the source position to the surrounding, the wave line direction is different, the wave front presents spherical shape, and the energy attenuation is inversely proportional to the square of the propagation distance. The surface wave induced will spread in the form of cylindrical wave to the surrounding, the wave line direction is also different, the wave front presents cylindrical shape, and the energy attenuation is inversely proportional to the square root of the propagation distance. In particular, the spherical wave and the cylindrical wave attenuate in all directions, and the equal energy line is a circle composed of the same spherical radius and cylindrical radius. In comparison, the energy of the plane wave only attenuates along its propagation direction, and the equal energy line is a line perpendicular to the propagation direction. The above differences in attenuation characteristics lead to different free-field responses in the near-fault and far-fault regions, so that the far-fault free-field inversion method based on the assumption of plane wave is difficult to reflect the free-field characteristics in the near-fault region. More importantly, the near-fault region will generate special pulse-type earthquakes due to the directionality effect and the slip-off effect. The long-period component of the pulse-type earthquake is significant, and the energy carried by the pulse-type earthquake is mainly concentrated in the low frequency band. This low frequency characteristic makes the long-period structure in the near-fault region prone to "resonance effect" with the pulse-type earthquake, thereby exacerbating the damage and destruction of the structure. The SSI system usually has a long period, and the existing far-fault inversion method is difficult to accurately reflect the long-period pulse characteristics of the near-fault pulse-type ground motion, which will directly affect the seismic response results of the structure. Therefore, it is urgent to propose a free-field inversion method suitable for the near-fault site, to provide accurate excitation for the seismic analysis of the SSI system located in the near-fault region, and to provide reliable guidance for the structure response evaluation and seismic design. SUMMARY

[0004] The application provides a near-fault site multi-phase free-field inversion method based on spherical / cylindrical diffusion characteristics, which perfects the research system in the near-fault inversion field.

[0005] The technical scheme of the application is as follows:

[0006] Based on the spherical diffusion characteristics of body waves and the cylindrical diffusion characteristics of surface waves, a near-fault site multi-phase free-field inversion framework is constructed by using the dynamic stiffness matrix method, and the inversion problem is converted into a multi-objective optimization problem of solving the "epicentral distance" and "focal depth". On this basis, the dynamic regularization path distance is used to quantify the spatial and temporal differences of free-field displacement waveforms, and the target function is constructed by minimizing the three-dimensional displacement waveform differences of two nearby points. At the same time, combined with the propagation and attenuation characteristics of seismic waves, the constraint condition is set, and a near-fault inversion multi-objective optimization model is established. Then, the Pareto optimal solution of the optimization variables in the model is solved by introducing the fast non-dominated sorting genetic algorithm (NSGA-II), and the near-fault site multi-phase free-field inversion is completed by substituting it into the inversion framework.

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

[0008] First step: Based on the spherical diffusion characteristics of body waves and the cylindrical diffusion characteristics of surface waves, a near-fault layered medium site multi-phase free-field inversion framework is established in the cylindrical coordinate system;

[0009] Second step: The unknown "epicentral distance r HO " and "focal depth H" in the near-fault layered medium site multi-phase free-field inversion framework are taken as optimization variables, and then the inversion problem is converted into a multi-objective optimization problem. The multi-objective function is constructed by minimizing the three-dimensional displacement waveform differences of two nearby points, and the constraint condition is set by the propagation and attenuation characteristics of seismic waves to establish a near-fault inversion multi-objective optimization model;

[0010] Third step: The Pareto optimal solution of the optimization variables is solved by using the fast non-dominated sorting genetic algorithm, and the near-fault site multi-phase free-field inversion is completed by substituting it into the inversion framework.

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

[0012] Step 1.1, let the source be located in the underlying bedrock half space, establish a cylindrical coordinate system r-θ-Z with the epicenter position as the origin, the direction of the line connecting the epicenter and the ground station as the radial direction, and the epicenter pointing to the source direction as the vertical direction;

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

[0014]

[0015] In the formula, k is the horizontal radial wave number; and are the wave number domain displacements of the jth layer of soil medium particles induced by unit amplitude P wave along the radial and vertical directions, respectively; with are the wave number domain displacements of the jth layer of soil medium particle along the radial and vertical directions induced by unit wave amplitude SV wave, respectively; is the wave number domain displacement of the jth layer of soil medium particle along the circumferential direction induced by unit wave amplitude SH wave; are the spatial domain displacement Green's functions of the above wave number domain displacements obtained by inverse Hankel transform;

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

[0017]

[0018] In the formula, the Green's function solution depends on the unknown quantity "epicentral distance r HO " and "focal depth H";

[0019] Step 1.4, the body wave displacement components collected at the surface station are extracted by using the normalized inner product method and the incident wave amplitudes of P wave, SV wave and SH wave are inverted based on formula (5)

[0020]

[0021] Step 1.5, the incident wave amplitudes obtained by formula (5) are multiplied by the Green's function shown in formula (1) respectively, and the displacement free field at any position in the spatial domain induced by spherical diffuse body waves is constructed

[0022]

[0023] Step 1.6, the time partial differential and space partial differential of the body wave displacement free field are calculated respectively to construct the corresponding velocity free field and stress free field, and the body wave free field inversion is completed;

[0024] Step 1.7, the body wave dynamic stiffness matrix method in cylindrical coordinate system is extended and applied to surface wave, and the spatial domain free field displacement induced by cylindrical diffuse Rayleigh surface wave is derived with

[0025]

[0026] In the formula, with are the wave number domain displacements of the jth layer of soil medium particle along the radial and vertical directions induced by Rayleigh wave in cylindrical coordinate system, respectively; and respectively, are the spatial-domain displacements induced by Rayleigh waves in the planar coordinate system; i is the imaginary unit, satisfying i 2 = -1; J α (kr) is the first kind of Bessel function of order a; and and and and and and is the wave velocity ratio of P wave to Rayleigh surface wave in the jth layer of medium; is the wave velocity ratio of S wave to Rayleigh surface wave in the jth layer of medium;

[0027] Step 1.8, similarly, the free-field displacement induced by Love surface waves diffused by cylindrical surface is derived as

[0028]

[0029] where, is the wave-number-domain displacement of the jth layer of medium induced by Love waves along the circumferential direction in the cylindrical coordinate system; is the spatial-domain displacement induced by Love waves in the planar coordinate system; and and is the vertical wave number of non-uniform SH wave;

[0030] Step 1.9, the modal decoupling of frequency-dispersive surface waves is performed, and based on the single-mode displacement modal shape, the attenuation factor of single-mode surface wave free field along the vertical direction is derived by using formula (12):

[0031]

[0032] wherein subscript m represents the mth mode of frequency-dispersive surface waves;

[0033] Step 1.10, the single-mode surface wave free field at the station position (r HO , 0, 0) on the ground surface and at any reference point (r, 0, 0) on the ground surface is constructed according to formulas (7)-(11); further based on the consistency criterion of two-point vertical motion, the attenuation factor of single-mode surface wave free field displacement along the horizontal radial direction is derived by 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] In the formula, n is the total number of discrete points of two time periods; the smaller the value of WPD is, the shorter the traversal path is, and the smaller the displacement waveform difference between the two adjacent points is;

[0049] Step 2.4, based on formula (16), the three-dimensional displacement waveform difference between the two adjacent points in the field is minimized as the criterion to construct a multi-objective function:

[0050]

[0051] In the formula, and is the multi-seismic phase displacement of the two adjacent points in the field in the radial direction; 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 the seismic wave in the propagation process, the displacement amplitude of the propagation lag point is less than the displacement amplitude of the first arrival point as a constraint condition:

[0053]

[0054] Step 2.6, the values of the two optimization variables are kept greater than 0, and the feasible domain range is further narrowed:

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

[0056] Step 2.7, the multi-objective optimization model of near-fault inversion is established by combining the objective function in formula (17) and the constraint conditions shown in formula (18)-(19):

[0057]

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

[0059] (1) The present application provides a near-fault site multi-seismic phase free field inversion method based on spherical / cylindrical diffusion characteristics, which fully considers the spherical diffusion characteristics of body waves and the cylindrical diffusion characteristics of surface waves, makes up for the applicability limitations of the existing far-fault inversion method in near-fault sites, and greatly improves the inversion accuracy. The comparative analysis results show that the free field results obtained by inversion are highly consistent with the true borehole records, thereby providing more reliable seismic excitation for the seismic response analysis of the SSI system, and having important engineering application for revealing the real seismic response behavior, damage mechanism and failure mode of the structure.

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

[0061] The long-period pulse component in pulse-type seismic motion is the main cause of the damage of medium and long-period structures in the near-fault region. In particular, when the pulse period is close to the natural period of the structure, resonance effect will be easily triggered. Therefore, the determination of pulse characteristics and the calculation of pulse period are of great significance to the seismic design of medium and long-period structures, especially long-period SSI systems. When the near-fault ground 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 evaluation and seismic design of SSI systems in the near-fault region. The method can truly reflect the pulse characteristics of underground seismic motion according to the pulse-type ground motion and local site conditions, thereby providing accurate excitation for the seismic analysis of SSI systems, and has important engineering application value. The formula derivation result is reliable, and the calculation process is efficient.

[0062] (3) The method converts the overdetermined problem in the free field inversion of near-fault sites into an optimization problem, and by introducing a multi-objective optimization theory framework, the solution of the overdetermined equation set without exact solution is converted into a parameter optimization problem with clear physical meaning, and the balanced solution of the optimization variables under multi-objective control is obtained.

[0063] In view of the problem that the overdetermined equation set has no exact solution, the mathematical field usually constructs a constrained least squares problem, and then converts the least squares problem into an optimization problem, and finds an optimal solution to approximate the non-existent exact solution. By analogy with this solution idea, the method converts the overdetermined problem of near-fault free field inversion into an optimization problem of finding the optimal solution of "epicentral distance r HO " and "source depth H". This idea has clear physical meaning, sufficient theoretical basis, and accurate and objective results, and provides a scientific analysis method for near-fault seismic inversion research. BRIEF DESCRIPTION OF DRAWINGS

[0064] Figure 1 The flowchart of the method is shown in the figure. DETAILED DESCRIPTION

[0065] The application will be further described below in combination with the drawings in the embodiments of the application.

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

[0067] According to the wave theory, the frequency domain dynamic stiffness matrix of the site in the cylindrical coordinate system is consistent with that in the plane coordinate system, so the frequency domain dynamic stiffness matrix method based on the plane wave assumption can be extended to the free field induced by the cylindrical wave. In addition, the free field induced by the spherical wave can be transformed into the superposition of the cylindrical wave free field in the frequency-wave number domain, so as to realize the construction of the spherical wave and the cylindrical wave free field respectively. Based on the above principle, Ballentine et al. proposed the Green function solution of spherical P wave, SV wave and SH wave in three-dimensional layered medium under the assumption of point source. The general idea is as follows: firstly, the fixed layer solution and the fixed end surface reaction are solved, then the fixed end surface reaction is applied to other layers of the layered medium site to solve the reaction solution, and finally the total Green function solution is obtained by superimposing the fixed layer solution and the fixed end surface reaction solution.

[0068] The method of the application learns from the idea, and the source is fixed in the underlying bedrock half space of the layered medium site, and the body wave free field and the surface wave free field inversion framework are constructed. In the underlying bedrock half space, the spatial domain potential function of the spherical P wave, SV wave and SH wave induced by the point source can be expanded into the integral superposition of the cylindrical wave potential function in the frequency-wave number domain by using Hankel transformation:

[0069]

[0070] In the formula, respectively, the potential function of the spatial domain spherical P wave, SV wave and SH wave; respectively, the incident wave amplitude of the three types of spherical body waves; k P and k S respectively, the wave number of P wave and S wave, related to the complex compression wave velocity and the complex shear wave velocity of the half space site, satisfy and R is the seismic wave propagation distance, that is, the distance between the propagation position and the source. When transformed from the spatial domain to the frequency-wave number domain, the seismic wave propagation distance R can be expressed as the square root of the square sum of the horizontal radial propagation distance r and the vertical propagation distance z; k represents the horizontal radial wave number in the frequency-wave number domain; J0(kr) is the 0th Bessel function (the first kind); and respectively, the vertical wave number of P wave and S wave in the frequency-wave number domain, satisfy:

[0071]

[0072] According to formula (21), the cylindrical wave component potential function of the spherical P wave, SV wave and SH wave under the unit incident wave amplitude is respectively:

[0073]

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

[0075]

[0076] where, and denote the free-field displacement amplitudes in radial, vertical and circumferential directions in frequency-wave number domain, respectively. It is worth noting that P wave and SV wave are coupled in P-SV plane, thus when calculating the free-field displacement induced by P wave and SV wave respectively, the potential function of the other one should be set to 0. For example, when solving the free-field displacement induced by P wave independently, let On the contrary, when solving the free-field displacement induced by SV wave independently, let Based on the free-field displacement obtained from the above formula, the free-field stress induced by unit amplitude body wave can be further calculated according to the stress-displacement relationship, as shown below:

[0077]

[0078] where, λ and μ are the Lame constants of the underlying half-space. In particular, when the distance between the upper interface of the half-space and the source depth is d R , substituting z = -d R into formula (24)-(25) can calculate the displacement amplitude and stress amplitude of the particles in the soil layer on the upper interface of the half-space. This part of the solution is a special solution, and the superscript's' represents:

[0079]

[0080] At this time, the amplitude of the external load on the upper interface of the half-space satisfies and Substituting the special solution of the displacement amplitude of the upper interface of the half-space into the local dynamic equilibrium equation of the half-space, the homogeneous solution of the amplitude of the external load on the upper interface of the half-space, i.e. the general solution, can be obtained, as shown in formula (27) and (28). The homogeneous solution is denoted by the superscript 'h'.

[0081]

[0082] where, [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. Superimposing the special solution and the homogeneous solution of the amplitude of the external load on the upper interface of the half-space, 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] On the basis of the displacement Green's function, the free-field displacement of body waves at any spatial position near the fault is obtained by multiplying the displacement Green's function with the amplitude in the frequency domain, as shown in equation (6). For the unknown amplitude, the displacement record collected at the surface station can be used to obtain the unknown amplitude by inversion, as shown in equation (5). Finally, the free-field displacement of body waves is differentiated with respect to time and space to calculate the free-field velocity and stress, and the inversion of the free-field displacement of body waves near the fault is completed.

[0090] The dynamic stiffness matrix method of body waves in the cylindrical coordinate system is extended to surface waves, and the free-field displacement induced by the Rayleigh / Love surface waves diffused from the cylindrical surface is derived, as shown in equations (7)-(11). Based on the independence of the horizontal radial and vertical motion of the free-field displacement of surface waves, the inversion is performed in two directions. In the vertical direction, the surface waves have multi-modal dispersion characteristics, and the displacement modal shape corresponding to each single mode is calculated by decoupling the modal participation coefficient. The displacement modal shape represents the displacement amplitude ratio at any vertical position to the surface, and therefore the ratio of other elements to the first element can be defined as the vertical attenuation factor, as shown in equation (12). In the horizontal radial direction, based on the consistency of the vertical motion, the relationship expression between the horizontal radial displacement at the surface station and at any reference point is established, which can be defined as the horizontal radial attenuation factor, as shown in equation (13). Further, the single-mode surface wave displacement component collected at the surface station is multiplied by the horizontal radial attenuation factor and the vertical attenuation factor, respectively, to establish the single-mode surface wave free-field displacement at any spatial position, as shown in equation (14). The single-mode velocity free-field and stress free-field are constructed by differentiating the single-mode surface wave free-field displacement with respect to time and space, respectively. Finally, according to the modal superposition principle, all single-mode surface wave free-fields are superimposed to form a multi-modal surface wave free-field, and the inversion of the surface wave free-field is completed.

[0091] The total free-field of the multi-phase near the fault site is formed by superimposing the above-mentioned body wave free-field and surface wave free-field. However, in the inversion process of the body wave free-field, there are two unknown quantities at the position of the surface station, i.e. "epicentral distance r HO " and "focal depth H". The motion balance equation shown in equation (5) is used to solve these two unknown quantities. However, the number of equations (N f = 3) is greater than the number of unknowns (N x = 2), which constitutes an over-determined equation in mathematics. In addition, in the inversion process of the surface wave free-field, the unknown "epicentral distance r HO " also needs to be solved. In the motion balance equation shown in equation (14), the number of equations (N f = 3) is also greater than the number of unknowns (Nx = 1), which is also a system of over-determined equations in mathematics. This means that the near-fault site multi-phase free-field inversion is essentially an over-determined problem, and there is no exact analytical solution for the unknowns "epicentral distance r HO " and "focal depth H".

[0092] The present application is not limited by the specific embodiments described herein. Figure 1 The second step is to transform the over-determined problem existing in the first step inversion framework into a multi-objective optimization problem, and to establish a multi-objective optimization model by constructing a reasonable objective function and constraint condition. The specific process is as follows:

[0093] In view of the problem that there is no exact solution to the system of over-determined equations, the mathematical field usually constructs a constrained least squares problem by using the least squares method, and then converts the least squares problem into an optimization problem, and finds an optimal solution to approximate the non-existent exact solution. By using the solving idea, the present application transforms the over-determined problem of the near-fault free-field inversion into an optimization problem of finding the optimal solution of the optimization variables ("epicentral distance r HO " and "focal depth H"). In the optimization problem, the optimization objective is represented by the objective function, and the convergence of iteration is regulated by imposing constraint conditions.

[0094] Under the framework of near-fault free-field inversion, the objective function is defined as: based on the optimization variables, the multi-phase seismic record at any two points of the site has the smallest error with respect to the target variable. It can be seen that the key to constructing the objective function is to determine the target variable representing the error of the seismic motion of the two points. The target variable should be selected as the seismic motion index with similar characteristics at the two points. In addition, the selected target variable should have high sensitivity to the two optimization variables. Based on the above criteria, the present application selects the frequency spectrum as the target variable from the three elements (amplitude, frequency spectrum, and duration) representing the characteristics of the seismic motion. The frequency spectrum is reflected in the time domain as a waveform, so the objective function can be constructed by minimizing the difference between the waveforms at the two points. In view of the problem that there is a phase difference between the waveforms at any two points, the traditional difference measurement index, such as root mean square error, cannot quantitatively calculate the difference between the two time series with phase difference. Based on this, the present application uses the dynamic time warping algorithm (Dynamic Time Warping, DTW) to define the warping path distance (Warp Path Distance, WPD) as an index for quantifying the difference between the waveforms, which is used to represent the difference between the two waveforms with phase difference characteristics.

[0095] Dynamic time warping (DTW) is an effective method for measuring the difference between two time series, especially when the time series have different time scales and time shifts. The algorithm finds the optimal match between the two sequences by allowing the time series to be warped in time, overcoming the limitations of traditional difference measures when dealing with time series. The core of the algorithm is to match corresponding points in the two sequences by warping the time axis, and then find the best correspondence between them. Then, based on the best correspondence, a "warping" path is found that minimizes the sum of distances between all points on the path. This warping path is composed of non-overlapping elements in the cumulative distance matrix D, and is quantified by the regularized path distance. In the calculation 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 points in the field as the measure, the time history The displacement of each discrete point is arranged from left to right along the horizontal dimension, and the time history The displacement of each discrete point is arranged from top to bottom along the vertical dimension. When the number of time discrete points of the two displacement time histories is n, an n x n two-dimensional matrix is formed. The value of each element in the matrix is the cumulative distance of the corresponding points in the two displacement time histories. The DTW warping path starts from the top left corner element position (1,1), then traverses the adjacent elements and diagonal elements, and follows the continuity criterion. The continuity criterion requires that any two adjacent points (p1, q1) and (p2, q2) on the path must satisfy the condition 0 ≤ |p1-p2| ≤ 1 or 0 ≤ |q1-q2| ≤ 1. In addition, the traversal process also follows the monotonicity criterion, i.e. the two points (p1, q1) and (p2, q2) on the path must satisfy p2-p1 ≥ 0 or q2-q1 ≥ 0. Finally, the warping path will stop at the bottom right corner element (n, n).

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

[0097]

[0098] Starting from the element position (1,1), there are three possibilities for the first traversal step, namely the horizontal adjacent element position (2,1), the vertical adjacent element position (1,2), and the diagonal adjacent element position (2,2). The final traversal direction will be along the path direction of the minimum value of 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 calculated by and Base distance d between two points 21 and the accumulated distance D(1,1) of the previous traversal step:

[0099]

[0100] Similarly, the calculation of D(1,2) is shown in equation (36):

[0101]

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

[0103]

[0104] where,

[0105]

[0106] When the path is traversed along the internal element positions of the matrix, the accumulated distance to reach any element position (p,q) is shown in equation (15). Next, the traversal path will continue to travel to one of the matrix elements (p+1,q), (p,q+1), and (p+1,q+1) until reaching the end position (n,n) of the path. Accordingly, the element value D(n,n) is the final accumulated distance along the regular path, that is, the two displacement time history regular path distance WPD. As shown in equation (16), the smaller the value of WPD, the shorter the traversal path, and thus the smaller the difference between the displacement waveforms of the two nearby points. Based on this, a multi-objective function is established according to the minimization of the displacement waveform difference of the two nearby points in three dimensions, as shown in equation (17).

[0107] The constraint condition is a limitation on the solution domain range of the optimization variable in the optimization problem, which is usually expressed in the form of an equation or an inequality. Considering the attenuation characteristics in the propagation process of body waves and surface waves, the present method takes the displacement peak value at the propagation lag point to be less than that at the first arrival point as a constraint condition, thereby forming the inequality constraint shown in equation (18). In addition, the values of the two optimization variables should both be greater than 0, as shown in equation (19). Finally, the constraint condition and the objective function together form the multi-objective optimization model shown in equation (20).

[0108] Appendix Figure 1 The third step uses the fast non-dominated multi-objective genetic algorithm (NSGA-II) to solve the Pareto optimal solution of the optimization variables ("epicentral distance r HO " and "focal depth H") in the model, and substitutes it into the inversion framework to complete the near-fault site multi-phase free-field inversion.

[0109] The principle is as follows:

[0110] In a multi-objective optimization problem, each objective is usually mutually constrained, and there can be a large negative correlation between any two objectives. This means that in the optimization process, the improvement of the performance of one objective is often at the expense of the performance of other objectives, and it is usually impossible to have a solution that optimizes all objectives. Therefore, for a multi-objective optimization problem, the optimal solution is usually a balanced solution under multi-objective control, and this solution cannot be improved, that is, it cannot be changed without making other objectives worse. Based on this, the method of the application uses a non-dominated sorting genetic algorithm based on simulated biological evolution (NSGA-II) for solution. This algorithm is a classic algorithm for solving multi-objective optimization problems, and 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, and the set of all Pareto optimal solutions forms a non-inferior solution set, the Pareto solution set.

Claims

1. A multi-phase free-field inversion method for near-fault sites based on spherical and 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: A fast non-dominated multi-objective genetic algorithm is used to find 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; 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: (1) ； 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: , P wave (2); , SV wave (3); , SH wave (4); 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 、 、 , and invert the incident amplitudes of P waves, SV waves, and SH waves based on formula (5) 、 、 : (5) ； 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. 、 、 : (6) ； 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 : (7) ； (8) ； (9) ； 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 amplitudes of the inhomogeneous SV wave in the medium; ks̃ Ray,j and kt̃ Ray,j 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. : (10) ； (11) ； 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 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; kt̃ Love,j 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: (12) ； 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): (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: (14) ； 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.

2. The multi-phase free-field inversion method for near-fault sites based on spherical and 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: (15) ； 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: (16) ； 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: (17) ； 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 converted to | | max Less than the displacement amplitude at the first point| | max As constraints: (18) ； Step 2.6: Keep the values ​​of the two optimization variables greater than 0 to further narrow the feasible region: (19) ； Step 2.7: Combine the objective function in formula (17) with the constraints shown in formulas (18) and (19) to establish a multi-objective optimization model for near-fault inversion: (20)。

Images