A Method and System for Constructing a 3D Model of Seismogenic Fault Based on Aftershock Precision Location Data

Abstract

This invention belongs to the field of seismic information processing technology and discloses a method and system for constructing a three-dimensional model of seismogenic faults based on aftershock precision location data. The method includes: acquiring seismic geological data of the study area; establishing a high-resolution three-dimensional crustal velocity model of the study area and determining precise source parameters; obtaining tectonic stress field characteristics through source mechanism solutions and GNSS observation data; and studying the seismogenic mechanism of the study area and obtaining information on potential strong earthquake hazards. This invention combines two grid intervals, 0.2°×0.2° and 0.1°×0.1°, to obtain both the velocity structure of the non-source area with a 0.2°×0.2° grid interval and the velocity structure of the source area with a 0.1°×0.1° grid interval. This is suitable for studying the velocity structure of source areas with a relatively small number of stations.

Description

Technical Field

[0001] This invention belongs to the field of earthquake information processing technology, and in particular relates to a method and system for constructing a three-dimensional model of the seismogenic fault based on aftershock precision location data. Background Technology

[0002] A certain city in a certain province is located at the southwestern end of the Shaoguang Depression in the Earth's dome system, at the southern end of the western edge of the second giant uplift zone in the tectonic system. Its seismic geological structure is relatively complex, making it one of the most seismically active areas along the coastal belt. [The remaining text appears to be unrelated and refers to a specific region, M.] L Earthquakes of magnitude 3.0 or higher can serve as a seismic "window" into coastal seismic zones. Therefore, studying the subsurface structure and seismic activity in a region is of great significance for assessing the seismic hazard of economically developed coastal areas.

[0003] Based on the above analysis, the problems and shortcomings of existing technologies are as follows: In existing seismic analysis, the prediction results for the three-dimensional crustal P-wave velocity structure and the precise location of earthquakes are low. Therefore, they cannot provide a theoretical basis for seismic hazard assessment. Summary of the Invention

[0004] To overcome the problems existing in related technologies, the present invention discloses a method and system for constructing a three-dimensional model of the seismogenic fault based on aftershock precision location data.

[0005] The technical solution is as follows: A method for constructing a three-dimensional model of a seismogenic fault based on aftershock precision location data, comprising the following steps:

[0006] S1, collecting earthquake and geological data information of the study area;

[0007] S2, establish a high-resolution three-dimensional crustal velocity model of the study area and determine the precise source parameters;

[0008] S3, the characteristics of the tectonic stress field are obtained through focal mechanism solutions and GNSS observation data;

[0009] S4 predicts the earthquake-generating mechanism, seismogenic structure, and source rupture process in the study area. By combining velocity structure, fault distribution, source mechanism, velocity field characteristics, and post-earthquake field investigation data, it predicts the influence of faults on the earthquake-generating process, thereby obtaining information on potential strong earthquake hazards.

[0010] In step S1, earthquake geological data information of the study area is collected, including: collecting earthquake phase reports, earthquake waveform data, GNSS observation data to establish a basic database, and conducting field surveys and reviewing previous data to collect geological structure data of earthquake-prone areas in the study area.

[0011] In step S2, a high-resolution three-dimensional crustal velocity model of the study area is established to determine precise source parameters, including:

[0012] S201, using a grid interval of 0.2°×0.2°, double-difference seismic tomography tomoDD inversion was performed to obtain the velocity structure MODEL_0.2 of the study area and surrounding areas;

[0013] S202, linear interpolation is performed on the velocity structure MODEL_0.2 of the study area and surrounding areas to obtain a velocity model MODEL_0.1 with a grid spacing of 0.1°×0.1°;

[0014] S203 uses the velocity model MODEL_0.1 with a grid interval of 0.1°×0.1° as the initial model. The tomoDD double-difference seismic tomography method is used for inversion calculation to obtain the final three-dimensional crustal P-wave velocity structure MODEL_final and the location results of multiple earthquakes in the study area.

[0015] In step S3, the characteristics of the tectonic stress field are obtained through focal mechanism solutions and GNSS observation data, including:

[0016] S301, obtain the focal mechanism solutions for multiple earthquakes;

[0017] S302, based on focal mechanism deconstruction, studies the regional tectonic stress field;

[0018] S303, obtains velocity field and surface dilatation rate field through GNSS data;

[0019] S304, fault plane parameter inversion: Based on the precise source parameter information of 6255 relocated earthquakes obtained by the double-difference seismic tomography method, and in accordance with the principle that cluster earthquakes occur near faults, the fault plane parameters of a seismic strip near a certain study area are inverted by combining the simulated annealing algorithm and the Gauss-Newton algorithm, to obtain the fault strike, dip angle, dip direction (NW), fault length, and fault depth.

[0020] In step S302, based on the focal mechanism deconstruction, the regional tectonic stress field is studied, including:

[0021] Based on magnitude, a weighted matching process is performed on the focal mechanism solutions, combining the focal mechanisms of large and small earthquakes. The given weights are derived from the following formula:

[0022] ω=e r / D 2

[0023] Set the weight ω of the earthquake with the smallest magnitude in ML3.0. min The weight ω of the largest earthquake, ML6.6, is 1.0. max The magnitude is 5, where r is the relative magnitude, r = mm. min ;

[0024] The magnitude attenuation coefficient D is derived from the maximum and minimum regional magnitudes, and its expression is:

[0025] D = (MM) min ) / ln(ω max )

[0026] In the formula, ω is the given weight, e r A given coefficient for source attenuation, D is the magnitude attenuation coefficient, ω max The maximum magnitude weight is M, where M is the current magnitude source. min The smallest magnitude earthquake source;

[0027] The D value varies in different regions, and the D value in the study area is 2.2368. In the actual inversion, the search interval of the three rotation angles of the stress field rotation axis is set to 1.0°, the search interval of the stress ratio is set to 0.1, and the confidence level of the stress field parameters is set to 90%. The optimal stress state parameters and the confidence intervals of each parameter with a confidence level of 90% under the F test are obtained.

[0028] The inversion results of the tectonic stress field are projected onto an equal area to obtain a graphical representation of the tectonic stress field in the study area, including the strike, dip, and slip angle of nodal surface 1 and the strike, dip, and slip angle of nodal surface 2 under the optimal stress state.

[0029] In step S303, the velocity field and surface dilatation rate field are obtained through GNSS data, including:

[0030] Data from GNSS stations in the study area were collected, and high-precision data processing was performed using GAMIT / GLOBK software to obtain daily relaxation solutions. These solutions were then merged with those published by IGS for overall adjustment to obtain the velocity field within the ITRF framework. This ITRF-based velocity field was then transformed into a velocity field between adjacent continents. A reference frame between adjacent continents was defined by differentiating the velocity from the ITRF with the rotation of the Eurasian Plate. During the reference frame transformation, the Euler parameters of the plates between adjacent continents obtained from the global geological model NUVELL-1A were used for Euler transformation to remove the background field of Eurasian Plate motion, thus obtaining the velocity field within the plate reference frame between adjacent continents.

[0031] The surface dilatation field of the study region is obtained from the velocity field calculation, including:

[0032] Step 1: Calculate the mean, randomness, and stability of the velocity field dataset under the plate reference frame between two adjacent continents; calculate the thresholds of the three features based on the Euler parameters and the daily relaxation solution parameters.

[0033] Step 2: Determine whether the data is surface dilation based on the threshold, remove the background field data, and retain the surface dilation data;

[0034] Step 3: Using the velocity change rate under the plate reference frame between two adjacent continents, the feature matrix in the three-dimensional space is mapped to the low-dimensional surface dilatation rate field embedding space RD, and the three feature quantities of each sample are reduced in dimensionality and fused in parameters.

[0035] Step 4: Input the fused feature data into the Deep Belief Network (DDN) and optimize the input parameters of the DDN using the fault change optimization method.

[0036] Step 5: During the training of DDN, the biases a and b of the visible and hidden layers of DDN, as well as the parameters of the number of network layers and nodes, are optimized using the fault change optimization method. After data training, mathematical models corresponding to the three types of surface dilatation are obtained, and the construction of the velocity field data self-confirmation system under the plate reference frame between two adjacent continents is completed.

[0037] Step 6: Input the signal from the velocity field data acquisition system under the plate reference frame between the two adjacent continents that needs to be diagnosed to the velocity field data self-confirmation system under the plate reference frame between the two adjacent continents. After the system processes the data, the processing result is input into the mathematical model of the velocity field data self-confirmation system under the plate reference frame between the two adjacent continents to complete the determination of the surface dilatation type of the velocity field signal acquisition system under the plate reference frame between the two adjacent continents.

[0038] In step 1, the mean, randomness, and stability of the velocity field dataset under the plate reference frame between two adjacent continents are calculated, including:

[0039] (1) Calculate the mean using the following formula:

[0040]

[0041] In the formula, F1 is the mean of the velocity field dataset, S is the total number of sample points, and u i Let there be i sample points;

[0042] (2) Calculate randomness: Arrange the collected signals from largest to smallest, merge all data of the same size, and the number of data retained is the data randomness index, denoted as F2;

[0043] (3) Calculate the stability using the following formula:

[0044]

[0045] In the formula, Let \(x(n)\) be the stable value of the current \(n\) continent sample points of the normal signal, \(u(n)\) be the stable value of the total \(n\) continent sample points of the normal signal, \(n\) be the \(n\)th continent, \(u(\eta)\) be the difference in sample points between two adjacent continents, \(\tau\) be the previous continent, \(\eta\) be a certain pair of adjacent continents, \(Q(n)\) be the positive characteristic quantity obtained by taking the square root of the stable value of the normal signal, and \(F3\) be the data stability index;

[0046] Calculate the mean threshold \(G1\), and the expression is:

[0047]

[0048] Calculate the randomness threshold \(G2\), and the expression is:

[0049]

[0050] Calculate the stability threshold \(G3\), and the expression is:

[0051]

[0052] In the formula, \(v\) i is the \(i\)th speed sample point, \(f\) s [[ID=2�]]is the data retained after arranging the collected signals in descending order of numerical value and merging all data with the same size, is the stability characteristic quantity of the normal signal.

[0053] In step 3, use the global linear embedding to calculate the rate of change of speed under the plate reference frame between two adjacent continents, including: mapping the three-dimensional space feature matrix to the low-dimensional surface dilation rate field embedding space \(O\) E and performing data dimensionality reduction and parameter fusion on the three characteristic quantities of each sample; the specific steps are:

[0054] (1) Select the local neighborhood. For a given data set, there is:

[0055] \(U = \{U1, U2 \cdots U\) N \}, U\) i \(\in O\) E

[0056] In the formula, \(U\) is the data set of the given local neighborhood, \(U\) N is the \(N\)th local neighborhood data, and \(U\) i is the \(i\)th local neighborhood data set;

