Earthquake-generating fault three-dimensional model construction method and system based on aftershock fine positioning data

Abstract

The invention belongs to the technical field of seismic information processing, and discloses a method and a system for constructing a three-dimensional model of a seismic fault based on aftershock fine positioning data. The method comprises the following steps: acquiring seismic geological data information of a research area; establishing a high-resolution earth crust three-dimensional speed model in a research area and determining accurate seismic source parameters; obtaining tectonic stress field characteristics through a seismic source mechanism solution and GNSS observation data; a regional pregnancy and earthquake mechanism is researched, and potential strong earthquake risk information is obtained. According to the method, the grids with the interval of 0.2 degrees * 0.2 degrees and the grids with the interval of 0.1 degrees * 0.1 degrees are combined, the speed structure of the interval of the grids with the interval of 0.2 degrees * 0.2 degrees in the non-seismic source area can be obtained, the speed structure of the interval of the grids with the interval of 0.1 degrees * 0.1 degrees in the seismic source area can also be obtained, and the method is suitable for research on the speed structure of the seismic source area with a relatively small number of stations.

Description

Technical Field

[0001] The present invention belongs to the technical field of earthquake information processing, and in particular relates to a method and system for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data. Background Art

[0002] A certain city in a certain province is located at the southwestern end of the Shaoguang depression series in the earth dome system and at the southern end of the western edge of the second giant uplift belt in the tectonic system. The earthquake geological structure is relatively complex and it is one of the most active areas in the coastal seismic belt. L Earthquake activities above magnitude 3.0 can serve as earthquake "windows" in coastal seismic zones. Therefore, studying the underground medium structure and seismic activity in a certain area is of great significance for the assessment of earthquake risks in economically developed coastal areas.

[0003] The above analysis reveals the following problems and drawbacks of the existing technology: The three-dimensional crustal P-wave velocity structure and earthquake location prediction results in the existing earthquake analysis technology are low, and cannot provide a theoretical basis for earthquake risk assessment. Summary of the Invention

[0004] To overcome the problems existing in the related art, the disclosed embodiments of the present invention provide a method and system for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data.

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

[0006] S1, collect earthquake and geological data information in the study area;

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

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

[0009] S4, predict the earthquake-pregnancy mechanism, seismogenic structure and source rupture process in the study area, combine the velocity structure, fault distribution, focal mechanism, velocity field characteristics and post-earthquake field investigation data to predict the impact of the combined action of faults on the earthquake-pregnancy process, and thus obtain information on potential strong earthquake hazards.

[0010] In step S1, seismic and geological data information of the study area is collected, including: collecting earthquake phase reports, seismic waveform data, GNSS observation data to establish a basic database, and conducting field surveys and reviewing previous data to collect geological and structural 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 the precise earthquake source parameters, including:

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

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

[0014] S203: Using the velocity model MODEL_0.1 with a grid interval of 0.1°×0.1° as the initial model, the double-difference seismic tomography method tomoDD 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 tectonic stress field characteristics are obtained through the focal mechanism solution and GNSS observation data, including:

[0016] S301, obtaining focal mechanism solutions of multiple earthquakes;

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

[0018] S303, obtaining the velocity field and the surface expansion rate field through GNSS data;

[0019] S304, fault plane parameter inversion, based on the precise source parameter information of 6255 relocated earthquakes obtained by double-difference seismic tomography, and based on the principle that clustered earthquakes occur near faults, combined with the simulated annealing algorithm and the Gauss-Newton algorithm, fault plane parameter inversion is performed on the seismic belt near a certain study area to obtain the fault strike, dip, NW dip, fracture 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 the magnitude, the focal mechanism solutions are weighted and the focal mechanisms of large earthquakes and small earthquakes are combined. The given weights are obtained by the following formula:

[0022] ω=e r / D 2

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

[0024] The magnitude attenuation coefficient D is calculated based on the maximum and minimum regional magnitudes, and the expression is:

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

[0026] In the formula, ω is the given weight, e r is the source attenuation coefficient, D is the magnitude attenuation coefficient, ω max is the maximum magnitude weight, M is the current magnitude source, and M min is the minimum 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 intervals of the three rotation angles of the stress field rotation axis are all 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 on 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 plane 1 and the strike, dip, and slip angle of nodal plane 2 under the optimal stress state.

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

[0030] Data from GNSS stations in the study area were collected and processed with high precision using GAMIT / GLOBK software to obtain a single-day relaxation solution. The single-day relaxation solution files published by the IGS were then merged for overall adjustment to obtain the velocity field in the ITRF framework. The velocity field in the ITRF framework was converted into a velocity field relative to the plates between adjacent continents. The reference frame between adjacent continents was defined by differentially comparing the velocity in the ITRF with the rotation of the Eurasian Plate. When performing the reference frame conversion, the Euler parameters of the plates between adjacent continents obtained from the global geological model NUVELL-1A were used to perform the Euler transformation, remove the background field of the Eurasian Plate motion, and obtain the velocity field in the plate reference frame relative to the adjacent continents.

[0031] The surface expansion rate field of the study area is obtained from the velocity field calculation, including:

[0032] Step 1: Calculate the mean, randomness, and stability of the velocity field dataset between two adjacent continents in the plate reference frame; calculate the thresholds of the three characteristics based on the Euler parameters and the single-day relaxation solution parameters;

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

[0034] Step 3: Use global linear embedding to calculate the velocity change rate in the plate reference frame between two adjacent continents, map the feature matrix of the three-dimensional space to the low-dimensional surface expansion rate field embedding space RD, and perform data dimensionality reduction and parameter fusion on the three feature quantities of each sample;

[0035] Step 4: Input the fused feature data into the deep belief network (DDN) and use the fault variation optimization method to optimize the input parameters of the DDN.

[0036] Step 5: During the DDN training process, the fault variation optimization method is used to optimize the parameters of the DDN's visible and hidden layer biases a and b, as well as the number of network layers and nodes. After data training, mathematical models corresponding to the three surface expansions are obtained, completing the construction of a self-verification system for velocity field data in the plate reference frame between two adjacent continents.

[0037] In step 6, the signals of the velocity field data acquisition system in the plate reference frame between the two adjacent continents to be diagnosed are input into the velocity field data self-verification system in the plate reference frame between the two adjacent continents. After the system processes the data, the processing results are input into the mathematical model of the velocity field data self-verification system in the plate reference frame between the two adjacent continents, thereby completing the surface expansion type discrimination of the velocity field signal acquisition system in the plate reference frame between the two adjacent continents.

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

[0039] (1) Use the following formula to calculate the mean:

[0040]

[0041] Where F1 is the mean of the velocity field data set, S is the total number of sample points, and u i is i sample points;

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

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

[0044]