[0057] Find the \(k\), \(k < N\) nearest neighbor points of each sample point \(U\) i neighborhood, and calculate using the Euclidean distance formula:

[0058]

[0059] In the formula, \(d\) ijLet u be the Euclidean distance between the i-th sample point and the j-th sample point. ik Let u be the velocity value of the i-th sample point in its neighborhood k. jk Let k be the velocity value of the j-th sample point in its neighborhood, where k is the k-th neighborhood, N is the N-th neighborhood, and D is the magnitude attenuation coefficient.

[0060] (2) Calculate the reconstruction weights of the neighborhood of the sample points, construct the local reconstruction weight matrix, and define the error function as:

[0061]

[0062] In the formula, ε(R) is the weight error, u i For i sample points, u j For j sample points, R is the second-order absolute value. ij For u i with u ij The weights between, u ij For u i The k nearest neighbors, j = 1, 2, ..., k;

[0063] As a constraint, the error function expression is transformed into:

[0064]

[0065] In the formula, R i R is the local reconstruction weight for the i-th sample point. i =[R1,R2,…,R k ] G R k Z is the local reconstruction weight of the k-th neighborhood, G is the iteration value; i For the i-th sample point, the local covariance matrix Z is obtained. i =(u i -u ij ) G (u i -u ij );

[0066] Introducing the Lagrange multipliers to solve the constrained problem, the expression is:

[0067]

[0068] In the formula, λ is the constraint weight, L is the Lagrange multiplier, and W... ij Solve for the constraint weights between the i-th sample point and the j-th sample;

[0069] Let Z i R i=1, readjust the weights so that the sum of the weights is 1, and finally obtain R. i ;

[0070] (3) Find the low-dimensional surface dilatation field embedding Z of the sample set using the obtained weight matrix R, and minimize the reconstruction error and function, expressed as:

[0071]

[0072] In the formula, φ(Z) represents the reconstruction error and the minimum value of the function, z j Let J be the local covariance matrix of j sample points;

[0073] With respect to Z, the expression is:

[0074]

[0075] In the formula, I is an N-dimensional identity matrix;

[0076] The optimization problem is transformed into a constrained optimization problem, expressed as:

[0077]

[0078] In the formula, I i Let R be the i-th identity matrix. i ZMZ represents the local reconstruction weights for the i-th sample point. G Z is the constrained optimization value for the Gth iteration. G Let Z be the dataset after G iterations;

[0079] Given a dataset Z, it is equivalent to finding the eigenvectors of a symmetric, positive semi-definite, sparse matrix M, resulting in:

[0080] M = (IW) G (IW)

[0081] In the formula, W is the constraint weight matrix for solving;

[0082] Using the Lagrange multiplier method, we get:

[0083] L(Z) = ZMZ G -λ(ZZ G -SI)

[0084] In the formula, L(Z) is the constrained optimization value obtained by the Lagrange multiplier method, and λ is the constraint weight;

[0085] Taking the partial derivative with respect to Z and minimizing L(Z) yields:

[0086]

[0087] In the formula, MZG =λZ G Finding Z is equivalent to finding the eigenvectors of M, thus obtaining MZ. G =λZ G The resulting embedded coordinates are the eigenvectors of M; the eigenvectors corresponding to the smallest d non-zero eigenvalues ​​are used as the values ​​of M, and the coordinates of the low-dimensional surface dilatation field Z are obtained. The eigenvectors corresponding to the eigenvalues ​​are the output results.

[0088] In fact, the entire process of the rate of change of velocity between two adjacent continents under the plate reference frame is shown by the following formula:

[0089]

[0090] In the formula, This is a set of velocities within a plate reference frame between two adjacent continents. Let i be the set of local reconstruction velocity weights for the i-th sample point. To embed a low-dimensional surface dilatation rate field into a set of velocity change data, → represents the rate of change.

[0091] In step 4, the fused feature data is input into the Deep Belief Network (DDN), and the input parameters of the DDN are optimized using a fault change optimization method, including:

[0092] (i) Construct a deep belief neural network (DDN);

[0093] Based on the established energy function formula, the joint probability distribution of (v,h) is calculated as follows:

[0094]

[0095] p(v,h|θ)=l -C(v,h|θ) / t(θ)

[0096] In the formula, C(v,h) is the joint probability distribution value, I is the N-dimensional identity matrix, J represents the J samples, and e i Let f be the attenuation coefficient of the i-th source, v be the velocity sample point of the i-th point, and f be the velocity sample point of the i-th source. i Let h be the randomness value of i sample points. j Let ξ be the thickness of j sample points. ji Let p(v,h|θ) be the energy coefficient between the j-th sample point and the i-th sample, and let l be the calculated value of the energy function. -C(v,h|θ) Let t(θ) be the energy value calculated under the joint probability distribution, and t(θ) be the partition function.

[0097]

[0098] In the formula, v is the velocity and h is the thickness;

[0099] The likelihood function defined by RBM, p(v|θ) is the marginal distribution of the joint probability p(v,h|θ);

[0100] Since neurons within the RBM layer are unconnected, the operating state and activation of the hidden layers are relatively independent based on the visible layer state. The activation probability of the j-th hidden layer node is:

[0101]

[0102] In the formula, F'() is the current activation probability function, h j Let j be the thickness of the hidden layers, θ be the likelihood coefficient, σ() be the sigmoid function, and f j Let j be the randomness values ​​of the sample points;

[0103] Where σ(u)=1 / (1+l) -u Let be the sigmoid function; given the states of the hidden layer nodes, the activation probability of the i-th visible layer node is:

[0104]

[0105] In the formula, e i Let be the attenuation coefficient of the i-th earthquake source;

[0106] RBMs are trained and run iteratively, and the surface dilatation field during the run is studied and calculated using the parameter θ = (ξ). ij ,e i ,f j The value of ) is obtained from the given training data and training samples; the maximum log-likelihood function is calculated using the parameter θ, with the number of training samples set to T, and the expression is:

[0107]

[0108] In the formula, θ * To calculate the maximum log-likelihood function value using parameter θ, F'(v (T) |θ) is the activation probability function for the joint probability given the number of samples T in the training set;

[0109] (ii) Training process of DDN: The velocity change rate parameter passing through the plate reference frame between two adjacent continents is fused with the velocity field signal data feature quantity under the plate reference frame between two adjacent continents as input, and an RBM is trained from bottom to top each time; in each layer, the parameter space ωk is constructed by the numbers calculated in the (k-1)th layer, and the weights are updated according to the following formula:

[0110] ξ ij =θξ ij +η( <v i h j >data - <v i h j > model )

[0111] Where ξ ij is the weight update value, θ is the update coefficient, η() is the characteristic function between two adjacent continents, <v i h j > data is the current (v, h) value, <v i h j > model is the (v, h) standard value.

[0112] In step 5, the fault change optimization method is used to optimize the parameters of the biases a and b, the number of network layers and the number of nodes of the visible layer and the hidden layer of the DDN, including:

[0113] (a) Initialization of the parameters of the fault change optimization method and the optimization problem. The optimization problem is described as:

[0114] Minimize f(u) subject to u j ∈E j

[0115] Where Minimize f(u) is the initialization value of the optimization surface expansion rate field, u j is the decision vector velocity of the jth sample, f(u) is the optimization surface expansion rate field function, j is the signal length of the jth sample, j = 1, 2, 3... N, E j is the constraint interval of u j after the decision vector;

[0116] (b) Set the initial plate and calculate the enthalpy value or entropy value. Use the unified population method to uniformly initialize the initial plate in the feasible search interval; Set two initial plates O0 = (w1, w2... w n ​​​​​​​​(d) Plate update: Plates are updated by testing the equilibrium of changes; the state when all surface expansion no longer changes is the equilibrium state; in the equilibrium state, the plate that increases the entropy of surface expansion or decreases the enthalpy value is the new plate, and abnormal plates are excluded.

[0119] (e) Determine the termination condition of the iteration: if the condition is met, the algorithm terminates; otherwise, proceed to step (c). The termination condition of the fault change optimization method is to satisfy the maximum number of iterations or the minimum enthalpy and maximum entropy.

[0120] In step S4, depth profiles are established along the NEE and NNW directions, including the vertical velocity V. p and relative velocity disturbance velocity ΔV p .

[0121] Another objective of this invention is to provide a system for constructing a three-dimensional model of a seismogenic fault based on aftershock precision location data. This system implements the method for constructing a three-dimensional model of a seismogenic fault based on aftershock precision location data. The system includes:

[0122] The information acquisition module is used to collect seismic and geological data information of the study area;

[0123] The crustal three-dimensional velocity model building module is used to build a high-resolution crustal three-dimensional velocity model of the study area and determine accurate seismic source parameters.

[0124] The tectonic stress field characteristic acquisition module is used to obtain tectonic stress field characteristics through focal mechanism solutions and GNSS observation data.

[0125] The study module predicts earthquake-generating mechanisms, seismogenic structures, and focal rupture processes in the study area. It combines velocity structure, fault distribution, focal mechanism, velocity field characteristics, and post-earthquake field survey data to predict the influence of faults on the earthquake-generating process, thereby obtaining information on potential strong earthquake hazards.

[0126] Combining all the above technical solutions, the beneficial effects of this invention are as follows:

[0127] This invention first uses a 0.2°×0.2° grid interval for double-difference seismic tomography (tomoDD) inversion to obtain the velocity structure (MODEL_0.2) of the study area and surrounding areas. Then, linear interpolation is performed on MODEL_0.2 to obtain a velocity model (MODEL_0.1) with a 0.1°×0.1° grid interval. Using MODEL_0.1 as the initial model, double-difference seismic tomography (tomoDD) is used again for inversion calculation to obtain the final velocity structure (MODEL_final) of the study area.

[0128] This invention combines two grid intervals, 0.2°×0.2° and 0.1°×0.1°, to obtain velocity structures with a 0.2°×0.2° grid interval in the non-source region and a 0.1°×0.1° grid interval in the source region. It is suitable for studying velocity structures in source regions with a relatively small number of stations.

[0129] Using the velocity model MODEL_0.1 as the initial model, imaging inversion calculations were performed with a grid interval of 0.1°×0.1°. After multiple iterations, the three-dimensional crustal P-wave velocity structure and precise seismic location results for the study area were finally obtained. The root mean square of the seismic travel time residuals was significantly reduced, and the travel time residuals before and after the inversion followed a Gaussian distribution. Attached Figure Description

[0130] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this disclosure and, together with the description, serve to explain the principles of this disclosure;

[0131] Figure 1 This is a flowchart of the method for constructing a three-dimensional model of a seismogenic fault based on aftershock precision location data provided in an embodiment of the present invention;

[0132] Figure 2 The following are the test results diagrams of the test plate with a grid interval of 0.2° provided for the prior art; wherein, (a) is the test result diagram of the test plate with Z=0km, (b) is the test result diagram of the test plate with Z=5km, (c) is the test result diagram of the test plate with Z=9km, (d) is the test result diagram of the test plate with Z=14km, (e) is the test result diagram of the test plate with Z=20km, and (f) is the test result diagram of the test plate with Z=25km;

[0133] Figure 3 The following are the test results diagrams of the test plate with a grid interval of 0.1° provided for the prior art; wherein, (a) is the test result diagram of the test plate with Z=5km, (b) is the test result diagram of the test plate with Z=9km, (c) is the test result diagram of the test plate with Z=14km, and (d) is the test result diagram of the test plate with Z=20km.

[0134] Figure 4 A flowchart of velocity structure inversion provided for this invention;

[0135] Figure 5 Histogram of residual distribution of travel time before (hollow rod) and after (solid rod) positioning provided by the present invention;

[0136] Figure 6The present invention provides velocity distribution maps for six depth levels ranging from 5 to 31 km; wherein, (a) is the velocity distribution map for the depth level of Z=5 km, (b) is the velocity distribution map for the depth level of Z=9 km, (c) is the velocity distribution map for the depth level of Z=14 km, (d) is the velocity distribution map for the depth level of Z=20 km, (e) is the velocity distribution map for the depth level of Z=25 km, and (f) is the velocity distribution map for the depth level of Z=31 km.

[0137] Figure 7 The azimuth distribution diagrams of the P-axis and T-axis in areas A and B provided by this invention;

[0138] Figure 8 This is a diagram showing the stress field inversion results for a certain city area according to the present invention;

[0139] Figure 9 The velocity V on the NEE profile of this invention p and relative velocity disturbance ΔV p Distribution diagram; where (a1) shows the velocity V on the NEE profile. p The A-AA cross-section diagram, and (a2) shows the relative velocity disturbance ΔV. p The A-AA cross-section diagram, (b1) shows the velocity V on the NEE profile. p The B-BB profile, (b2) shows the relative velocity disturbance ΔV. p The B-BB cross-section diagram, (c1) shows the velocity V on the NEE profile. p The C-CC profile, (b2) shows the relative velocity disturbance ΔV. p The C-CC profile diagram, (d1) shows the velocity V on the NEE profile. p The D-DD profile, (d2) shows the relative velocity disturbance ΔV. p D-DD cross-section;

[0140] Figure 10 The velocity V on the NNW-direction cross-section of this invention p and relative velocity disturbance ΔV p Distribution diagram; where (a1) shows the velocity V on the NEE profile. p The E-EE profile, (a2) shows the relative velocity disturbance ΔV. p The E-EE profile diagram, (b1) shows the velocity V on the NEE profile. p The F-FF profile, (a2) shows the relative velocity disturbance ΔV p The F-FF cross-section diagram, (c1) shows the velocity V on the NEE profile. p The G-GG profile, (c2) shows the relative velocity disturbance ΔV. pThe G-GG cross-section diagram, (d1) shows the velocity V on the NEE profile. p The H-HH profile, (d2) shows the relative velocity disturbance ΔV. p H-HH cross-section diagram. Detailed Implementation

[0141] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0142] Example 1, such as Figure 1 As shown, the method for constructing a three-dimensional model of the seismogenic fault based on aftershock precision location data includes:

[0143] S1, collecting earthquake and geological data information of the study area;

[0144] S2, establish a high-resolution three-dimensional crustal velocity model of the study area and determine the precise source parameters;

[0145] S3, the characteristics of the tectonic stress field are obtained through focal mechanism solutions and GNSS observation data;

[0146] S4 predicts the earthquake-generating mechanism, seismogenic structure, and source rupture process in the study area. By combining velocity structure, fault distribution, source mechanism, velocity field characteristics, and post-earthquake field investigation data, it predicts the influence of faults on the earthquake-generating process, thereby obtaining information on potential strong earthquake hazards.

[0147] For example, the present invention also provides a system for constructing a three-dimensional model of a seismogenic fault based on aftershock precision location data, comprising:

[0148] The information acquisition module is used to collect seismic and geological data information of the study area;

[0149] The crustal three-dimensional velocity model building module is used to build a high-resolution crustal three-dimensional velocity model of the study area and determine accurate seismic source parameters.

[0150] The tectonic stress field characteristic acquisition module is used to obtain tectonic stress field characteristics through focal mechanism solutions and GNSS observation data.

[0151] The study module predicts earthquake-generating mechanisms, seismogenic structures, and focal rupture processes in the study area. It combines velocity structure, fault distribution, focal mechanism, velocity field characteristics, and post-earthquake field survey data to predict the influence of faults on the earthquake-generating process, thereby obtaining information on potential strong earthquake hazards.

[0152] This invention establishes a high-resolution three-dimensional velocity model of the Earth's crust and obtains accurate source parameter information; combined with the source mechanism solutions of small and medium-sized earthquakes and GNSS observation data, it determines the tectonic stress field of the study area; based on the coupling mode between the three-dimensional velocity structure of the crust, earthquake distribution and tectonic stress field, and combined with field geological survey data, it studies the earthquake-generating mechanism, seismogenic structure and source rupture process, and judges the potential strong earthquake hazard.

[0153] This invention obtained the three-dimensional crustal velocity structure with a resolution of 0.1°×0.1° in the source area and 0.2°×0.2° in the non-source area of ​​a certain city, as well as the source parameters of 6255 earthquakes. Based on the focal mechanism solutions of 64 earthquakes in the city and GNSS data from 7 stations in the city, the direction of the maximum principal compressive stress in the main tectonic stress field of the city was determined to be NW-SE. Based on the information of the three-dimensional crustal velocity structure and source parameters, as well as the coupling mode between the tectonic stress field, and combined with field geological survey data, the relationship between the crustal velocity structure characteristics and seismic geology and faults in the city was analyzed, and four major earthquakes in the history of the city were identified. L The seismogenic structures and earthquake-generating mechanisms of earthquakes of magnitude 5.0 and above provide important evidence for assessing the seismic hazard of the region.

[0154] Example 2, as another implementation of the present invention, further describes the method for constructing a three-dimensional model of the seismogenic fault based on aftershock precision location data, based on the theoretical basis of the present invention.

[0155] The characteristics of the tectonic stress field were obtained through focal mechanism solutions and GNSS observation data: The focal mechanism solutions for small and medium-sized earthquakes in a certain city were obtained by combining the P-wave initial motion with the amplitude ratio of P-waves and S-waves. Based on continuous GNSS observation data, the regional motion field and surface dilatation rate field were calculated through fitting. Combined with the focal mechanism solutions, the characteristics of the tectonic stress field in the city were studied to understand the stress accumulation and crustal movement within the region. Taking into account the accurate focal parameter information after relocation, and based on the principle of small earthquake clustering, mathematical methods were used to re-invert the parameters of the fault plane that caused the earthquake in the city. On this basis, the slip angle of the fault plane was calculated with reference to the tectonic stress field.

[0156] A Study on Earthquake-Pregnancy Mechanism, Seismogenic Structure, and Focal Rupture Process in a Certain City: A field geological survey was conducted in the earthquake-prone area of ​​the city, investigating the geomorphological features, lithological distribution, planar distribution of exposed strata, location and occurrence of surface fault zones, and other geological conditions. Combining the three-dimensional crustal wave velocity fine structure, earthquake distribution, tectonic stress field, and field geological survey data, the location of the seismogenic structures in the city was identified. The relationship between velocity anomalies, seismic activity, and tectonic stress field was analyzed, and the earthquake-pregnancy mechanism, seismogenic structure, and focal rupture process were studied. Comparative analysis was conducted with representative earthquakes from the Southeast Coastal Seismic Belt and the eastern edge of the Tibetan Plateau, focusing on crustal velocity structure, seismic activity, fault structure characteristics, and crustal movement. This improved the understanding of the earthquake-pregnancy mechanism, seismogenic structure, and focal rupture process in the study area, and ultimately analyzed the potential strong earthquake hazard in the region.

[0157] Example 3, as another embodiment of the present invention, mainly includes the following research progress, important results, key data, and their scientific significance or application prospects.

[0158] (2.1) Detailed Seismic Geological Conditions of a Certain City. A field geological survey was conducted in the city, and previous data on the seismic geology of the area were collected, ultimately yielding relatively detailed seismic geological data for the city. The city is located at the southwestern end of the Shaoguang Depression in the Zhejiang-Guangdong Dome System, situated at the southern end of the western edge of the second giant uplift zone of the Neo-Cathaysian tectonic system. The seismic geological structure is relatively complex, with four main faults: the Yangchun-Zhizhen Fault (F1), the eastern branch of the Wuchuan-Sihui Fault; the southwestern segment of the Cangcheng-Hailing Fault (F2), the western branch of the Enping-Xinfeng Fault; the Yangbianhai Fault (F3); and the Pinggang Fault (F4). The Yangchun-Jizhen Fault (F1) generally strikes NE or NNE, starting near Xinxing in the north, extending southwest through Yangchun and Yangxi to the Shaba area before entering the sea. It is inferred to continue southwestward, with an onshore length of 150 km. It primarily dips NW, with local dips SE, at angles of 60-80°. This fault originated in the Indosinian period and exhibited thrust characteristics, reactivating and transforming into a normal fault during the Yanshanian orogeny. The southwestern segment of the Cangcheng-Hailing Fault (F2) generally strikes NE, starting southwest of Enping in the northeast, passing through Dahai and Nalong to the southeast of a certain city, crossing Hailing Island and extending into the South China Sea. It is approximately 40 km long, dips NW, and at angles exceeding 60°. The Cangcheng-Hailing Fault is largely composed of high-angle reverse faults, cutting through the Lower Paleozoic, Upper Paleozoic, Yanshanian granites, Cretaceous, and Paleogene strata, controlling the sedimentary boundaries of Cretaceous and Paleogene rift basins. The Yangbianhai Fault, also known as the Fengtouhe Fault (F3), exposes on the northeastern beach of Yangbianhai, trending NW and dipping NE or SW at an angle of 65-75°. It is approximately 25 km long and can be divided into northern and southern segments. The northern segment mainly extends along Yangbianhai, with most of it submerged underwater. It is a left-lateral normal fault, only exposed above water on the right bank of Yangbianhai. The southern segment is found in the villages of Shajiao, Beiting, and Guangling on Hailing Island. The fault slip surface is clear, and its striations indicate a left-lateral strike-slip nature. The Pinggang Fault (F4) trends NE to NEE, starting from Beiguan southwest and extending southwest into Yangbianhai. It is approximately 25 km long and is a typical right-lateral normal fault of the South China Sea system. The northeastern segment dips SE at an angle of 60-70°, while the southwestern segment has a steeper dip, even becoming vertical.