[0045] Where, Let \(S_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 \(F_3\) be the data stability index;

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

[0047]

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

[0049]

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

[0051]

[0052] In the formula, \(v\) i is the \(i\)th speed sample point, \(f\) s is the data obtained by arranging the collected signals from large to small and retaining the data after merging all data with the same size, is the stability characteristic quantity of the normal signal.

[0053] In step 3, use the rate of change of velocity in the plate reference frame between two adjacent continents with global linear embedding, including: mapping the three-dimensional space feature matrix to the low-dimensional surface expansion 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: ​​​​​​​​​​​​​​​​​​​​​​​​​​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 in the neighborhood k, u jk is the velocity value of the jth sample point in the kth neighborhood, k is the kth neighborhood, N is the Nth neighborhood, and D is the magnitude attenuation coefficient;

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

[0061]

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

[0063] To constrain the conditions, the error function expression is transformed into:

[0064]

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

[0066] The Lagrange multiplier is introduced to solve the constraint problem, and the expression is:

[0067]

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

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

[0070] (3) The low-dimensional surface expansion rate field of the sample set is found through the obtained weight matrix R to embed Z, and the reconstruction error and function are minimized. The expression is:

[0071]

[0072] Where φ(Z) is the reconstruction error and function minimization value, z j is the local covariance matrix of j sample points;

[0073] With restrictions on Z, the expression is:

[0074]

[0075] Where I is the N-dimensional identity matrix;

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

[0077]

[0078] Where, I i is the i-th unit matrix, R i is the local reconstruction weight of the i-th sample point, ZMZ G is the constrained optimization value of G iterations, Z G is the data set after G iterations, Z is the data set;

[0079] The dataset Z is equivalent to finding the eigenvectors of the symmetric, semi-positive definite, sparse matrix M, and we get:

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

[0081] Where W is the constraint weight matrix to be solved;

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

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

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

[0085] Taking the partial derivative of Z to minimize L(Z) yields:

[0086]

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

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

[0089]

[0090] Where, is the velocity set between two adjacent continents in the plate reference frame, is the local reconstruction speed weight set of the i-th sample point, is the velocity change data set after embedding the low-dimensional surface expansion rate field, → is 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 the fault variation optimization method, including:

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

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

[0094]

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

[0096] Where C(v,h) is the joint probability distribution value, I is the N-dimensional unit matrix, J is the J samples, e i is the attenuation coefficient of the i-th earthquake source, v is the i-th velocity sample point, f i is the randomness value of the i sample point, h j is the thickness of j sample points, ξ ji is the energy coefficient between the jth sample point and the ith sample, p(v,h|θ) is the calculated value of the energy function, l -C(v,h|θ) is the energy value calculated under the joint probability distribution, t(θ) is the partition function;

[0097]

[0098] Where v is velocity and h is thickness;

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

[0100] Since neurons in the RBM layer are not connected, the operating state and activation of the hidden layer are relatively independent according to the state of the visible layer. The activation probability of the jth hidden layer node is:

[0101]

[0102] Where F'() is the current activation probability function, h j is the thickness of j hidden layers, θ is the likelihood coefficient, σ() is the sigmoid function, f j is the randomness value of j sample points;

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

[0104]

[0105] Where, e i is the attenuation coefficient of the ith earthquake source;

[0106] RBM is trained and operated in an iterative manner. The field of the operating surface expansion rate is to study the calculation parameter θ=(ξ ij ,e i ,f j ), given training data and training samples; the maximum log-likelihood function is calculated by parameter θ, the number of training set samples is set to T, and the expression is:

[0107]

[0108] Where θ * To calculate the maximum log-likelihood function value through the parameter θ, F'(v (T) |θ) is the activation probability function of the joint probability under the number of training set samples T;

[0109] (ii) DDN training process: The velocity change rate parameter in the plate reference frame between two adjacent continents is integrated with the velocity field signal data feature in the plate reference frame between two adjacent continents as input, and one RBM is trained at a time from the bottom up. At each layer, the parameter space ωk is constructed by the values ​​calculated in the k-1th 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 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 in 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, and 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 the decision vector after u j ;

[0116] (b) Set the initial plate and calculate the enthalpy or entropy value, and use the unified population method to uniformly initialize the initial plate in the feasible search interval; set two initial plates O0 = (w1, w_{2}... w n ), O1 = (a1, a_{2}... a n ), n is the number of data points of the initial plate; initialize the segmentation factor, let k = 1; let the segmentation factor be incremented 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 2^{k - 2} different plates; if the initial plate category is R and the population size is F, then only need to satisfy R < F, continuously add new plates to the set of the initial plates; stop when R > F, so as to initialize and obtain an initial group containing F plates;

[0117] (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;

[0118] (d) Plate update: The plate is updated by testing the equilibrium of changes; the state when all surface expansions no longer change is the equilibrium state; in the equilibrium state, the plate that increases the entropy or decreases the enthalpy of surface expansion is considered a new plate, and abnormal plates are excluded;

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

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

[0121] Another object of the present invention is to provide a system for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data, the system implementing the method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data, the system comprising:

[0122] Information collection module, used to collect earthquake and geological data information in the study area;

[0123] The module for establishing a 3D crustal velocity model is used to establish a high-resolution 3D crustal velocity model of the study area and determine the precise earthquake source parameters.

[0124] Tectonic stress field feature acquisition module, used to obtain tectonic stress field features through focal mechanism solutions and GNSS observation data;

[0125] The module for predicting earthquake-pregnancy mechanism, earthquake-generating structure and source rupture process in the research area is used to combine velocity structure, fault distribution, source mechanism, velocity field characteristics and post-earthquake field investigation data to predict the impact of the combined action of faults on the earthquake-pregnancy process, thereby obtaining information on potential strong earthquake hazards.

[0126] In combination with all the above technical solutions, the beneficial effects of the present invention are as follows:

[0127] The present invention first uses double-difference seismic tomography (tomoDD) inversion with a grid interval of 0.2°×0.2° to obtain the velocity structure of the study area and surrounding areas (MODEL_0.2). MODEL_0.2 is linearly interpolated to obtain a velocity model with a grid interval of 0.1°×0.1° (MODEL_0.1). MODEL_0.1 is used as the initial model and the double-difference seismic tomography method (tomoDD) is used again for inversion calculation to obtain the final velocity structure of the study area (MODEL_final).

[0128] The present invention combines two grid spacings, 0.2°×0.2° and 0.1°×0.1°, to obtain both the velocity structure of the 0.2°×0.2° grid spacing in the non-seismic area and the velocity structure of the 0.1°×0.1° grid spacing in the seismic area. This method is suitable for studying the velocity structure in seismic area with a relatively small number of stations.

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

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

[0131] Figure 1 This is a flow chart of a method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data provided by an embodiment of the present invention;

[0132] Figure 2 0.2° grid interval detection board inspection result diagram provided by the prior art; wherein, (a) is the detection board inspection result diagram of Z = 0km, (b) is the detection board inspection result diagram of Z = 5km, (c) is the detection board inspection result diagram of Z = 9km, (d) is the detection board inspection result diagram of Z = 14km, (e) is the detection board inspection result diagram of Z = 20km, and (f) is the detection board inspection result diagram of Z = 25km;

[0133] Figure 3 0.1° grid interval detection plate test result diagram provided by the prior art; wherein (a) is the detection plate test result diagram for Z = 5km, (b) is the detection plate test result diagram for Z = 9km, (c) is the detection plate test result diagram for Z = 14km, and (d) is the detection plate test result diagram for Z = 20km;

[0134] Figure 4 Velocity structure inversion flow chart provided by the present invention;

[0135] Figure 5 The travel time residual distribution histogram before (hollow bar) and after (solid bar) positioning provided by the present invention;

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

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

[0138] Figure 8 This is a diagram showing the inversion results of the stress field in a certain city area of ​​the present invention;

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

[0140] Figure 10 is the velocity V on the NNW section of the present invention p and relative velocity disturbance ΔV p Distribution diagram; among them, (a1) diagram is the velocity V on the NEE profile p E-EE cross-section diagram, (a2) is the relative velocity disturbance ΔV p E-EE cross-section diagram, (b1) shows the velocity V on the NEE cross-section p F-FF cross-section diagram, (a2) is the relative velocity disturbance ΔV p F-FF cross-section diagram, (c1) is the velocity V on the NEE cross-section p G-GG cross-section diagram, (c2) is the relative velocity perturbation ΔV pG-GG cross-section diagram, (d1) shows the velocity V on the NEE cross-section p H-HH profile, (d2) is the relative velocity perturbation ΔV p H-HH cross-section diagram. DETAILED DESCRIPTION

[0141] To make the above-mentioned objects, features, and advantages of the present invention more readily apparent, specific embodiments of the present invention are described in detail below with reference to the accompanying drawings. The following description sets forth numerous specific details to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art may make similar modifications without departing from the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0142] Example 1, as Figure 1 As shown in FIG, a method for constructing a three-dimensional model of a seismogenic fault based on aftershock precise positioning data includes:

[0143] S1, collect earthquake and geological data information in the study area;

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

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

[0146] S4, predict the earthquake-pregnancy mechanism, seismogenic structure and source rupture process in the study area, combine the velocity structure, fault distribution, focal mechanism, velocity field characteristics and post-earthquake field investigation data to predict the impact of the combined action of faults on the earthquake-pregnancy process, and thus obtain information on potential strong earthquake hazards.

[0147] Exemplarily, the present invention further provides a system for constructing a three-dimensional model of a seismogenic fault based on aftershock precise positioning data, comprising:

[0148] Information collection module, used to collect earthquake and geological data information in the study area;

[0149] The module for establishing a 3D crustal velocity model is used to establish a high-resolution 3D crustal velocity model of the study area and determine accurate earthquake source parameters.

[0150] Tectonic stress field feature acquisition module, used to obtain tectonic stress field features through focal mechanism solutions and GNSS observation data;

[0151] The module for predicting earthquake-pregnancy mechanism, earthquake-generating structure and source rupture process in the research area is used to combine velocity structure, fault distribution, source mechanism, velocity field characteristics and post-earthquake field investigation data to predict the impact of the combined action of faults on the earthquake-pregnancy process, thereby obtaining information on potential strong earthquake hazards.

[0152] The present invention establishes a high-resolution three-dimensional crustal velocity model and obtains accurate source parameter information; combines the focal mechanism solutions of small and medium earthquakes with GNSS observation data to determine the tectonic stress field in the study area; and based on the coupling pattern between the three-dimensional crustal velocity structure, earthquake distribution and tectonic stress field, combines field geological survey data to study the earthquake-pregnancy mechanism, seismogenic structure and source rupture process, and judge the potential risk of strong earthquakes.

[0153] The present invention obtains the three-dimensional velocity structure of the crust with a resolution of 0.1°×0.1° in the earthquake source area and 0.2°×0.2° in the non-seismic source area, as well as the source parameters of 6255 earthquakes in a certain city. Based on the focal mechanism solutions of 64 earthquakes in the certain city and the GNSS data of 7 stations in the certain city, it is obtained that the direction of the maximum principal compressive stress in the main tectonic stress field in the certain city is NW-SE. Based on the three-dimensional velocity structure of the crust, the source parameter information and the coupling mode between the tectonic stress field, combined with the field geological survey data, the relationship between the velocity structure characteristics of the crust in the certain city and earthquake geology and faults in the certain city was analyzed, and the four M-type earthquakes in the history of the certain city were determined. L The seismogenic structure and earthquake-pregnancy mechanism of earthquakes above magnitude 5.0 provide important basis for assessing the earthquake risk in the region.

[0154] Example 2, as another embodiment of the present invention, further describes the method of constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data of the present invention from a theoretical basis.

[0155] Obtaining tectonic stress field characteristics through focal mechanism solutions and GNSS observation data: The focal mechanism solution for small and medium-sized earthquakes in a certain city is obtained by combining the P-wave first motion with the P-wave and S-wave amplitude ratio. Based on continuous GNSS observation data, the regional motion field and surface expansion rate field are fitted and calculated. Combined with the focal mechanism solution, the characteristics of the tectonic stress field in the city are studied to understand the stress accumulation and crustal movement within the area. Taking into account the precise source parameter information after relocation and based on the principle of small earthquake clustering, mathematical methods are used to re-invert the parameters of the earthquake-causing fault plane in the city. On this basis, the fault plane slip angle is calculated by referring to the tectonic stress field.

[0156] Study on the earthquake-pregnancy mechanism, seismogenic structure, and focal rupture process in a certain city: Field geological surveys were conducted in the earthquake-prone areas of a certain city, and geological conditions such as the geomorphological characteristics, lithologic distribution, planar distribution characteristics of exposed strata, and the location and occurrence of surface fault zones were investigated and studied. Combining the three-dimensional crustal wave velocity structure, earthquake distribution, tectonic stress field, and field geological survey data in the city, the focal structure in the city was located, and the relationship between velocity anomaly areas, seismicity, and tectonic stress fields was analyzed. The earthquake-pregnancy mechanism, seismogenic structure, and focal rupture process were studied. Comparisons were made with representative earthquakes in the southeastern coastal seismic belt and the eastern edge of the Qinghai-Tibet Plateau, and comparative analyses were conducted from the perspectives of crustal velocity structure, seismicity, fault structure characteristics, and crustal movement. This improved understanding of the earthquake-pregnancy mechanism, seismogenic structure, and focal rupture process in the study area, and further analyzed the potential strong earthquake hazard in the area.

[0157] Example 3, as another embodiment of the present invention, the main research progress, important results, key data, etc. and their scientific significance or application prospects are as follows.

[0158] (2.1) Detailed seismic geological conditions in a certain city. A field geological survey was conducted in the city, and previous data on the seismic geology of the city were collected. Ultimately, relatively detailed seismic geological data for the city were obtained. The city is located at the southwestern end of the Shaoguang depression series within the Zhejiang-Guangdong dome system and at the southern end of the western margin of the second giant uplift belt of the New Cathaysian tectonic system. The seismic geological conditions are relatively complex, with four major faults: the Yangchun-Zhilang fault (F1), the eastern branch of the Wuchuan-Sihui fault; the southwestern section 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-Jiling Fault (F1) generally strikes northeast or north-northeast, originating near Xinxing in the north and extending southwestward through Yangchun, Yangxi, and Shapa to estuary. It is presumed to continue southwestward. Its onshore section is 150 km long and dips primarily north-west, with areas dipping southeast at 60-80°. This fault originated during the Indosinian period and was a thrust fault. It was reactivated during the Yanshanian movement and converted into a normal fault. The southwestern section of the Cangcheng-Hailing Fault (F2) generally strikes northeast, originating northeast from southwest Enping, passing through Dahuai and Nalong to southeast of a certain city, crossing Hailing Island and venturing into the South China Sea. It is approximately 40 km long and dips northwest at a dip exceeding 60°. The Cangcheng-Hailing Fault is primarily composed of high-angle reverse faults, dissecting the Lower Paleozoic, Upper Paleozoic, Yanshanian granites, Cretaceous, and Paleogene, controlling the sedimentary boundaries of the Cretaceous and Paleogene rift basins. The Yangbianhai Fault, also known as the Fengtouhe Fault (F3), is exposed on the northeastern shore of the Yangbianhai Sea. It strikes NW and dips NE or SW at a dip of 65-75°. It is approximately 25 km long and can be divided into two sections: the northern and southern sections. The northern section extends primarily along the Yangbianhai Sea, largely submerged underwater. It is a sinistral normal fault, emerging only above water on the right shore of the Yangbianhai Sea. The southern section, found in villages such as Shajiao, Beiting, and Guangling on Hailing Island, features a well-defined fault slip surface, with scratches indicating a sinistral strike-slip pattern. The Pinggang Fault (F4) strikes NE to NEE, originating in the northeast and extending southwest into the Yangbianhai Sea. It is approximately 25 km long and is a typical right-lateral normal fault of the South China Sea system. The northeastern section dips SE at a dip of 60-70°, while the southwestern section steepens and even becomes vertical.

[0159] (2.2) Establishment of a high-resolution three-dimensional crustal velocity model for a certain city area and determination of precise earthquake source parameters.

[0160] In order to improve the resolution of the three-dimensional velocity structure of a certain city area, the present invention uses as many reliable seismic data as possible to participate in the inversion, and collects the P-wave arrival time data recorded by a certain provincial seismic network from January 1990 to August 2019. Among them, the accuracy of simulated earthquake records before June 2007, especially before 2000, is low. Therefore, the double-difference earthquake positioning method is used to relocate the earthquakes before June 2007 to obtain the relocated phase report. The above seismic data are further screened to eliminate the phase data that deviates from the theoretical travel time curve by more than 2s, and to ensure that each earthquake is recorded by at least 3 stations. Finally, 43,225 P-wave absolute arrival time data and 422,956 relative arrival time data of 6,390 earthquakes (736 earthquakes before June 2007) recorded by 49 stations are obtained for inversion calculation. The focal depth range is 0-32km, and the magnitude range of earthquakes before and after June 2007 is M, M, and M, respectively. L 0.9-4.4 and M L 0.2-4.2.

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

[0162] Table 1 One-dimensional P-wave velocity structure model of a city and its adjacent 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 of the detection board show that using a grid interval of 0.2°×0.2°, the restoration effects of land and sea areas in a certain city are significantly different, and the restoration effect of land areas at a depth of 5-31 km is better ( Figure 2 The test results of the test board with a 0.2° grid interval are shown in the figure), and the sea area has a good recovery effect only in the area close to the land. Using a 0.1°×0.1° grid interval, the land area has a good recovery effect at depths of 9km and 14km ( Figure 3(The test results of the detection board with a 0.1° grid spacing are shown in the figure). At depths of 5km and 20km, good recovery effects can only be obtained near the epicenter of a certain city (black box). This shows that using a 0.2°×0.2° grid spacing can achieve good resolution at all depth levels in the study area. However, given that the epicenter of a certain city covers an area of ​​approximately 0.2°×0.2°, using only a 0.2°×0.2° grid spacing cannot accurately reflect the situation in the epicenter. Using a 0.1°×0.1° grid spacing, a good resolution effect can be achieved in the epicenter of a certain city (black box), but the resolution effect outside the epicenter is not ideal.

[0165] To solve the above problem, the present invention first uses a 0.2°×0.2° grid interval to perform double-difference seismic tomography (tomoDD) inversion to obtain the velocity structure (MODEL_0.2) of a certain city and its surrounding areas, performs linear interpolation on MODEL_0.2, obtains a velocity model (MODEL_0.1) with a 0.1°×0.1° grid interval, uses MODEL_0.1 as the initial model, and again uses the double-difference seismic tomography method (tomoDD) for inversion calculation to obtain the final velocity structure (MODEL_final) of a certain city area. The above solution ( Figure 4 The velocity structure inversion flow chart shown in the figure combines the two grid intervals of 0.2°×0.2° and 0.1°×0.1°, which can obtain the velocity structure of the 0.2°×0.2° grid interval in the non-seismic area and the velocity structure of the 0.1°×0.1° grid interval in the seismic area. It is suitable for the study of the velocity structure in the seismic source area with a relatively small number of stations.

[0166] Using velocity model MODEL_0.1 as the initial model, an imaging inversion calculation with a grid spacing of 0.1°×0.1° was performed. After 20 iterations, the three-dimensional crustal P-wave velocity structure and the precise location of 6,255 earthquakes in a certain city were finally obtained. The root mean square residual of the earthquake travel time decreased from 365ms to 0.1ms, and the travel time residuals before and after inversion satisfied the Gaussian distribution ( Figure 5 As shown in the distribution histogram of travel time residuals before (hollow bars) and after (solid bars) positioning, the travel time residuals before inversion are between -0.6s and 0.6s, while the residuals after inversion are significantly reduced, mainly between -0.1s and 0.1s, among which the amount of data with residuals within 0.05s reaches 93%.

[0167] Velocity distribution in horizontal planes: After double-difference seismic tomography inversion calculation, the present invention finally obtains a high-precision three-dimensional P-wave velocity model (MODEL_final) of a certain city area. Figure 6 The velocity distribution diagram of 6 depth levels from 5 to 31 km is shown. Figure 6 In the figure (a), the five-pointed star represents the 6.4 magnitude earthquake in a certain city in 1969; Figure 6 In the figure (b), the five-pointed stars indicate the two magnitude 5.0 earthquakes in 1986 and 1987; Figure 6 The five-pointed star in the figure (c) indicates the 2004 4.9 magnitude earthquake; F1: Yangchun-Zhiluo fault; F2, southwestern section of Cangcheng-Hailing fault; F3, Yangbianhai fault; F4, Pinggang fault;

[0168] The velocity structure in the city and its surrounding areas exhibits significant lateral heterogeneity. The velocity distribution at a depth of 5 km clearly correlates with the fault structure and geological structure. The Mashui-Hekou-Tangkou-Pupai region lies within a high-to-low velocity transition zone. This transition zone coincides with the strike of the Yangchun-Zhilu fault north of Tangkou, while the fault deviates toward the high-velocity zone south of Tangkou. To the west of the transition zone, the Sanjia-Xinwei-Jueshan region exhibits a distinct low-velocity anomaly, corresponding to the Yanshanian granitic series in the area. This extends to 20 km below ground, before gradually transitioning to a high-velocity anomaly at a depth of 25 km in the Xinwei-Jueshan region. To the east of the transition zone, the Yangchun-Yangxi region exhibits a northeast-trending high-velocity anomaly corresponding to Cambrian metamorphic rocks. The southern segment of the Yangchun-Zhilu fault lies within this high-velocity anomaly. Two earthquake swarms occurred on either side of the fault zone: the western swarm lies within the high-to-low velocity transition zone, while the eastern Xinhu Reservoir swarm lies within the high-velocity zone. To the east of the NE-trending high-velocity anomaly lies a smaller, low-velocity anomaly corresponding to Yanshanian and Caledonian granites. However, a small area of ​​Cambrian metamorphic rocks is present near Chengcun, so the low-velocity zone appears to be interrupted by the high-velocity anomaly located on the Yangbianhai Fault. The Pinggang Fault lies at the junction of the high- and low-velocity anomalies, tending toward the high-velocity anomaly. Earthquakes are concentrated at the intersection of the Yangbianhai and Pinggang Faults, including the 1969 magnitude 6.4 earthquake.

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

[0170] 14km ​​depth ( Figure 6(c) in the figure), the Yangbianhai fault is located at the junction of high and low velocity anomalies, and turns into a low velocity anomaly near Pinggang. With the increase of depth, the low velocity anomaly near Pinggang expands and connects with the low velocity anomaly in the Sanjia-Xinwei-Jueshan area at 20km. Figure 6 (d) in the figure), the overall velocity structure shows a low-velocity anomaly, the NEE earthquake belt phenomenon of Yangbianhai-Pinggang disappears, and the earthquake swarm activity on both sides of the Yangchun-Zhiluo fault disappears. 25km to 31km depth ( Figure 6 In Figures (e) and (f), the 6.5 km / s contour line gradually becomes parallel to the coastline.

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

[0172] ① Focal mechanism solution results. This paper obtains focal mechanism solutions for 64 earthquakes, including 9 earthquakes before 1999 (M L The results (upper hemisphere projection) given by the focal mechanism solution of the southeastern coastal earthquakes with a magnitude of ≥3.6 are obtained. The data used are simulated network data. The focal mechanism solutions of these earthquakes are recalculated and the corresponding lower hemisphere projection parameters are obtained. L ≥3.0) focal mechanism solution for a certain province, the results obtained by using the snoke method. From 2004 to 2006, a total of 6 focal mechanism solutions were obtained, which were obtained by performing stress field inversion on a certain province and its neighboring areas. For the 42 M since 2007, L The focal mechanism solution of earthquakes with a magnitude of 3.6 or above is solved by an interactive program (Focmec-Interface) based on the Focmec method for inverting the focal mechanism solution.

[0173] According to the distribution map of focal mechanism solutions in a certain city, it can be seen that ML3.0 earthquakes in the city are mainly distributed in two areas, namely the intersection of the Yangbianhai Fault and the Pinggang Fault (46 focal mechanisms in Area A) and near the Xinhu Reservoir on the east side of the Yangchun-Zhiluo Fault (18 focal mechanisms in Area B).

[0174] The focal mechanism solution of Area A is relatively complex. The main earthquake in this area is right-lateral slip, accounting for 46%. Among them, 10 earthquakes have positive slip characteristics and 2 earthquakes have reverse slip characteristics. In addition, 11 earthquakes have reverse slip characteristics, accounting for 26%. The principal compressive stress axis P axis and the principal stress axis T axis do not have obvious convergence properties ( Figure 7(Figure 2 shows the distribution of P-axis and T-axis azimuths in Areas A and B). The four earthquakes with a magnitude of 4.9 or greater were all located in Area A, but the nature of the fault dislocations differed significantly. The largest earthquake was a magnitude 6.4 (ML 6.6) earthquake on July 26, 1969. Its focal mechanism was strike-slip, with two nodal planes trending N74°E and N20°W, respectively, consistent with the Yangbianhai Fault and the Pinggang Fault. The focal mechanism solutions for the 5.0 magnitude earthquake on January 28, 1986, and the 5.0 magnitude earthquake on February 25, 1987, are similar, both exhibiting strike-slip characteristics. The 1987 earthquake also exhibited a small amount of normal faulting. The two nodal planes for these two earthquakes trend NE and NW, respectively. The post-earthquake isoseismal lines for both earthquakes exhibit elliptical shapes with a NE-trending major axis. This suggests that the seismogenic structure for these two strong aftershocks was still the NEE-trending Pinggang Fault. However, if the Pinggang Fault were the seismogenic structure for both earthquakes, their left-lateral strike-slip characteristics would contradict the right-lateral strike-slip normal faulting of the Pinggang Fault. Seventeen years later, on September 17, 2004, another 4.9 magnitude earthquake occurred near Pinggang, northeast of the three aforementioned earthquakes. This earthquake exhibited reverse faulting, and the focal mechanism solution indicates a NEE-trending nodal plane, consistent with the major axis of the isoseismals and the strike of the Pinggang Fault.

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

[0176] The 18 earthquakes in Area B mostly occurred on the west side of the Yangchun-Zhiluo fault, especially near the Xinhu Reservoir. The focal mechanism solutions in this area are quite consistent, with 94% of the earthquakes exhibiting right-lateral strike-slip or normal-dip slip properties. The T-axis azimuth and dip are relatively stable, with the T-axis azimuth mainly located between NE16-40° and the dip mainly between 0-30°. The P-axis azimuth is mainly between 140±50° and 208±20°, but the dip varies greatly ( 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 was used to solve the tectonic stress field in a certain city. Because the number and quality of earthquake receiving stations vary for different magnitudes, the accuracy of the focal mechanism solution varies. This paper assumes that the larger the earthquake magnitude, the more accurate the focal mechanism solution. Therefore, the focal mechanism solution is weighted based on magnitude, thereby rationally combining the focal mechanisms of large and small earthquakes. The weights given by this paper are given by the following formula:

[0179]

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

[0181] The magnitude attenuation coefficient D is calculated based on the maximum and minimum regional magnitudes, and the expression is:

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

[0183] In the formula, ω is the given weight, e r is the source attenuation coefficient, D is the magnitude attenuation coefficient, ω max is the maximum magnitude weight, M is the current magnitude source, and M min is the minimum 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 intervals for the three rotation angles of the stress field rotation axis are all set to 1.0°, the search interval for 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 (see Table 2).

[0185] Table 2 Inversion results of tectonic stress field in a certain city

[0186]

[0187] Among them, the numerical range in brackets in Table 2 is the uncertainty range of each parameter under the 90 confidence level F test, and the value above it is the optimal stress state parameter.

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

[0189] Analysis suggests that the city has experienced multiple tectonic movements, such as the Caledonian, Indosinian, and Yanshanian movements, which produced strong NW-SE compression in the region. These residual compressive stresses still have a certain impact on the region's current tectonic stress field. According to geological survey results, the Yangchun-Jiling Fault generally strikes NE or NNE, predominantly 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 of over 60°. It is also composed of many high-angle reverse faults, indicating that the region was subject to NW-SE compression, consistent with the direction of the maximum NW-SE principal compressive stress in the tectonic stress field.

[0190] The Yangbianhai Fault strikes NW and dips NE. It is divided into two sections, south and north. The northern section mainly extends along the Yangbianhai (Fengtou River) and is a left-lateral normal fault. The scratches on the fault sliding surface show the left-lateral strike-slip nature. It can be seen that the Yangbianhai Fault is affected by the NE-SW tensile force, which is consistent with the NE-SW maximum dominant stress direction of the tectonic stress field.

[0191] Based on the tectonic stress field, it can be inferred that the right-lateral strike-slip dislocation of the NE-trending Pinggang Fault and the sinistral strike-slip of the NW-trending Yangbianhai Fault are the combined effects of NW-SE compression and NE-SW extension of tectonic stress. However, the normal faulting properties of these faults are inconsistent with the NW-SE compression of tectonic stress in this area, requiring further analysis.

[0192] ③ GNSS data to obtain velocity field and surface expansion rate field. The present invention collected data from 7 stations in a certain city area from 2018 to 2019, used GAMIT / GLOBK software for overall high-precision data processing to obtain a single-day relaxation solution, merged the single-day relaxation solution file published by IGS for overall adjustment, and finally obtained the velocity field under the ITRF framework. The velocity field under the ITRF framework is converted into a velocity field relative to the Eurasian plate, and the Eurasian reference frame is defined by differentially defining the rate in ITRF and the rotation of the Eurasian plate. When performing the reference frame conversion, the present invention uses the Euler parameters of the Eurasian plate obtained by the global geological model NUVELL-1A to perform Euler transformation, remove the background field of the Eurasian plate movement, and obtain the velocity field under the reference frame relative to the Eurasian plate.

[0193] The velocity field of the city reveals recent crustal movement characteristics in the study area. The movement direction at all seven stations is east-southeast, with a relatively small southward component, consistent with the NW-SE direction of the maximum principal compressive stress derived 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 in South China. This is primarily due to the combined effects of the SE-directed lateral compressive component derived from the strong uplift of the plateau under the compression of the Indian plate and the NW-SE-directed compression generated by the subduction of the oceanic plate under the Eurasian plate. The overall movement velocity of the city ranges from 1.02 to 11.26 mm / yr. The YCSJ, YJJZ, and YDDG stations in the north exhibit relatively high relative movement rates, potentially indicating localized internal locking, which favors strain accumulation. The surface expansion rate field for the study area, calculated from the velocity field, shows an overall low value and indicates the presence of a nearly N-S-trending contraction-expansion transition zone in the city, favoring stress accumulation.

[0194] Exemplarily, the calculation of the surface expansion rate field of the study area from the velocity field includes:

[0195] Step 1: Calculate the mean, randomness, and stability of the velocity field dataset between two adjacent continents in the plate reference frame; calculate the thresholds of the three characteristics based on the Euler parameters and the single-day relaxation solution parameters;

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

[0197] Step 3: Use the global linear embedding of the velocity change rate in the plate reference frame between two adjacent continents to map the three-dimensional feature matrix into the low-dimensional surface expansion rate field embedding space R D Above, the three feature quantities of each sample are subjected to data dimension reduction and parameter fusion;

[0198] Step 4: Input the fused feature data into the deep belief network (DDN) and use the fault variation optimization method to optimize the input parameters of the DDN.

[0199] Step 5: During the DDN training process, the fault variation optimization method is used to optimize the parameters of the DDN's visible and hidden layer biases a and b, as well as the number of network layers and nodes. After data training, mathematical models corresponding to the three surface expansions are obtained, completing the construction of a self-verification system for velocity field data in the plate reference frame between two adjacent continents.

[0200] In step 6, the signals of the velocity field data acquisition system in the plate reference frame between the two adjacent continents to be diagnosed are input into the velocity field data self-verification system in the plate reference frame between the two adjacent continents. After the system processes the data, the processing results are input into the mathematical model of the velocity field data self-verification system in the plate reference frame between the two adjacent continents, thereby completing the surface expansion type discrimination of the velocity field signal acquisition system in the plate reference frame between the two adjacent continents.

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

[0202] (1) Use the following formula to calculate the mean:

[0203]

[0204] Where F1 is the mean of the velocity field data set, S is the total number of sample points, and u i is i sample points;

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

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

[0207]

[0208] Where, is the stability value of the current n continent sample points of the normal signal, u(n) is the stability value of the total n continent sample points of the normal signal, n is the nth continent, u(η) is the difference between the sample points of two adjacent continents, τ is the previous continent, η is the difference between two adjacent continents, Q(n) is the positive characteristic value of the square root of the stability value of the normal signal, and F3 is the data stability index;

[0209] Calculate the mean threshold G1, the expression is:

[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 velocity sample point, f s is the data obtained by arranging the collected signals from large to small in numerical value and retaining the data after merging all the data with 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 The

[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] Where ε(R) is the weight error, u i is i sample points, u j are j sample points, is the second-order absolute value, R ij for u i with u ij The weight between ij for u i k nearest neighbors of , j = 1, 2, ..., k;

[0226] To constrain the conditions, the error function expression is transformed into:

[0227]

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

[0229] The Lagrange multiplier is introduced to solve the constraint problem, and the expression is:

[0230]

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

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

[0233] (3) The low-dimensional surface expansion rate field of the sample set is found through the obtained weight matrix R to embed Z, and the reconstruction error and function are minimized. The expression is:

[0234]

[0235] Where φ(Z) is the reconstruction error and function minimization value, z jis the local covariance matrix of j sample points;

[0236] With restrictions on Z, the expression is:

[0237]

[0238] Where I is the N-dimensional unit matrix; the optimization problem is transformed into a constrained optimization problem, expressed as:

[0239]

[0240] Where, I i is the i-th unit matrix, R i is the local reconstruction weight of the i-th sample point, ZMZ G is the constrained optimization value of G iterations, Z G is the data set after G iterations, Z is the data set;

[0241] The dataset Z is equivalent to finding the eigenvectors of the symmetric, semi-positive definite, sparse matrix M, and we get:

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

[0243] Where W is the constraint weight matrix to be solved;

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

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

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

[0247] Taking the partial derivative of Z to minimize L(Z) yields:

[0248]

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

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

[0251]

[0252] Where, is the velocity set between two adjacent continents in the plate reference frame, is the local reconstruction speed weight set of the i-th sample point, is the velocity change data set after embedding the low-dimensional surface expansion rate field, → is 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 the fault variation optimization method, including:

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

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

[0256]

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

[0258] Where C(v,h) is the joint probability distribution value, I is the N-dimensional unit matrix, J is the J samples, e i is the attenuation coefficient of the ith earthquake source, v i is the i velocity sample point, f i is the randomness value of the i sample point, h j is the thickness of j sample points, ξ ji is the energy coefficient between the jth sample point and the ith sample, p(v,h|θ) is the calculated value of the energy function, l -C(v,h|θ) is the energy value calculated under the joint probability distribution, t(θ) is the partition function;

[0259]

[0260] Where v is velocity and h is thickness;

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

[0262] Since neurons in the RBM layer are not connected, the operating state and activation of the hidden layer are relatively independent according to the state of the visible layer. The activation probability of the jth hidden layer node is:

[0263]

[0264] Where F'() is the current activation probability function, h jis the thickness of j hidden layers, θ is the likelihood coefficient, σ() is the sigmoid function, f j is the randomness value of j sample points;

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

[0266]

[0267] Where, e i is the attenuation coefficient of the ith earthquake source;

[0268] RBM is trained and operated in an iterative manner. The field of the operating surface expansion rate is to study the calculation parameter θ=(ξ ij ,e i ,f j ), given training data and training samples; the maximum log-likelihood function is calculated by parameter θ, the number of training set samples is set to T, and the expression is:

[0269]

[0270] Where θ * To calculate the maximum log-likelihood function value through the parameter θ, F'(v (T) |θ) is the activation probability function of the joint probability under the number of training set samples T;

[0271] (ii) DDN training process: The velocity change rate parameter in the plate reference frame between two adjacent continents is integrated with the velocity field signal data feature in the plate reference frame between two adjacent continents as input, and one RBM is trained at a time from the bottom up. At each layer, the parameter space ωk is constructed by the values ​​calculated in the k-1th layer, and the weights are updated according to the following formula:

[0272]

[0273] 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 standard value of (v,h).

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

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

[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 within 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; [[ID=X]]

[0280] (d) Plate update: Update the plate by testing the change balance; The state when all surface expansions no longer change is the balance state; In the balance state, the one that makes the surface expansion entropy increase or the enthalpy value decrease is used as the new plate, and 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 meet 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 earthquakes obtained using double-difference seismic tomography, and based on the principle that clustered earthquakes occur near faults, a simulated annealing algorithm combined with the Gauss-Newton algorithm was used to invert the fault plane parameters of the seismic belt near Pinggang. The final result was a fault strike of approximately N78°E, a dip of approximately 85°, a NW dip, a length of approximately 16 km, and a depth of 4-13 km.

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

[0284] The present invention establishes 8 depth profiles along the NEE and NNW directions, and the vertical velocity Vp and relative velocity disturbance △Vp are distributed as follows: Figure 9 As shown, it can be seen that the P-wave velocity in the shallow layer (above 5 km) of a certain city area is generally high, reaching 5.5-6.0 km / 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 the bay, and Yanshanian granite, Indosinian granite and Cambrian metamorphic rocks are generally exposed.

[0285] From the A-AA and B-BB sections, we can see that ( Figure 9 The velocity Vp and relative velocity disturbance △Vp distribution diagram on the NEE section shown in Figure a, Figure 9 (Figure b) A fault (F1) in a certain city is located at the junction of high and low velocity anomalies. The low velocity anomaly on the west side corresponds to the Yanshanian granite system, and the high velocity anomaly on the east side corresponds to the Cambrian metamorphic rock system. Two clusters of earthquakes appeared on both sides of the fault zone. The distribution of earthquakes showed a nearly vertical (slightly inclined to SWW) downward extension on the profile. The earthquake cluster on the west side of the fault is a new earthquake cluster after 2014, with a focal depth range of 2-10km. It is located at the junction of high and low velocity anomalies. The low velocity anomaly on the west side corresponds to the Yanshanian granite system, and the high velocity anomaly on the east side corresponds to the Cambrian metamorphic rock system. The earthquake cluster on the east side of the fault is located near the reservoir, with a focal depth range of 0-15km. The reservoir area is 17km 2 The annual water level change is only 2 meters, which excludes the influence of the reservoir on the seismic activity of the earthquake swarm. The earthquakes near the reservoir are mainly located in the high-speed body where the Cambrian metamorphic rock system is located, and below it is a low-speed anomaly. This velocity structure is conducive to energy accumulation. However, given that the amplitude of the high-speed anomaly is not large (about 2%), the ultimate stress that the rock can withstand is not high, and the possibility of a destructive earthquake is small. The largest earthquake that has occurred so far is M on March 20, 2018. L 4.2 magnitude earthquake.

[0286] From the C-CC and D-DD sections ( Figure 9 Figure C, Figure 9(Figure d in the middle) The 1969 6.4 magnitude earthquake occurred in a certain city. The 6.4 magnitude earthquake was located in a high-velocity anomaly body sandwiched by low-velocity anomalies in the Yanshanian granite system, with a low-velocity anomaly below the high-velocity body. This structure also appeared in the 1976 earthquake in a certain city, the 1995 earthquake in a certain country, the 2001 earthquake in a certain country, the 2008 earthquake in a city under the jurisdiction of a certain province, and the 2012 double earthquake in a certain district of a certain province. After the 6.4 magnitude earthquake, aftershocks were also concentrated in the high-velocity body surrounded by the above-mentioned low-velocity body. The low-velocity layer below the highly 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 the partial melting of the rocks in this layer to show low-velocity anomalies, and the hot springs exposed on the surface also provide evidence for this. According to the focal mechanism analysis, the 6.4-magnitude earthquake in a certain city was a strike-slip type, with two nodal planes in the NNW and NEE directions, which were consistent with the Yangbianhai Fault and the Pinggang Fault respectively. Combining the velocity structure, fault distribution, focal mechanism, velocity field characteristics and post-earthquake field investigation data, it is speculated that the seismogenic process of the 6.4-magnitude earthquake was affected by the combined effects of the Pinggang Fault and the Yangbianhai Fault. The difference between the velocity of the medium below the seismogenic layer and the velocity of the seismogenic layer is mostly between 0.1-0.9km / s. The smaller the difference, the higher the seismogenic capacity. Due to the existence of the low-velocity layer, the difference between the velocity of the medium below the seismogenic layer and the velocity of the seismogenic layer in the 6.4-magnitude earthquake in a certain city was <0.0km / s, and the seismogenic capacity was strong. The intersection of the Pinggang Fault and the Yangbianhai Fault was conducive to the occurrence of earthquakes, so a strong earthquake of magnitude 6.4 occurred. In addition, along the C-CC profile ( Figure 9 (Figure c in the middle) shows obvious differences in the distribution characteristics of earthquakes. With a depth of 40 km as the boundary, the focal depth range near the western Yangbianhai fault is relatively large, extending from the surface to about 25 km underground, while the shallow earthquakes along the eastern Pinggang fault are relatively rare, with a depth range of about 4-13 km.

[0287] From the E-EE and F-FF sections ( Figure 10 The velocity Vp and relative velocity disturbance △Vp distribution diagrams on the NNW section shown in Figure a and Figure b) are combined with the B-BB section ( Figure 9 In the seismic projection of Figure b), the small earthquake swarms on the west and east sides of the Yangchun-Zhilu (F1) fault both appear as high-angle seismic stripes dipping toward the SE. The earthquake swarm on the west side is located at the junction of high and low velocity anomalies, and its SE disk has high-speed material uplift, indicating that this area has a reverse fault nature dipping toward the SE, corresponding to the Indosinian thrust movement of the Yangchun-Zhilu fault.

[0288] From the G-GG and H-HH sections ( Figure 10As shown in the distribution of velocity Vp and relative velocity perturbation ΔVp on the NNW-striking profile (Figures c and d), earthquakes near Pinggang are primarily concentrated below 3 km, forming a NW-dipping seismic stripe located at the boundary between high- and low-velocity anomalies and biased toward the high-velocity anomaly to the northwest. Inversion of fault plane parameters indicates a strike of approximately N78°E, a dip of approximately 85°, a NW-dipping fault length of approximately 16 km, and a depth of 4-13 km. Based on field investigations, no surface faults have been observed within the stripe. It is speculated that the fault underlying the Pinggang seismic stripe is likely a buried reverse fault, belonging to the southwest buried segment of the Pinggang Fault. This buried segment dips NW, exhibits uplift on the northwest side, and subsidence on the southeast side. This is consistent with the tectonic characteristics of the Pinggang Fault (F4) since the Miocene, characterized by uplift on the northwest side and subsidence on the southeast side (Zhong Yijun and Ren Zhenhuan, 2003). Three earthquakes of magnitude 5 or greater have occurred at the intersection of the hidden 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 speculated that all three earthquakes were the result of the combined effects of the Yangbianhai and Pinggang Faults. The 2004 magnitude 4.9 earthquake in a certain city, however, originated from the hidden southwest thrust segment of the Pinggang Fault, which was less affected by the Yangbianhai Fault. This may explain the difference in focal mechanism between the three earthquakes of magnitude 5 or greater and the 2004 magnitude 4.9 earthquake in a certain city.

[0289] Furthermore, the northeastern section of the Pinggang Fault, bounded by Pinggang, dips SE at an angle of 60-90°, exhibiting normal fault properties. Near Pinggang, it steepens, even becoming upright, extending southwestward through Pinggang but not exposed. The southwest hidden section dips NW and dips approximately 85°, exhibiting reverse fault properties. The normal fault properties of the northeastern section of the Pinggang Fault are inconsistent with the NW-SE compression of tectonic stress in the area. It is speculated that the fault planes of the northeastern and southwest hidden sections of the Pinggang Fault twisted into a "twist-shaped" shape near Pinggang. The northeastern section of the Pinggang Fault was affected by NW-SE extension in the early stages of neotectonic movement, resulting in a normal fault and the formation of a NE-trending trough. However, in the late stages of tectonic movement, the fault shows left-lateral slip motion, as seen from the scratch surface of the Pinggang Fault, indicating that the northeastern section was later affected by NW-SE compression.

[0290] The above description is only a preferred specific implementation method 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 any technician familiar with this technical field within the technical scope disclosed by the present invention and within the spirit and principles of the present invention should be covered by the scope of protection of the present invention.

Claims

1. A method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data, characterized in that: The method comprises the following steps: S1, collect earthquake and geological data information in the study area; S2, establish a high-resolution three-dimensional crustal velocity model of the study area and determine the precise earthquake source parameters; S3, obtain the tectonic stress field characteristics through focal mechanism solutions and GNSS observation data; S4, predict the earthquake-pregnancy mechanism, seismogenic structure and source rupture process in the study area, combine the velocity structure, fault distribution, focal mechanism, velocity field characteristics and post-earthquake field investigation data to predict the impact of the combined action of faults on the earthquake-pregnancy process, and thus obtain information on potential strong earthquake hazards.

2. The method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data according to claim 1, characterized in that: In step S1, seismic and geological data information of the study area is collected, including: collecting earthquake phase reports, seismic waveform data, GNSS observation data, establishing a basic database, and collecting geological and structural data of earthquake-prone areas in the study area through field surveys and reviewing previous data; In step S2, a high-resolution three-dimensional crustal velocity model of the study area is established to determine the precise earthquake source parameters, including: S201, double-difference seismic tomography tomoDD inversion was performed using a grid interval of 0.2° × 0.2° to obtain the velocity structure MODEL_0.2 of the study area and surrounding areas; S202, linear interpolation of the velocity structure MODEL_0.2 in the study area and surrounding areas was performed to obtain the velocity model MODEL_0.1 with a grid interval of 0.1°×0.1°; S203: Using the velocity model MODEL_0.1 with a grid interval of 0.1°×0.1° as the initial model, the double-difference seismic tomography method tomoDD is used for inversion calculation to obtain the final three-dimensional crustal P-wave velocity structure MODEL_final and the earthquake location results in the study area.

3. The method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data according to claim 1, characterized in that: In step S3, the tectonic stress field characteristics are obtained through the focal mechanism solution and GNSS observation data, including: S301, obtaining focal mechanism solutions of multiple earthquakes; S302, study the regional tectonic stress field based on focal mechanism deconstruction; S303, obtaining the velocity field and the surface expansion rate field through GNSS data; S304, fault plane parameter inversion, uses the double-difference seismic tomography method to obtain accurate source parameter information after relocation. Based on the principle that clustered earthquakes occur near faults, combined with the simulated annealing algorithm and the Gauss-Newton algorithm, the fault plane parameter inversion is performed on the seismic strip near a certain study area to obtain the fault strike, dip, NW dip, fracture length, and fault depth.

4. The method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data according to claim 3, characterized in that: In step S302, based on the focal mechanism deconstruction, the regional tectonic stress field is studied, including: Based on the magnitude, the focal mechanism solutions are weighted and the focal mechanisms of large earthquakes and small earthquakes are combined. The given weights are obtained by the following formula: D=(M-M min ) / ln(ω max ) In the formula, ω is the given weight, e r is the source attenuation coefficient, D is the magnitude attenuation coefficient, ω max is the maximum magnitude weight, M is the current magnitude source, and M min The smallest magnitude earthquake source.

5. The method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data according to claim 3, characterized in that: In step S303, the velocity field and the surface expansion rate field are obtained through GNSS data, including: Step 1: Calculate the mean, randomness, and stability of the velocity field dataset between two adjacent continents in the plate reference frame; calculate the thresholds of the three characteristics based on the Euler parameters and the single-day relaxation solution parameters; Step 2: determine whether it is surface expansion data based on the threshold, remove the background field data, and retain the surface expansion data; Step 3: Use global linear embedding to calculate the velocity change rate in the plate reference frame between two adjacent continents, map the feature matrix of the three-dimensional space to the low-dimensional surface expansion rate field embedding space RD, and perform data dimensionality reduction and parameter fusion on the three feature quantities of each sample; Step 4: Input the fused feature data into the deep belief network (DDN) and use the fault variation optimization method to optimize the input parameters of the DDN. Step 5: During the DDN training process, the fault variation optimization method is used to optimize the parameters of the DDN's visible and hidden layer biases a and b, as well as the number of network layers and nodes. After data training, mathematical models corresponding to the three surface expansions are obtained, completing the construction of a self-verification system for velocity field data in the plate reference frame between two adjacent continents. Step 6: Input the signal of the velocity field data acquisition system into the velocity field data self-verification system. The self-verification system processes the data and inputs the processing results into the mathematical model of the velocity field data self-verification system to complete the surface expansion type discrimination of the velocity field signal acquisition system under the plate reference frame between the two adjacent continents.

6. The method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data according to claim 5, characterized in that: In step 1, the mean, randomness, and stability of the velocity field dataset between two adjacent continents are calculated in the plate reference frame, including: (1) Calculate the mean F1 using the following formula: Where S is the total number of sample points, u i is i sample points; (2) Calculate randomness: Arrange the collected signals from large to small according to their values, merge all data of the same size, and the number of retained data is the data randomness index F2; (3) Calculate the stability according to the following formula: Where, is the stability value of the current n continent sample points of the normal signal, u(n) is the stability value of the total n continent sample points of the normal signal, n is the nth continent, u(η) is the difference between the sample points of two adjacent continents, τ is the previous continent, η is the difference between two adjacent continents, Q(n) is the positive characteristic value of the square root of the stability value of the normal signal, and F3 is the data stability index; Calculate the mean threshold G1, the expression is: Calculate the randomness threshold G2, the expression is: Calculate the stability threshold G3, the expression is: Where, v i is the i velocity sample point, f s Arrange the collected signals from large to small in value, merge all the data of the same size and retain the data. is the stability characteristic of the normal signal.

7. The method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data according to claim 5, characterized in that: In step 3, the velocity change rate in the plate reference frame between two adjacent continents is embedded using global linear embedding, including: mapping the three-dimensional spatial feature matrix to the low-dimensional surface expansion rate field embedding space O E The three feature quantities of each sample are subjected to data dimension reduction and parameter fusion; the specific steps are: (1) Select a local neighborhood. For a given data set, we have: U={U1,U2…U N },U i ∈O E Where U is the dataset given by the local neighborhood, U N is the Nth local neighborhood data, U i is the i-th local neighborhood dataset; Find each sample point U i The k, k<N nearest neighbor points in the neighborhood are calculated using the Euclidean distance formula: Where, 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 in the neighborhood k, u jk is the velocity value of the jth sample point in the kth neighborhood, k is the kth neighborhood, N is the Nth neighborhood, and D 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: Where ε(R) is the weight error, u i is i sample points, u j are j sample points, is the second-order absolute value, R ij for u i with u ij The weight between ij for u i k nearest neighbors of , j = 1, 2, ..., k; To constrain the conditions, the error function expression is transformed into: Where R i is the local reconstruction weight of the i-th sample point, R i =[R1,R2,…,R k ] G , R k is the local reconstruction weight of the kth neighborhood, G is the iteration value; Z i Get the local covariance matrix for the i-th sample point, Z i =(u i -u ij ) G (u i -u ij ); The Lagrange multiplier is introduced to solve the constraint problem, and the expression is: Where λ is the constraint weight, L is the Lagrange multiplier, and W ij Solve the constraint weights for the i-th sample point and the j-th sample; Order Z i R i =1, readjust the weights so that the sum of the weights is 1, and finally get R i ; (3) The low-dimensional surface expansion rate field of the sample set is found through the obtained weight matrix R to embed Z, and the reconstruction error and function are minimized. The expression is: Where φ(Z) is the reconstruction error and function minimization value, z j is the local covariance matrix of j sample points; With restrictions on Z, the expression is: Where I is the N-dimensional identity matrix; The optimization problem is transformed into a constrained optimization problem, which is expressed as: Where, I i is the i-th unit matrix, R i is the local reconstruction weight of the i-th sample point, ZMZ G is the constrained optimization value of G iterations, Z G is the data set after G iterations, Z is the data set; The dataset Z is equivalent to finding the eigenvectors of the symmetric, semi-positive definite, sparse matrix M, and we get: M=(I-W) G (I-W) Where W is the constraint weight matrix to be solved; Using the Lagrange multiplier method, we get: L(Z)=ZMZ G -λ(ZZ G -SI) Where L(Z) is the constraint optimization value obtained by the Lagrange multiplier method, and λ is the constraint weight; Taking the partial derivative of Z to minimize L(Z) yields: Where, MZ G =λZ G , finding Z is equivalent to finding the eigenvector of M, thus obtaining MZ G =λZ G , the embedded coordinates obtained are the eigenvectors of M; the eigenvectors corresponding to the smallest d non-zero eigenvalues ​​are used as the values ​​of M, and the low-dimensional surface expansion rate field coordinates Z are obtained. The eigenvectors corresponding to the eigenvalues ​​are the output results; The entire process of velocity change rate between two adjacent continents in the plate reference frame is expressed as follows: Where, is the velocity set between two adjacent continents in the plate reference frame, is the local reconstruction speed weight set of the i-th sample point, is the velocity change data set after embedding the low-dimensional surface expansion rate field, → is the rate of change.

8. The method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data according to claim 5, characterized in that: 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 the fault variation optimization method, including: (i) Constructing a deep belief neural network (DDN); Based on the energy function formula for the calculation, the joint probability distribution of (v,h) is calculated as: p(v,h|θ)=l -C(v,h|θ) / t(θ) Where C(v,h) is the joint probability distribution value, I is the N-dimensional unit matrix, J is the J samples, e i is the attenuation coefficient of the ith earthquake source, v i is the i velocity sample point, f i is the randomness value of the i sample point, h j is the thickness of j sample points, ξ ji is the energy coefficient between the jth sample point and the ith sample, p(v,h|θ) is the calculated value of the energy function, l -C(v,h|θ) is the energy value calculated under the joint probability distribution, t(θ) is the partition function; Where v is velocity and h is thickness; The likelihood function defined by RBM, p(v|θ) is the marginal distribution of the established joint probability p(v,h|θ); Since neurons in the RBM layer are not connected, the operation state and activation of the hidden layer are relatively independent according to the state of the visible layer. The activation probability of the jth hidden layer node is: Where F'() is the current activation probability function, h j is the thickness of j hidden layers, θ is the likelihood coefficient, σ() is the sigmoid function, f j is the randomness value of j sample points; Where, σ(u)=1 / (1+l -u ) is the sigmoid function; Given the state of the hidden layer nodes, the activation probability of the i-th visible layer node is: Where, e i is the attenuation coefficient of the ith earthquake source; RBM is trained and operated in an iterative manner. The field of the operating surface expansion rate is to study the calculation parameter θ=(ξ ij ,e i ,f j ), given training data and training samples; the maximum log-likelihood function is calculated by parameter θ, the number of training set samples is set to T, and the expression is: Where θ * To calculate the maximum log-likelihood function value through the parameter θ, F'(v (T) |θ) is the activation probability function of the joint probability under the number of training set samples T; (ii) DDN training process: The velocity change rate parameter in the plate reference frame between two adjacent continents is integrated with the velocity field signal data feature in the plate reference frame between two adjacent continents as input, and one RBM is trained at a time from the bottom up. At each layer, the parameter space ωk is constructed by the values ​​calculated in the k-1th layer, and the weights are updated according to the following formula: Where, ξ ij Update the value for the weight, 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 standard value of (v,h).

9. The method for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data according to claim 5, characterized in that: In step 5, the fault change optimization method is used to optimize the parameters of the DDN's visible layer and hidden layer bias a, b, the number of network layers and the number of nodes, including: (a) Initialization of the parameters of the fault change optimization method and the optimization problem. The optimization problem is described as: Minimize f(u)subject to u j ∈E j Where, Minimize f(u) is the initialization value of the optimized surface expansion rate field, u j is the decision vector velocity of j samples, f(u) is the optimized surface expansion rate field function, j is the signal length of j samples, j=1,2,3…N, E j After the decision vector u j The constraint interval of (b) Set the initial plate and calculate the enthalpy or entropy value, and 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 ), n is the number of data points of the initial plate; initialize the partition factor, set k = 1; set the partition factor to 1, if the partition factor k = 2, then generate two more different plates; if the partition factor value is k, then generate 2k-2 different plates; if the initial plate category is R and the population size is F, then just satisfy R < F and continuously add new plates to the initial plate set; stop when R > F, and initialize to obtain an initial group containing F plates; (c) Simulating the surface expansion rate field change process and encoding the surface expansion rate field change process; (d) Plate update: The plate is updated by testing the equilibrium of changes; the state when all surface expansions no longer change is the equilibrium state; in the equilibrium state, the plate that increases the entropy or decreases the enthalpy of surface expansion is considered a new plate, and abnormal plates are excluded; (e) Determine the iteration termination condition: if the condition is met, the algorithm terminates, otherwise it proceeds to step (c); the termination condition of the fault change optimization method is to meet the maximum number of iterations or the minimum enthalpy value and the maximum entropy value.

10. A system for constructing a three-dimensional model of a seismic fault based on aftershock precise positioning data, characterized in that: The system implements the method for constructing a three-dimensional model of a seismogenic fault based on aftershock precise positioning data as described in any one of claims 1 to 9, and the system comprises: Information collection module, used to collect earthquake and geological data information in the study area; The module for establishing a 3D crustal velocity model is used to establish a high-resolution 3D crustal velocity model of the study area and determine accurate earthquake source parameters. Tectonic stress field feature acquisition module, used to obtain tectonic stress field features through focal mechanism solutions and GNSS observation data; The module for predicting earthquake-pregnancy mechanism, earthquake-generating structure and source rupture process in the research area is used to combine velocity structure, fault distribution, source mechanism, velocity field characteristics and post-earthquake field investigation data to predict the impact of the combined action of faults on the earthquake-pregnancy process, thereby obtaining information on potential strong earthquake hazards.

Images