[0159] (2.2) Establishment of a high-resolution three-dimensional velocity model of the crust in a certain city and determination of precise source parameters.

[0160] To improve the resolution of the three-dimensional velocity structure in a certain city, this invention utilizes as much reliable seismic data as possible for inversion. P-wave arrival data from January 1990 to August 2019, recorded by a provincial seismic network, were collected. The accuracy of simulated seismic records before June 2007, especially before 2000, was relatively low. Therefore, a double-difference seismic location method was used to relocate earthquakes before June 2007, resulting in relocated phase reports. Further screening of the above seismic data removed phase data deviating from the theoretical travel time curve by more than 2 seconds, ensuring that each earthquake was recorded by at least 3 stations. Ultimately, 43,225 absolute P-wave arrival data and 422,956 relative arrival data from 6,390 earthquakes (736 before June 2007) recorded by 49 stations were used for inversion calculations. The focal depth ranged from 0 to 32 km. The magnitude range of earthquakes before and after June 2007 was M... L 0.9-4.4 and M L 0.2-4.2.

[0161] To determine the optimal resolution of underground velocity structures using existing seismic data, this invention uses grid intervals of 0.2°×0.2° and 0.1°×0.1° to perform a test panel on a certain city area. The vertical grid nodes are set at 0km, 5km, 9km, 14km, 20km, 25km, 31km and 40km, respectively. The initial model used is a one-dimensional P-wave velocity model of the city and its neighboring areas (MODEL_1d, see Table 1).

[0162] Table 1. One-dimensional P-wave velocity structure model of a certain city and its neighboring areas.

[0163] Depth (km) 0 5 9 14 20 25 31 40 Vp(km / s) 4.97 5.70 6.00 6.04 6.15 6.31 6.72 6.96

[0164] The test results show that, using a grid interval of 0.2° × 0.2°, the restoration effects on land and sea areas in a certain city are significantly different, with the land area showing better restoration at depths of 5-31 km. Figure 2 The test results for the detection plate with a 0.2° grid interval are shown. In the marine area, the recovery effect is better only in areas close to land. However, using a 0.1°×0.1° grid interval, the recovery effect is better in the land area at depths of 9km and 14km. Figure 3The results of the test panel with a 0.1° grid interval are shown. At depths of 5km and 20km, only the vicinity of the epicenter area (black box) in a certain city showed good reconstruction. This indicates that using a 0.2°×0.2° grid interval provides good resolution across all depth levels in the study area. However, given that the epicenter area of ​​the certain city is approximately 0.2°×0.2°, a grid interval of only 0.2°×0.2° is insufficient to accurately reflect the epicenter area. Using a 0.1°×0.1° grid interval, the epicenter area (black box) in the certain city achieves good resolution, but the resolution outside the epicenter area is not ideal.

[0165] To address the aforementioned issues, this invention first employs double-difference seismic tomography (tomoDD) inversion with a 0.2°×0.2° grid interval to obtain the velocity structure (MODEL_0.2) of a city and its surrounding areas. Linear interpolation is then performed on MODEL_0.2 to obtain a velocity model (MODEL_0.1) with a 0.1°×0.1° grid interval. Using MODEL_0.1 as the initial model, double-difference seismic tomography (tomoDD) is used again for inversion calculations to obtain the final velocity structure (MODEL_final) of the city area. The above scheme ( Figure 4 The velocity structure inversion flowchart shown combines two grid intervals, 0.2°×0.2° and 0.1°×0.1°, to obtain both the velocity structure of the non-source region with a 0.2°×0.2° grid interval and the velocity structure of the source region with a 0.1°×0.1° grid interval. This is suitable for the study of velocity structure in the source region with a relatively small number of stations.

[0166] Using the velocity model MODEL_0.1 as the initial model, imaging inversion calculations were performed with a grid interval of 0.1°×0.1°. After 20 iterations, the three-dimensional crustal P-wave velocity structure and the precise location results of 6,255 earthquakes in a certain city area were finally obtained. The root mean square of the earthquake travel time residuals decreased from 365ms to 0.1ms, and the travel time residuals before and after the inversion followed a Gaussian distribution. Figure 5 The histograms showing the distribution of travel time residuals before (hollow rod) and after (solid rod) positioning are shown. Before inversion, the travel time residuals are located in the range of -0.6s to 0.6s, while after inversion, the residuals are significantly reduced and mainly located in the range of -0.1s to 0.1s. Among them, the amount of data with residuals within 0.05s reaches 93%.

[0167] Horizontal velocity distribution: Through inversion calculations using double-difference seismic tomography, this invention ultimately yields a high-precision three-dimensional P-wave velocity model (MODEL_final) for a certain city area. Figure 6 The map shows the velocity distribution at six depth levels, ranging from 5 to 31 km. Figure 6 In the middle (a) diagram, the five-pointed star represents a magnitude 6.4 earthquake that occurred in a certain city in 1969; Figure 6 In the middle (b) diagram, the five-pointed stars represent the two magnitude 5.0 earthquakes that occurred in 1986 and 1987; Figure 6 In the middle (c) diagram, the pentagram represents the 4.9 magnitude earthquake in 2004; F1: Yangchun-Zhizhen Fault; F2: Southwest segment of Cangcheng-Hailing Fault; F3: Yangbianhai Fault; F4: Pinggang Fault;

[0168] It can be seen that the velocity structure of a certain city and its neighboring areas exhibits significant lateral heterogeneity. The velocity distribution at a depth of 5 km shows a clear correspondence with fault structures and geological formations. The Mashui-Hekou-Tangkou-Pupai area is located on a high-velocity / low-velocity anomaly transition zone. This transition zone north of Tangkou follows the trend of the Yangchun-Zhizhen fault, while south of Tangkou, the fault deviates towards the high-velocity side. The Sanjia-Xinwei-Jueshan area west of the transition zone shows a distinct low-velocity anomaly, corresponding to the Yanshanian granite system in the region, extending to a depth of 20 km. The Xinwei-Jueshan area then gradually transitions to a high-velocity anomaly at a depth of 25 km. The Yangchun-Yangxi area east of the transition zone shows a NE-trending high-velocity anomaly, corresponding to the Cambrian metamorphic rock system. The southern segment of the Yangchun-Zhizhen fault is located within this high-velocity anomaly. Two earthquake swarms appear on either side of the fault zone; the western swarm is located on the high-velocity / low-velocity transition zone, while the eastern Xinhu Reservoir swarm is located within the high-velocity body. The eastern side of the NE-trending high-velocity anomaly exhibits a smaller-amplitude low-velocity anomaly, corresponding to the Yanshanian and Caledonian granitic systems. However, a small area of ​​Cambrian metamorphic rocks is distributed near Chengcun, thus the low-velocity zone near Chengcun is truncated by the high-velocity anomaly located at the Yangbianhai Fault. The Pinggang Fault is situated at the boundary between the high and low-velocity anomalies, biased towards the high-velocity anomaly. Earthquakes are relatively concentrated at the intersection of the Yangbianhai Fault and the Pinggang Fault; the 6.4 magnitude earthquake in 1969 occurred here.

[0169] 9km depth ( Figure 6 (Figure b) shows that the range of high-velocity anomalies increases. The Gushan-Hekou-Tangkou-Pupai area is a transition zone between high and low velocity anomalies, with low velocity anomalies to the west and high velocity anomalies to the east. The Hailing-Yashao-Beiguan area is a low velocity anomaly. The NEE-trending seismic strips of Yangbianhai-Pinggang... Figure 6 (b) The yellow dashed line in the figure intersects with the Pinggang Fault and the Yangbianhai Fault. The strip shows the boundary between high and low velocity anomalies. The northwest is a high velocity anomaly and the southeast is a low velocity anomaly. The two magnitude 5.0 earthquakes in 1986 and 1987 both occurred at the intersection of the NEE-trending seismic strip and the Yangbianhai Fault. The magnitude 5.0 earthquake in 1986 was located at the junction of the north and south segments of the Yangbianhai Fault, while the magnitude 5.0 earthquake in 1987 was located in the north segment of the Yangbianhai Fault.

[0170] 14km ​​depth ( Figure 6(Figure (c)) shows that the Yangbianhai Fault is located at the boundary between high and low velocity anomalies. Near Pinggang, it transitions to a low velocity anomaly. With increasing depth, the low velocity anomaly near Pinggang expands and merges with the low velocity anomaly in the Sanjia-Xinwei-Jueshan area at a depth of 20 km. (20 km depth) Figure 6 (Figure d) shows that the overall velocity structure exhibits a low-velocity anomaly, the NEE-trending seismic zonation of the Yangbianhai-Pinggang fault has disappeared, and the seismic swarm activity on both sides of the Yangchun-Zhizhen fault has ceased. (25km to 31km depth) Figure 6 In Figures (e) and (f), the 6.5 km / s contour lines gradually become parallel to the coastline.

[0171] (2.3) Obtain the characteristics of tectonic stress field through focal mechanism solutions and GNSS observation data.

[0172] ① Focal mechanism solutions. This invention obtained focal mechanism solutions for 64 earthquakes, including 9 earthquakes prior to 1999 (M... L The focal mechanism solutions (≥3.6) given in the analysis of typical earthquakes along the southeast coast (upper hemisphere projection) used simulated seismic network data. The focal mechanism solutions for these earthquakes were recalculated to obtain the corresponding lower hemisphere projection parameters. Seven earthquakes (M ≥3.6) from 1999 to 2003... L The results of the snoke method were obtained when performing characteristic analysis on focal mechanism solutions for a certain province with focal mechanism values ​​≥3.0. Six focal mechanism solutions were obtained from 2004 to 2006, which were obtained during stress field inversion of the province and its adjacent areas. However, for the 42 focal mechanism solutions since 2007... L For the focal mechanism solutions of earthquakes of magnitude 3.6 and above, this invention uses an interactive program (Focmec-Interface) based on the Focmec method to invert the focal mechanism solution.

[0173] According to the distribution map of focal mechanism solutions in a certain city, the ML3.0 earthquake in that city was mainly distributed in two areas: the intersection of the Yangbianhai Fault and the Pinggang Fault (Area A, 46 focal mechanisms) and the area near Xinhu Reservoir on the east side of the Yangchun-Zhizhen Fault (Area B, 18 focal mechanisms).

[0174] The focal mechanism solutions in Region A are relatively complex, with right-lateral slip predominating (46%). Ten earthquakes in this region also exhibit normal slip characteristics, and two exhibit reverse slip characteristics. Additionally, eleven earthquakes (26%) show reverse slip characteristics. The principal compressive stress axis (P-axis) and principal stress axis (T-axis) do not show obvious convergence properties. Figure 7(This is a map showing the azimuth distribution of the P-axis and T-axis in areas A and B). All four earthquakes of magnitude 4.9 or higher were located in area A, but the nature of the fault slippage differed significantly. The largest earthquake was a magnitude 6.4 (ML6.6) earthquake on July 26, 1969, with a strike-slip focal mechanism. The two nodal planes trended N74°E and N20°W, respectively, consistent with the trends of the Yangbianhai Fault and the Pinggang Fault. The focal mechanism solutions of the magnitude 5.0 earthquakes of January 28, 1986, and February 25, 1987, are quite similar, both being strike-slip earthquakes. The 1987 magnitude 5.0 earthquake also exhibited some normal fault characteristics. The two nodal planes of both earthquakes trend NE and NW, respectively. The post-earthquake isoseismal lines of both earthquakes are elliptical with a major axis trending NE. It is believed that the seismogenic structure of the two strong aftershocks is still the Pinggang Fault trending NEE. However, if the seismogenic structure of the two earthquakes is indeed the Pinggang Fault, then its left-lateral strike-slip nature contradicts the right-lateral strike-slip normal fault nature of the Pinggang Fault. Seventeen years later, on September 17, 2004, a magnitude 4.9 earthquake occurred near Pinggang, northeast of the aforementioned three earthquakes. This earthquake was a reverse fault, and the focal mechanism solution shows that the nodal plane trending NEE, consistent with the direction of the major axis of the isoseismal lines and also consistent with the trend of the Pinggang Fault.

[0175] In summary, the focal mechanism of the magnitude 6.4 mainshock in a certain city was strike-slip, with the two nodal planes trending NNW and NEE respectively, consistent with the strike of the Yangbianhai Fault and the Pinggang Fault. The focal mechanism solutions of the three strong aftershocks of magnitude 4.9 and above all trended NE-NEE, consistent with the strike of the Pinggang Fault, which is the main controlling fault of the earthquake in Area A. However, the slip modes revealed by the focal mechanisms of the four earthquakes are quite different, and further analysis is needed in conjunction with the characteristics of underground velocity structure and tectonic stress field.

[0176] Most of the 18 earthquakes in Zone B occurred on the west side of the Yangchun-Zhizhen Fault, especially near the Xinhu Reservoir. The focal mechanism solutions in this zone showed good consistency, with 94% of the earthquakes exhibiting right-lateral strike-slip or normal-dipping characteristics. The T-axis azimuth and dip angles were relatively stable, with the T-axis azimuth mainly between NE16-40° and the dip angle mainly between 0-30°. The P-axis azimuth mainly ranged between 140±50° and 208±20°, but the dip angle varied considerably. Figure 7 ).

[0177] ② Solve the tectonic stress field of a certain city based on the focal mechanism.

[0178] The grid search method proposed by Wan Yongge is used to solve the tectonic stress field of a certain city. Since the number and quality of seismic receiving stations vary for different magnitudes, the accuracy of the focal mechanism solutions differs. This invention assumes that the larger the earthquake magnitude, the higher the accuracy of the focal mechanism solution. Therefore, the focal mechanism solutions are weighted based on magnitude to reasonably combine the focal mechanisms of large and small earthquakes. The weights given in this invention are given by the following formula:

[0179]

[0180] Set the weight ω of the earthquake with the smallest magnitude in ML3.0. min The weight ω of the largest earthquake, ML6.6, is 1.0. max The magnitude is 5, where r is the relative magnitude, r = mm. min ;

[0181] The magnitude attenuation coefficient D is derived from the maximum and minimum regional magnitudes, and its expression is:

[0182] D = (<-M) min ) / ln(ω max )

[0183] In the formula, ω is the given weight, e r A given coefficient for source attenuation, D is the magnitude attenuation coefficient, ω max The maximum magnitude weight is M, where M is the current magnitude source. min The smallest magnitude earthquake source;

[0184] The D value varies in different regions, and the D value in the study area is 2.2368. In the actual inversion, the search interval of the three rotation angles of the stress field rotation axis is set to 1.0°, the search interval of the stress ratio is set to 0.1, and the confidence level of the stress field parameters is set to 90%. The optimal stress state parameters and the confidence intervals of each parameter with a confidence level of 90% under the F test are obtained (as shown in Table 2).

[0185] Table 2. Results of Tectonic Stress Field Inversion in a Certain City

[0186]

[0187] In Table 2, the numerical range in parentheses represents the uncertainty range of each parameter under the F-test at a 90% confidence level, and the values ​​above them represent the optimal stress state parameters.

[0188] The inversion results of the tectonic stress field are projected onto an equal area to obtain a graphical representation of the tectonic stress field in a certain city. Figure 8 Stress field inversion results for a certain city area. Figure 8 Figure (a) and Figure 8Figure (b) shows the equal-area projection of the stress field in regions A and B, respectively. The black arcs represent focal mechanism solution surfaces, and the green arcs represent the distribution of focal mechanism solution surfaces corresponding to the stress field with a 90% confidence level. Yellow arrows indicate the slip direction on the optimal surface, small red arrows indicate the theoretical slip direction, and small blue arrows indicate the actual slip direction. Large red arrows indicate the optimal direction of the S1 axis, and large blue arrows indicate the optimal direction of the S3 axis. The closed curves at points P, B, and T represent the range of principal stress parameters S1, S2, and S3 with a 90% confidence level. Under the optimal stress state, the strike, dip, and slip angle of surface 1 are 268.9°, 75°, and 179.2°, respectively, and the strike, dip, and slip angle of surface 2 are 359.1°, 89.2°, and 15.0°, respectively. It can be seen that the direction of the maximum principal compressive stress S1 axis in a certain city is N47°W, and the direction of the maximum principal compressive stress axis S3 axis is N45°E, indicating that the tectonic stress field is dominated by strike-slip displacement.

[0189] Analysis suggests that the area has experienced multiple tectonic movements, such as the Caledonian, Indosinian, and Yanshanian orogenies. These movements generated strong NW-SE trending compression in the region, and the residual compressive stresses still exert some influence on the current tectonic stress field. Geological surveys indicate that the Yangchun-Jizhen fault generally strikes NE or NNE, predominantly dipping NW, with some areas dipping SE, and is a high-angle thrust fault. The Cangcheng-Hailing fault strikes NE, dips NW, and has a dip angle exceeding 60°; it is also largely composed of high-angle thrust faults. This indicates that the area is influenced by NW-SE trending compression, consistent with the direction of the maximum principal compressive stress in the NW-SE trending tectonic stress field.

[0190] The Yangbianhai Fault strikes NW and dips NE, and is divided into two segments, north and south. The northern segment mainly extends along the Yangbianhai (Fengtou River) and is a left-lateral normal fault. The striations on the fault slip surface show that it is a left-lateral strike-slip fault. This indicates that the Yangbianhai Fault is affected by NE-SW extensional stress, which is consistent with the direction of the maximum assertive stress in the NE-SW tectonic stress field.

[0191] Based on the tectonic stress field, it can be inferred that the right-lateral strike-slip nature of the NE-trending Pinggang Fault and the left-lateral strike-slip nature of the NW-trending Yangbianhai Fault are a combination of NW-SE-trending compression and NE-SW-trending extensional tectonic stresses. However, its normal fault characteristics are not consistent with the NW-SE-trending compressional tectonic stress in this region, requiring further analysis.

[0192] ③ Obtaining the velocity field and surface dilatation rate field from GNSS data. This invention collected data from seven stations in a certain city from 2018 to 2019. High-precision data processing was performed using GAMIT / GLOBK software to obtain daily relaxation solutions. These solutions were then combined with daily relaxation solution files published by IGS for overall adjustment, finally obtaining the velocity field under the ITRF framework. The velocity field under the ITRF framework was converted into a velocity field relative to the Eurasian Plate. A Eurasian reference frame was defined by differentiating the velocity in the ITRF with the rotation of the Eurasian Plate. In this invention, the Euler parameters of the Eurasian Plate obtained from the global geological model NUVELL-1A were used for Euler transformation during the reference frame conversion to remove the background field of Eurasian Plate motion, thus obtaining the velocity field relative to the Eurasian Plate reference frame.

[0193] Based on the velocity field of a certain city, the recent crustal movement characteristics of the study area can be obtained. The movement direction of all seven stations is east-southeast, with a smaller southward component, consistent with the NW-SE direction of the maximum principal compressive stress obtained from the focal mechanism solution. Analysis suggests that the NW-SE direction of the maximum principal compressive stress in the city is consistent with the overall velocity field direction of South China, mainly due to the combined effects of the SE-direction lateral pressure derived from the strong uplift of a plateau under the compression of the Indian Plate and the NW-SE-direction compression generated by the subduction of a certain sea plate beneath the Eurasian Plate. The overall movement velocity of the city ranges from 1.02 to 11.26 mm / a. The relatively larger movement velocities at the northern stations YCSJ, YJJZ, and YDDG may indicate localized locking within the area, which is conducive to strain accumulation. The areal dilatation field of the study area, calculated from the velocity field, shows a relatively small overall areal dilatation value. The areal dilatation field indicates the existence of a near-NS-direction contraction-expansion transition zone in the city, which is conducive to stress accumulation.

[0194] For example, the surface dilatation rate field of the study area calculated from the velocity field includes:

[0195] Step 1: Calculate the mean, randomness, and stability of the velocity field dataset under the plate reference frame between two adjacent continents; calculate the thresholds of the three features based on the Euler parameters and the daily relaxation solution parameters.

[0196] Step 2: Determine whether the data is surface dilation based on the threshold, remove the background field data, and keep the surface dilation data for subsequent analysis;

[0197] Step 3: Map the feature matrix in three-dimensional space to a low-dimensional surface dilatation field embedding space R using the velocity change rate under a global linear embedding framework between two adjacent continental plates. D The above steps involve dimensionality reduction and parameter fusion of the three features of each sample.

[0198] Step 4: Input the fused feature data into the Deep Belief Network (DDN) and optimize the input parameters of the DDN using the fault change optimization method.

[0199] Step 5: During the training of DDN, the biases a and b of the visible and hidden layers of DDN, as well as the parameters of the number of network layers and nodes, are optimized using the fault change optimization method. After data training, mathematical models corresponding to the three types of surface dilatation are obtained, and the construction of the velocity field data self-confirmation system under the plate reference frame between two adjacent continents is completed.

[0200] Step 6: Input the signal from the velocity field data acquisition system under the plate reference frame between the two adjacent continents that needs to be diagnosed to the velocity field data self-confirmation system under the plate reference frame between the two adjacent continents. After the system processes the data, the processing result is input into the mathematical model of the velocity field data self-confirmation system under the plate reference frame between the two adjacent continents to complete the determination of the surface dilatation type of the velocity field signal acquisition system under the plate reference frame between the two adjacent continents.

[0201] The three characteristic quantities of mean, randomness, and stability of the velocity field dataset under the plate reference frame between each pair of adjacent continents in step 1 include:

[0202] (1) Calculate the mean using the following formula:

[0203]

[0204] In the formula, F1 is the mean of the velocity field dataset, S is the total number of sample points, and u i Let there be i sample points;

[0205] (2) Calculate randomness: Arrange the collected signals from largest to smallest, merge all data of the same size, and the number of data retained is the data randomness index, denoted as F2;

[0206] (3) Calculate the stability using the following formula:

[0207]

[0208] In the formula, Let u(n) be the stable value of the current n continent sample points of the normal signal, n be the nth continent, u(η) be the sample point difference between two adjacent continents, τ be the previous continent, η be a certain two adjacent continents, Q(n) be the positive characteristic quantity of the square root of the stable value of the normal signal, and F3 be the data stability index.

[0209] The mean threshold G1 is calculated using the following expression:

[0210]

[0211] Calculate the randomness threshold G2, and the expression is:

[0212]

[0213] Calculate the stability threshold G3, and the expression is:

[0214]

[0215] In the formula, v i is the i-th speed sample point, f s is the signal collected, arranged from large to small in value, and the data retained after combining all data of the same size, is the stability characteristic quantity of the normal signal.

[0216] In step 3, use the rate of change of velocity under the plate reference frame between adjacent continents by global linear embedding to map the three-dimensional space feature matrix to the low-dimensional surface expansion rate field embedding space O E On it, data dimensionality reduction and parameter fusion of the three characteristic quantities of each sample include:

[0217] (1) Select the local neighborhood. For a given data set, there is:

[0218] U = {U1, U2…U N}, U i ∈O E

[0219] In the formula, U is the data set of the given local neighborhood, U N is the N-th local neighborhood data, U i is the i-th local neighborhood data set;

[0220] Find the k, k < N nearest neighbor points of each sample point U i neighborhood, and calculate using the Euclidean distance formula:

[0221]

[0222] In the formula, d ij is the Euclidean distance between the i-th sample point and the j-th sample point, u ik is the velocity value of the i-th sample point at the k in the neighborhood, u jk is the velocity value of the j-th sample point at the k in the neighborhood, k is the k-th neighborhood, N is the N neighborhoods, and D is the magnitude attenuation coefficient;

[0223] (2) Calculate the reconstruction weight of the sample point neighborhood, construct the local reconstruction weight matrix, and define the error function as:

[0224]

[0225] In the formula, ε(R) is the weight error, u i For i sample points, u j For j sample points, R is the second-order absolute value. ij For u i with u ij The weights between, u ij For u i The k nearest neighbors, j = 1, 2, ..., k;

[0226] As a constraint, the error function expression is transformed into:

[0227]

[0228] In the formula, R i R is the local reconstruction weight for the i-th sample point. i =[R1,R2,…,R k ] G R k Z is the local reconstruction weight of the k-th neighborhood, G is the iteration value; i For the i-th sample point, the local covariance matrix Z is obtained. i =(u i -u ij ) G (u i -u ij );

[0229] Introducing the Lagrange multipliers to solve the constrained problem, the expression is:

[0230]

[0231] In the formula, λ is the constraint weight, L is the Lagrange multiplier, and W... ij Solve for the constraint weights between the i-th sample point and the j-th sample;

[0232] Let Z i R i =1, readjust the weights so that the sum of the weights is 1, and finally obtain R. i ;

[0233] (3) Find the low-dimensional surface dilatation field embedding Z of the sample set using the obtained weight matrix R, and minimize the reconstruction error and function, expressed as:

[0234]

[0235] In the formula, φ(Z) represents the reconstruction error and the minimum value of the function, z jLet J be the local covariance matrix of j sample points;

[0236] With respect to Z, the expression is:

[0237]

[0238] In the formula, I is an N-dimensional identity matrix; the optimization problem is transformed into a constrained optimization problem, expressed as:

[0239]

[0240] In the formula, I i Let R be the i-th identity matrix. i ZMZ represents the local reconstruction weights for the i-th sample point. G Z is the constrained optimization value for the Gth iteration. G Let Z be the dataset after G iterations;

[0241] Given a dataset Z, it is equivalent to finding the eigenvectors of a symmetric, positive semi-definite, sparse matrix M, resulting in:

[0242] M = (IW) G (IW)

[0243] In the formula, W is the constraint weight matrix for solving;

[0244] Using the Lagrange multiplier method, we get:

[0245] L(Z) = ZMZ G -λ(ZZ G -SI)

[0246] In the formula, L(Z) is the constrained optimization value obtained by the Lagrange multiplier method, and λ is the constraint weight;

[0247] Taking the partial derivative with respect to Z and minimizing L(Z) yields:

[0248]

[0249] In the formula, MZ G =λZ G Finding Z is equivalent to finding the eigenvectors of M, thus obtaining MZ. G =λZ G The resulting embedded coordinates are the eigenvectors of M; the eigenvectors corresponding to the smallest d non-zero eigenvalues ​​are used as the values ​​of M, and the coordinates of the low-dimensional surface dilatation field Z are obtained. The eigenvectors corresponding to the eigenvalues ​​are the output results.

[0250] In fact, the entire process of the rate of change of velocity between two adjacent continents under the plate reference frame is shown by the following formula:

[0251]

[0252] In the formula, This is a set of velocities within a plate reference frame between two adjacent continents. Let i be the set of local reconstruction velocity weights for the i-th sample point. To embed a low-dimensional surface dilatation rate field into a set of velocity change data, → represents the rate of change.

[0253] In step 4, the fused feature data is input into the Deep Belief Network (DDN), and the input parameters of the DDN are optimized using a fault change optimization method, including:

[0254] (i) Construct a deep belief neural network (DDN);

[0255] Based on the established energy function formula, the joint probability distribution of (v,h) is calculated as follows:

[0256]

[0257] p(v,h|θ)=l -C(v),h|θ) / t(θ)

[0258] In the formula, C(v,h) is the joint probability distribution value, I is the N-dimensional identity matrix, J represents the J samples, and e i v is the attenuation coefficient of the i-th source. i For i velocity sample points, f i Let h be the randomness value of i sample points. j Let ξ be the thickness of j sample points. ji Let p(v,h|θ) be the energy coefficient between the j-th sample point and the i-th sample, and let l be the calculated value of the energy function. -C(v,h|θ) Let t(θ) be the energy value calculated under the joint probability distribution, and t(θ) be the partition function.

[0259]

[0260] In the formula, v is the velocity and h is the thickness;

[0261] The likelihood function defined by RBM, p(v|θ) is the marginal distribution of the joint probability p(v,h|θ);

[0262] Since neurons within the RBM layer are unconnected, the operating state and activation of the hidden layers are relatively independent based on the visible layer state. The activation probability of the j-th hidden layer node is:

[0263]

[0264] In the formula, F'() is the current activation probability function, h jLet j be the thickness of the hidden layers, θ be the likelihood coefficient, σ() be the sigmoid function, and f j Let j be the randomness values ​​of the sample points;

[0265] Where σ(u)=1 / (1+l) -u Let be the sigmoid function; given the states of the hidden layer nodes, the activation probability of the i-th visible layer node is:

[0266]

[0267] In the formula, e i Let be the attenuation coefficient of the i-th earthquake source;

[0268] RBMs are trained and run iteratively, and the surface dilatation field during the run is studied and calculated using the parameter θ = (ξ). ij ,e i ,f j The value of ) is obtained from the given training data and training samples; the maximum log-likelihood function is calculated using the parameter θ, with the number of training samples set to T, and the expression is:

[0269]

[0270] In the formula, θ * To calculate the maximum log-likelihood function value using parameter θ, F'(v (T) |θ) is the activation probability function for the joint probability given the number of samples T in the training set;

[0271] (ii) Training process of DDN: The velocity change rate parameter passing through the plate reference frame between two adjacent continents is fused with the velocity field signal data feature quantity under the plate reference frame between two adjacent continents as input, and an RBM is trained from bottom to top each time; in each layer, the parameter space ωk is constructed by the numbers calculated in the (k-1)th layer, and the weights are updated according to the following formula:

[0272]

[0273] In the formula, ξ ij Here, θ represents the updated weight value, η() represents the update coefficient, and η() represents the characteristic function between two adjacent continents. <v i h j > data The current value of (v,h) <v i h j > model Let (v,h) be the standard value.

[0274] In step 5, the tortuosity variation optimization method is used to optimize the biases a and b of the visible and hidden layers of the DDN, as well as the parameters of the number of network layers and nodes, including:

[0275] (a) Initialization of the parameters of the fault change optimization method and the optimization problem. The optimization problem is described as follows:

[0276] Minimize f(u) subject to u j ∈E j

[0277] In the formula, Minimize f(u) is the initial value of the optimization surface expansion rate field, u j is the velocity of the decision vector for the j-th sample, f(u) is the function of the optimization surface expansion rate field, j is the signal length of the j-th sample, j = 1, 2, 3…N, and E j is the constraint interval of u after the decision vector j ;

[0278] (b) Set the initial plate and calculate the enthalpy value or entropy value. Use the unified population method to uniformly initialize the initial plate in the feasible search interval; Set two initial plates O0 = (w1, w2…w n ), O1 = (a1, a2…a n ), where n is the number of data points of the initial plate; Initialize the segmentation factor, let k = 1; Increment the segmentation factor by 1. If the segmentation factor k = 2, then generate two other different plates; If the value of the segmentation factor is k, then generate 2k - 2 different plates; If the initial plate category is R and the population size is F, then only need to satisfy R < F, and continuously add new plates to the set of the initial plates; Stop when R > F, thus initializing to obtain an initial population containing F plates;

[0279] (c) Simulate the process of the change of the surface expansion rate field, and encode the process of the change of the surface expansion rate field;

[0280] (d) Plate update: Update the plate by testing the change balance; The state when all surface expansions no longer change is the balanced state; In the balanced state, the plate with the increased entropy or decreased enthalpy of the surface expansion is used as the new plate, and the abnormal plates are excluded;

[0281] (e) Judge the iteration termination condition: If the condition is met, the algorithm terminates; otherwise, go to step (c); The termination condition of the fault change optimization method is to satisfy the maximum number of iterations or the minimum enthalpy value and the maximum entropy value.

[0282] Fault plane parameter inversion. Based on the precise source parameter information of 6,255 relocated seismic lines obtained by the double-difference seismic tomography method, and according to the principle that cluster earthquakes occur near faults, a combination of simulated annealing and Gauss-Newton algorithms was used to invert fault plane parameters of the seismic strip near Pinggang. The final results showed that the fault strike is approximately N78°E, dip angle is approximately 85°, dip direction is NW, fault length is approximately 16km, and fault depth is 4-13km.

[0283] (2.4) Study on the earthquake-generating mechanism, seismogenic structure and source rupture process in a certain city.

[0284] This invention established eight depth profiles along the NEE and NNW directions, and the distributions of their vertical velocity Vp and relative velocity perturbation ΔVp are as follows: Figure 9 As shown, the P-wave velocity in the shallow layer (above 5km) of a certain city area is generally high, reaching 5.5-6.0km / s. This may be because the Quaternary sedimentary layer in this area is relatively thin, with only 5-20m thick sedimentary layers along the river and bay, and Yanshanian granite, Indosinian granite and Cambrian metamorphic rocks are commonly exposed.

[0285] From the A-AA and B-BB cross sections, it can be seen that ( Figure 9 Figure a shows the velocity Vp and relative velocity disturbance ΔVp distribution on the NEE profile. Figure 9 (Figure b) A fault (F1) in a certain city is located at the boundary between high and low velocity anomalies. The western side has a low-velocity anomaly corresponding to the Yanshanian granite system, while the eastern side has a high-velocity anomaly corresponding to the Cambrian metamorphic rock system. Two earthquake swarms appear on either side of the fault zone. The earthquake distribution on the cross-section shows a near-vertical (slightly dipping SWW) downward extension. The earthquake swarm on the western side of the fault is a newly formed swarm after 2014, with a focal depth ranging from 2 to 10 km. The western swarm is located near a reservoir, with a focal depth ranging from 0 to 15 km. The reservoir area is 17 km². 2 The annual water level variation is only 2 meters, ruling out the reservoir's influence on the seismic activity of the earthquake swarm. Earthquakes near the reservoir are mainly located within a high-velocity body containing Cambrian metamorphic rocks, beneath which lies a low-velocity anomaly. This velocity structure is conducive to energy accumulation, but given the small amplitude of the high-velocity anomaly (approximately 2%), the ultimate stress that the rocks can withstand is not high, making a destructive earthquake unlikely. The largest earthquake to date was M on March 20, 2018. L A magnitude 4.2 earthquake.

[0286] From the C-CC and D-DD cross sections ( Figure 9 Figure C in the middle, Figure 9(See Figure d in the middle). The 6.4 magnitude earthquake in 1969 occurred in a certain city. The 6.4 magnitude earthquake was located within a high-velocity anomaly body sandwiched by a low-velocity anomaly in the Yanshanian granite system. The low-velocity anomaly is located below the high-velocity body. This structure has also appeared in the earthquakes in the city's jurisdiction in 1976, the first national earthquake in 1995, the second national earthquake in 2001, the urban area under the jurisdiction of a certain province in 2008, and the double earthquake in a certain district of a certain province in 2012. After the 6.4 magnitude earthquake, the aftershocks were mostly concentrated within the high-velocity body surrounded by the aforementioned low-velocity body. The low-velocity layer below the height anomaly body where the 6.4 magnitude earthquake occurred is a ductile shear layer. The melting of the mantle wedge and the underplating of basalt caused partial melting of the rocks in this layer, which showed the low-velocity anomaly. The hot springs exposed on the surface also provide evidence for this. Based on focal mechanism analysis, the 6.4 magnitude earthquake in a certain city was a strike-slip earthquake, with two subsurfaces trending NNW and NEE respectively. These subsurfaces align with the strikes of the Yangbianhai Fault and the Pinggang Fault, respectively. Combining velocity structure, fault distribution, focal mechanism, velocity field characteristics, and post-earthquake field investigation data, it is inferred that the earthquake's gestation process was influenced by the combined effects of the Pinggang Fault and the Yangbianhai Fault. The velocity difference between the medium below the seismogenic layer and the seismogenic layer itself is mostly between 0.1 and 0.9 km / s; the smaller the difference, the higher the seismogenic capacity. Due to the presence of a low-velocity layer, the velocity difference between the medium below the seismogenic layer and the seismogenic layer in the 6.4 magnitude earthquake in this city was <0.0 km / s, indicating a strong seismogenic capacity. The intersection of the Pinggang Fault and the Yangbianhai Fault is conducive to earthquake occurrence, thus leading to the 6.4 magnitude earthquake. Furthermore, along the C-CC profile (… Figure 9 As shown in Figure c), there are significant differences in the distribution characteristics of earthquakes. With 40km as the boundary, the focal depth range near the Yangbianhai Fault on the west side is relatively large, extending from the surface to about 25km underground, while the shallow earthquake distribution along the Pinggang Fault on the east side is less, with a depth range of about 4-13km.

[0287] From the E-EE and F-FF cross sections ( Figure 10 Figures a and b show the velocity Vp and relative velocity disturbance ΔVp distribution on the NNW profile, combined with the B-BB profile (Figure a). Figure 9 In the earthquake projection of Figure b), the small earthquake swarms on the west and east sides of the Yangchun-Zhizhen (F1) fault both show high-angle seismic bands dipping towards the SE. The earthquake swarm on the west side is located at the junction of high and low velocity anomalies, and the high-velocity material in its SE disk is uplifted, indicating that it has reverse fault characteristics dipping towards the SE, which corresponds to the Indosinian reverse movement of the Yangchun-Zhizhen fault.

[0288] From the G-GG and H-HH cross sections ( Figure 10Figures c and d show the velocity Vp and relative velocity disturbance ΔVp distribution maps on the NNW-trending profile. Earthquakes near Pinggang are mainly concentrated below 3 km, forming a seismic band dipping NW. This band is located at the boundary between high and low velocity anomalies, and is biased towards the high-velocity anomaly on the northwest side. Based on the fault plane parameter inversion results, the band's strike is approximately N78°E, dip angle is approximately 85°, dips NW, fault length is approximately 16 km, and fault depth is 4-13 km. Combined with the field investigation results, no fault exposure was found on the surface. It is speculated that the fault containing the Pinggang earthquake band may be a concealed reverse fault, belonging to the southwest concealed segment of the Pinggang fault. This concealed segment dips NW, with the northwest block uplifted and the southeast block subsided. This is consistent with the tectonic characteristics of the Pinggang fault (F4) since the Miocene, where the northwest block has uplifted and the southeast block has subsided (Zhong Yijun and Ren Zhenhuan, 2003). Three earthquakes of magnitude 5 or greater have occurred at the intersection of the concealed segment of the Pinggang Fault and the Yangbianhai Fault: a magnitude 6.4 earthquake in 1969, a magnitude 5.0 earthquake in 1986, and a magnitude 5.0 earthquake in 1987. It is inferred that all three earthquakes were the result of the combined action of the Yangbianhai Fault and the Pinggang Fault. In contrast, the magnitude 4.9 earthquake in a certain city in 2004 originated from the southwest thrust concealed segment of the Pinggang Fault, and was less influenced by the Yangbianhai Fault. This could explain the difference in focal mechanism between the three earthquakes of magnitude 5.0 or greater and the magnitude 4.9 earthquake in the same city in 2004.

[0289] Furthermore, taking Pinggang as the boundary, the northeastern segment of the Pinggang fault dips SE at an angle of 60-90°, exhibiting normal fault characteristics. Near Pinggang, the dip becomes steeper, even vertical, extending southwestward beyond Pinggang but not exposed above the surface. The southwestern concealed segment dips NW at approximately 85°, exhibiting reverse fault characteristics. The normal fault characteristics of the northeastern segment of the Pinggang fault are inconsistent with the NW-SE trending compressional stress in the region. It is speculated that the fault planes of the northeastern and southwestern concealed segments of the Pinggang fault are twisted into a "braided" shape near Pinggang. The northeastern segment of the Pinggang fault was influenced by NW-SE trending extensional forces in the early neotectonic movement, resulting in a normal fault and the formation of a NE-trending trough. However, in the later stages of tectonic movement, the fault slickensides on the Pinggang fault surface show left-lateral displacement, indicating that the northeastern segment was later influenced by NW-SE trending compression.

[0290] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions, and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention, and within the spirit and principles of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for constructing a three-dimensional model of a seismogenic fault based on refined aftershock location data, characterized by, The method comprises the following steps: S1, collecting seismic and geological data information of the study area; S2, establishing a high-resolution three-dimensional crustal velocity model of the study area to determine accurate source parameters; S3, obtaining the characteristics of the tectonic stress field through source mechanism solutions and GNSS observation data; S4, predicting the seismogenic mechanism, seismogenic structure and source rupture process of the study area, combining the velocity structure, fault distribution, source mechanism, velocity field characteristics and post-earthquake field investigation data, predicting the influence of the seismogenic process of the study area under the joint action of faults, and thus obtaining potential strong earthquake hazard information; In step S3, the characteristics of the tectonic stress field are obtained through source mechanism solutions and GNSS observation data, including: S301, obtaining source mechanism solutions of multiple earthquakes; S302, studying the tectonic stress field of the study area according to the source mechanism solutions; S303, obtaining the velocity field and surface dilatancy rate field through GNSS data; S304, fault plane parameter inversion, obtaining accurate source parameter information after relocation according to the double-difference seismic tomography method, according to the principle that clustered earthquakes occur near faults, combining simulated annealing algorithm and Gauss-Newton algorithm, performing fault plane parameter inversion on the earthquake strip near a certain study area, and obtaining fault strike, dip angle, tendency NW, fault length and fault depth. 2.The method according to claim 1, wherein, In step S1, collecting seismic and geological data information of the study area, including: collecting seismic phase reports, seismic waveform data and GNSS observation data, establishing a basic database, collecting geological structure data of earthquake-prone areas in the study area through field investigation and previous data, and collecting geological structure data of earthquake-prone areas in the study area; In step S2, a high-resolution three-dimensional crustal velocity model of the study area is established to determine accurate source parameters, including: S201, using a grid interval of 0.2°x0.2° to perform double-difference seismic tomography tomoDD inversion to obtain the velocity structure MODEL_0.2 of the study area and the surrounding area; S202, performing linear interpolation on the velocity structure MODEL_0.2 of the study area and the surrounding area to obtain a velocity model MODEL_0.1 with a grid interval of 0.1°x0.1°; S203, using the velocity model MODEL_0.1 with a grid interval of 0.1°x0.1° as an initial model, performing inversion calculation using the double-difference seismic tomography method tomoDD to obtain the final three-dimensional crustal P-wave velocity structure MODEL_final of the study area and the positioning results of the earthquake. 3.The method according to claim 1, wherein, In step S302, the tectonic stress field of the study area is studied according to the source mechanism solutions, including: According to the magnitude, the source mechanism solutions are matched with weights, the source mechanism solutions of large earthquakes and small earthquakes are combined, and the given weight is obtained by the following formula: ； wherein is a given weight, is a given coefficient for attenuation of the source, is a coefficient for magnitude attenuation, is a maximum magnitude weight, is a current magnitude source, is a minimum magnitude source. 4.The method according to claim 1, wherein, In step S303, the velocity field and surface dilatancy rate field are obtained through GNSS data, including: Step 1, calculating the mean, randomness and stability of the velocity field data set under the plate reference frame between the adjacent two continents; according to the Euler parameters and single-day relaxation solution parameters, calculating the threshold values of the three characteristic quantities; Step 2, judging whether it is a surface dilatancy data according to the threshold value, and removing the background field data to retain the surface dilatancy data; Step 3, using the global linear embedding velocity rate of change between the two adjacent continents under the plate reference frame, the feature matrix of three-dimensional space is mapped to the low-dimensional surface inflation rate field embedding space RD, and the three feature quantities of each sample are data dimension reduction and parameter fusion; Step 4, the fused feature data is input into the deep belief network DDN, and the input parameters of DDN are optimized using the fault change optimization method; Step 5, in the training process of DDN, the bias a, b of the visible layer and the hidden layer, the number of network layers and the number of nodes of DDN are optimized using the fault change optimization method; after data training, three mathematical models corresponding to the surface inflation are obtained, and the construction of the velocity field data self-checking system under the plate reference frame between the two adjacent continents is completed; Step 6, the signal of the velocity field data acquisition system is input into the velocity field data self-checking system, the data is processed by the self-checking system, the processing result is input into the mathematical model of the velocity field data self-checking system, and the surface inflation type of the velocity field signal acquisition system under the plate reference frame between the two adjacent continents is determined.

5. The method according to claim 4, wherein, In step 1, the mean, randomness and stability of the velocity field data set under the plate reference frame between the two adjacent continents are calculated, including: (1) Calculate the mean using the following formula : ； In the formula, is the total number of sample points, is sample points; (2) Calculation of randomness: the collected signals are arranged in descending order of numerical value, all data of the same size are combined, and the number of remaining data is the data randomness index ; (3) the stability is calculated according to the following formula: ； In the formula, is the current normal signal is the stable value of the sample point of the continent, is the total normal signal is the stable value of the sample point of the continent, is the first continent, is the difference value of the sample point between the adjacent two continents, is the previous continent, is a certain adjacent two continents, is the positive feature quantity of the square root of the stable value of the normal signal, is the data stability index; Computing the mean threshold with the expression ； Computing a randomness threshold , the expression is: ； Computing stability threshold , the expression is: ； In the formula, is a speed sample point, is the signal collected in numerical order from large to small, and the data retained after merging all data of the same size, is the stability characteristic quantity of the normal signal. 6.The method according to claim 4, wherein, In step 3, the velocity change rate between the adjacent two continents under the plate reference frame is embedded using a global linear embedding, including: mapping the three-dimensional spatial feature matrix to a low-dimensional surface inflation rate field embedding space In the foregoing, the three feature quantities of each sample are subjected to data dimension reduction and parameter fusion; the specific steps are as follows: (1) select a local neighborhood, for a given data set, then: ； wherein is a data set for a local neighborhood, is the data for the local neighborhood, is the data set for the local neighborhood; Find each sample point of the neighborhood nearest neighbors using the Euclidean distance formula: ； In the formula, For the first The sample point and the first Euclidean distance of each sample point For the first Each sample point in the neighborhood speed value, For the first Each sample point in the neighborhood speed value, For the first Neighborhood, for Each neighborhood, This is the magnitude attenuation coefficient; (2) calculate the reconstruction weight of the sample point neighborhood, construct the local reconstruction weight matrix, and define the error function as: ； wherein is the weight error, is sample points, is sample points, is the second order absolute value, is and the weight between is the nearest neighbor points, ; To limit the condition, the error function expression is transformed to: ； In the formula, For the first Local reconstruction weights for each sample point , For the first Local reconstruction weights of each neighborhood, For iteration values; For the first The local covariance matrix is ​​obtained from each sample point. ; Introduce the Lagrange multiplier to solve the constraint problem, the expression is: ； wherein is a constraint weight, is a Lagrange multiplier, is the i-th sample point, is the j-th sample point, and is a constraint weight. Let , readjust the weight, let the sum of the weight be 1, finally get ; (3) The weight matrix is obtained Finding a low-dimensional manifold embedding Z of the sample set that minimizes the reconstruction error and a function, expressed as: ； wherein is a reconstruction error and function minimization value, is is a local covariance matrix of the sample points; Restrict Z, the expression is: ； In the formula, is identity matrix; The optimization problem is transformed into a constrained optimization problem, the expression is: ； wherein is the is the is the is the local reconstruction weight of the is the constraint optimization value of the Gth iteration, is the data set after the Gth iteration, is the data set; Data set is equivalent to finding the eigenvectors of the symmetric, positive semi-definite, sparse matrix ​ ； In the formula, is a solution of the constraint weight matrix; Using the Lagrange multiplier method, we get: ； wherein is a constraint optimization value obtained by Lagrange multiplier method, is a constraint weight; To partial derivative minimization one obtains: ； where, , find i.e. find the eigenvectors of M, so that The obtained embedding coordinates are the eigenvectors of M; the eigenvector corresponding to the smallest nonzero eigenvalue is taken as the value of M, and the low-dimensional face inflation rate field coordinates The eigenvector corresponding to the eigenvalue is the output result; The whole process of the velocity rate of change between the two adjacent continents under the plate reference frame is shown in the following formula: ； wherein, is the velocity set under the plate reference frame between the two adjacent continents, is the local reconstruction velocity weight set of the th sample point, is the velocity change data set after embedding the low-dimensional surface inflation rate field, is the rate of change.

7. The method according to claim 4, wherein, In step 4, the fused feature data is input into the deep belief network DDN, and the input parameters of DDN are optimized using the fault change optimization method, including: (i) construct a deep belief neural network DDN; Based on the energy function formula set up by calculation, the joint probability distribution of is calculated as: ； In the formula, is a joint probability distribution value, is is a unit matrix, is is a sample, is the first is a source attenuation coefficient, is is a velocity sample point, is is a sample point randomness value, is is a sample point thickness, is the energy coefficient between the first is a sample point and the first is a sample, is an energy function calculation value, is an energy value calculated under a joint probability distribution, is a partition function; ； wherein is the velocity, is the thickness; the likelihood function defined by the RBM, the joint probability the marginal distribution; Since there is no connection between neurons in the RBM layer, the activation of the hidden layer is relatively independent of the activation of the visible layer, and the activation probability of the first hidden layer node is: ​ ； wherein, is the current activation probability function, is a hidden layer thickness, is a likelihood coefficient, is a sigmoid function, is a sample point randomness value; wherein, sigmoid function; Given the state of a hidden layer node, the activation probability of a visible layer node is: ​ ； In the formula, is the first source attenuation coefficient; The RBM is trained to run in an iterative manner, and the surface expansion rate field is calculated by studying the values of the calculation parameters through the given training data and training samples; the maximum log-likelihood function is calculated through the parameters , and the number of training set samples is set as , and the expression is: ； In the formula, is the parameter The maximum log-likelihood function value is calculated, is the number of training set samples The activation probability function of the lower joint probability; (ii) Training process of DDN: the rate of change of velocity parameters between adjacent two continents under the plate reference frame, and the signal data characteristics of the velocity field between adjacent two continents under the plate reference frame are fused as inputs, and an RBM is trained from bottom to top each time; in each layer, the parameter space is constructed by the number calculated by the first layer, and the weight is updated according to the following formula: ​​ ； wherein is a weight update value, is an update coefficient, is a characteristic function between two adjacent continents, is a current value, is a table reference value. 8.The method according to claim 4, wherein, In step 5, the bias a, b of the visible layer and the hidden layer, the number of network layers and the number of nodes of DDN are optimized using the fault change optimization method, including: (a) initialization of fault change optimization method parameters and optimization problem, the optimization problem is described as: ； wherein is the optimized face expansion rate field initialization value, is is the decision vector velocity of the sample, is the optimized face expansion rate field function, is is the signal length of the sample, , is the decision vector post constraint interval; (b) Set initial patches and calculate enthalpy or entropy, use uniform method to initialize initial patches in feasible search interval; set two initial patches , for the number of data points of initial patches; initialize the split factor, let ; let the split factor plus 1, if the split factor , then generate two different patches; if the value of the split factor is , then generate different patches; if the initial patch category is and the population size is , then only need to satisfy , continuously increase new patches to the initial patch set; when , stop, thereby the initial population containing patches is initialized; (c) simulate the surface inflation rate field change process, encode the surface inflation rate field change process; (d) plate update: update the plate by testing the change balance; the state when all surface inflations no longer change is the equilibrium state; in the equilibrium state, the surface inflation with increased entropy or decreased enthalpy is taken as the new plate, and the abnormal plate is excluded; (e) judge the iteration termination condition: if the condition is met, the algorithm terminates, otherwise go to step (c); the termination condition of the fault change optimization method is to meet the maximum number of iterations or the minimum enthalpy and the maximum entropy.

9. A system for constructing a three-dimensional model of a causative fault based on aftershock fine positioning data, characterized by, The system implements the method for constructing a three-dimensional model of a seismogenic fault based on precise positioning data of aftershocks according to any one of claims 1-8, and the system comprises: An information acquisition module for acquiring seismic and geological data information of a study area; The crust three-dimensional velocity model establishing module is configured to establish a high-resolution crust three-dimensional velocity model of a research area and determine accurate source parameters. The tectonic stress field characteristic acquisition module is configured to acquire tectonic stress field characteristics through source mechanism solutions and GNSS observation data. The tectonic stress field characteristics are acquired through source mechanism solutions and GNSS observation data, including: acquiring source mechanism solutions of multiple earthquakes; researching tectonic stress fields of the research area according to the source mechanism solutions; acquiring velocity fields and surface expansion rate fields through GNSS data; fault surface parameter inversion, acquiring accurate source parameter information after repositioning according to a double-difference seismic tomography method, according to the principle that clustered earthquakes occur near faults, combining a simulated annealing algorithm and a Gauss-Newton algorithm, performing fault surface parameter inversion on an earthquake strip near a research area, and obtaining a fault strike, a dip angle, a tendency NW, a fault length, and a fault depth; The research area seismogenic mechanism, seismogenic structure, and source rupture process prediction module is configured to combine velocity structures, fault distributions, source mechanisms, velocity field characteristics, and post-earthquake field investigation data to predict that a seismogenic process of an earthquake is affected by the joint action of faults, thereby obtaining potential strong earthquake risk information.

Images