A fault geometry inversion method, an electronic device, and a readable storage medium

Abstract

The present invention discloses a tomographic geometry inversion method, an electronic device and a readable storage medium. The method comprises the following steps: constructing an initial tomographic model: initializing the tomographic depth and n tomographic curves, and representing the initialized parameters by a vector m; obtaining the observation data of the target area; establishing a forward model and an inversion objective function based on the initial tomographic model and the observation data; using a gradient-based method to calculate the inversion objective function, obtaining the optimized vector m when the inversion objective function takes the minimum value, and further obtaining the updated tomographic model to complete the tomographic geometry inversion. The present invention obtains the geometric shape of the fault by inverting the long-term slip rate, has high inversion efficiency, strong practicability, and solves the problem that the existing methods cannot obtain the geometric shape of the fault without co-seismic observation data.

Description

Technical Field

[0001] The present invention relates to the field of geographic information technology, and particularly to a fault geometry inversion method, an electronic device, and a readable storage medium. Background Art

[0002] The geometry of a fault controls the dynamic behavior of fault slip, and is the basis for studies on the physical process of the earthquake source, numerical simulation of strong ground motion, seismic hazard assessment, seismogenic environment research, seismogenic structure research, and inversion of earthquake source dynamics and kinematics. Therefore, determining the geometry of a fault is of great significance for seismological research.

[0003] Currently, there are mainly three methods for obtaining the geometry of an earthquake: precise aftershock location, seismic reflection method, and fault geometry inversion. Precise aftershock location refers to using the seismic wave signals recorded by stations to obtain the early aftershock catalog of an event, and inferring the fault geometry through the spatial distribution of aftershocks. Precise aftershock location usually has low accuracy, sometimes cannot cover the entire fault area well, and is only applicable to earthquake events with abundant aftershocks. The seismic reflection method refers to collecting artificially excited seismic wave signals and recording them by an array seismograph, and obtaining the wave impedance distribution of the reflection profile through a series of seismic data processing methods, including numerical filtering, deconvolution, velocity analysis, static correction, dynamic correction, and migration, etc., so as to identify the spatial position of the fault. Although the seismic reflection method has high accuracy, it also has high costs, and it is difficult to accurately obtain fault information when affected by terrain and lateral medium inhomogeneity during mountain exploration.

[0004] Existing fault geometry inversion methods mainly invert co-seismic observation data. When an earthquake occurs, the ground surface will deform, and this deformation information can be recorded by satellites, including GPS, GNSS, InSAR. At the same time, it is also possible to perform constrained inversion by combining seismic wave records. The relationship between deformation data and fault slip can be represented by the Green's function, and finally a linear equation system for solution is formed, and the least squares method is used for solution.

[0005] However, co-seismic observation data mainly comes from strong earthquakes that have already occurred. However, the information acquisition technology for surface deformation developed in the past 20 years cannot obtain information on fault activities that have occurred or not occurred. Earthquakes that pose a huge safety hazard to human society are usually strong earthquakes above magnitude 7, and the recurrence period of these large earthquakes is usually more than several hundred years. This makes the existing inversion technology unable to obtain the geometric information of most faults in nature, such as earthquakes that occurred (more than 50 years ago) and earthquakes that have not occurred. Summary of the Invention

[0006] In view of the above deficiencies in the prior art, a fault geometry inversion method, an electronic device, and a readable storage medium provided by the present invention solve the problem that the existing fault inversion method cannot invert faults for which co-seismic observation data is not recorded.

[0007] To achieve the above object of the invention, the technical solution adopted by the present invention is as follows:

[0008] Provide a fault geometry inversion method, which includes the following steps:

[0009] S1. Construct an initial fault model: Initialize the fault depth and n fault curves, and represent the initialized parameters with a vector m;

[0010] S2. Obtain the observation data of the target area, including the maximum horizontal principal stress, the minimum horizontal principal stress, the vertical principal stress, and the long-term sliding rate at different positions on the fault;

[0011] S3. Establish a forward model and an inversion objective function based on the initial fault model and the observation data;

[0012] S4. Use a gradient-based method to calculate the inversion objective function, and obtain the optimized vector m when the inversion objective function takes the minimum value, and then obtain the updated fault model to complete the fault geometry inversion.

[0013] Further, the fault curve in step S1 includes at least one of a straight-line curve and an inverse trigonometric function curve; the expression of the straight-line curve is: z(y) = τ * y; the expression of the inverse trigonometric function curve is: z(y) represents the fault depth; τ * is the curvature parameter; y is the horizontal distance from a point on the fault to the surface fault trace; π is the pi.

[0014] Further, the specific method of step S3 includes the following sub-steps:

[0015] S3-1. Discretize the initial fault model into several fault units;

[0016] S3-2. Obtain the strike and dip angle θ of each fault unit through spline interpolation;

[0017] S3-3. According to the formula:

[0018]

[0019]

[0020]

[0021] And respectively obtain the unit vectors along the strike and the unit vector of dip as well as the normal vector

[0022] S3-4. According to the formula:

[0023]

[0024]

[0025] σ = R T τR

[0026] Obtain the horizontal traction force T1 ij (m) of the calculation point on the j-th fault unit and the traction force along the dip of the calculation point on the j-th fault unit where σ is the regional stress tensor; R is the rotation matrix; (·) T represents the transpose of the matrix; τ = diag{τ H , τ h , τ v}; diag{·} is the diag function; τ H is the horizontal maximum principal stress; τ h is the horizontal minimum principal stress; τ v is the vertical principal stress;

[0027] S3-5. Construct the forward model:

[0028] T1(m) = [T1 1 (m), T1 2 (m),..., T1 i (m),...]

[0029]

[0030]

[0031]

[0032]

[0033]

[0034] where T1(m) is the comprehensive traction force along the strike; T2(m) is the comprehensive traction force along the dip; T1 i (m) is the i-th comprehensive traction force along the strike; is the i-th comprehensive traction force along the dip; W j is the weight coefficient of the calculation point on the j-th fault unit; N is the number of weighted points within the calculation radius; is the nth fault curve; h0 is the fault depth; cos(·) is the cosine function; r j is the distance between the observation point and the fault unit; R j is the integration radius; π is the pi;

[0035] S3-6. Construct the inversion objective function:

[0036] ψ(m) = ||αT1(m) - V1|| 2 + ||αT2(m) - V2|| 2 + λ||Lm|| 2

[0037] where ψ(m) is the objective function, and the corresponding vector m is obtained when ψ(m) takes the minimum value; α is a constant; V1 is the long-term horizontal sliding rate of the observation point; V2 is the long-term vertical sliding rate of the observation point; λ is the regularization parameter, λ > 0; L is the first-order difference operator.

[0038] Further, the gradient-based methods in step S4 include the conjugate gradient method and the quasi-Newton method.

[0039] Provide an electronic device, which includes

[0040] a memory storing executable instructions; and

[0041] a processor configured to execute the executable instructions in the memory to implement the fault geometry inversion method.

[0042] Provide a readable storage medium, on which executable instructions are stored, and the fault geometry inversion method is implemented when the executable instructions are executed by a processor.

[0043] The beneficial effects of the present invention are as follows: The present invention obtains the geometric shape of the fault by inverting the long-term sliding rate, with high inversion efficiency and strong practicability, and solves the problem that the existing methods cannot obtain the geometric shape of the fault without co-seismic observation data. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 is the flow schematic diagram of the method;

[0045] Figure 2 is the schematic diagram of calculating the traction force by fault plane interpolation;

[0046] Figure 3 is the schematic diagram of the distribution of long-term sliding rate measurement points of the Wenmao fault in the embodiment;

[0047] Figure 4 is the schematic diagram of the initial fault inversion result; where Figure 4 (a) is from τ *The initial fault plane obtained with the initial values of =5 and h0=18 km; Figure 4 (b) The average τ of 100 inversions * The final fault plane obtained;

[0048] Figure 5 The comparison result between the observed long-term slip rate and the transformed traction force;

[0049] Figure 6 The distribution and average value of the inversion parameter τ for 100 times * ;

[0050] Figure 7 The objective function descent curve for 100 inversions. Specific implementation manners

[0051] The specific implementation manners of the present invention will be described below to facilitate the understanding of the present invention by those skilled in the art. However, it should be clear that the present invention is not limited to the scope of the specific implementation manners. For those skilled in the art, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions created using the concept of the present invention are within the scope of protection.

[0052] As Figure 1 shown, the fault geometry inversion method includes the following steps:

[0053] S1. Construct an initial fault model: Initialize the fault depth and n fault curves, and represent the initialized parameters with a vector m;

[0054] S2. Obtain the observation data of the target area, including the maximum horizontal principal stress, the minimum horizontal principal stress, the vertical principal stress, and the long-term slip rate at different positions on the fault;

[0055] S3. Establish a forward model and an inversion objective function based on the initial fault model and the observation data;

[0056] S4. Use a gradient-based method to calculate the inversion objective function. When the inversion objective function takes the minimum value, obtain the optimized vector m, and then obtain the updated fault model to complete the fault geometry inversion. The gradient-based methods include the conjugate gradient method and the quasi-Newton method, and the conjugate gradient method is preferably used.

[0057] In step S1, the fault curve includes at least one of a linear curve and an inverse trigonometric function curve; the expression of the linear curve is: z(y) = τ * y; the expression of the inverse trigonometric function curve is: z(y) represents the fault depth; τ *is the curvature parameter; y is the horizontal distance from a point on the fault to the surface fault track; and π is the circumference of a circle. In the inversion, all fault curves have the same fault depth, so a total of n+1 parameters need to be inverted: n curvature parameters and one fault depth.

[0058] The long-term slip rate at the observation point can be approximated as the combined effect of stress on the fault within a circle with a radius of R. Therefore, by calculating the weighted average stress across multiple fault units, the regional stress field can be projected onto the fault surface, yielding the combined shear stress near the surface. This allows the establishment of a forward model linking fault geometry with shear stress near the surface. Specifically, the method for step S3 includes the following sub-steps:

[0059] S3-1. Discrete the initial fault model into several fault units;

[0060] S3-2. Obtain the strike direction of each fault unit using spline interpolation and the inclination angle θ; the cubic spline interpolation method can be used in specific implementation;

[0061] S3-3, according to the formula:

[0062]

[0063]

[0064]

[0065] And get the unit vector along the direction respectively and the unit vector of the inclination angle and the normal vector

[0066] S3-4, according to the formula:

[0067]

[0068]

[0069] σ=R T τR

[0070] Get the horizontal traction T1 of the calculation point on the jth fault unit ij (m) and the traction force along the dip angle at the calculation point on the jth fault unit Where σ is the regional stress tensor; R is the rotation matrix; (·) T represents the transpose of the matrix; τ=diag{τ H ,τ h ,τ v}; diag{·} is the diag function; τH is the horizontal maximum principal stress; τ h is the horizontal minimum principal stress; τ v is the vertical principal stress;

[0071] S3-5. Construct a forward model:

[0072] T1(m) = [T1 1 (m), T1 2 (m),..., T1 i (m),...]

[0073]

[0074]

[0075]

[0076]

[0077]

[0078] where T1(m) is the comprehensive traction force along the strike; T2(m) is the comprehensive traction force along the dip angle; T1 i (m) is the i-th comprehensive traction force along the strike; is the i-th comprehensive traction force along the dip angle; W j is the weight coefficient of the calculation point on the j-th fault unit; N is the number of weighted points within the calculation radius; is the n-th fault curve; h0 is the fault depth; cos(·) is the cosine function; r j is the distance between the observation point and the fault unit; R j is the integration radius; π is the pi;

[0079] S3-6. Construct an inversion objective function:

[0080] ψ(m) = ||αT1(m) - V1|| 2 + ||αT2(m) - V2|| 2 + λ||Lm|| 2

[0081] where ψ(m) is the objective function, and the corresponding vector m is obtained when ψ(m) takes the minimum value; α is a constant; V1 is the long-term horizontal sliding rate of the observation point; V2 is the long-term vertical sliding rate of the observation point; λ is the regularization parameter, λ > 0; L is the first-order difference operator. When performing the first iteration, λ can be set to 1, and the value of λ is halved in each iteration process. This can make the forward model move as a whole in the early stage and change locally in the later stage to approach the optimal solution.

[0082] In the specific implementation process, the schematic diagram of interpolating and calculating the traction force on the fault plane is as Figure 2 shown. The circles are equally spaced control points on the fault curve, the squares represent equally spaced interpolation points on the curve, and the inverted triangles represent the positions of the observation points. Figure 2 In which, X is the horizontal interval of the fault curve; fault trace represents the position of the fault on the ground surface; interpolation curve represents the fault interpolation curve; dip curve represents the fault curve; the i-th measured point represents the i-th observation point; T ij represents the j-th traction force used for weighted calculation at the i-th observation point.

[0083] In an embodiment of the present invention, the geometry of the Wenchuan - Maoxian fault (abbreviated as the Wenmao fault) is used to test this method. In the past decade, there have been many studies on the fault geometry of the Beichuan fault and the Pengguan fault, including aftershock location, seismic reflection wave exploration, and joint inversion of GPS and InSAR data, etc. However, due to the rugged terrain, the distribution of aftershocks is very limited, and the fine geometry of the Wenmao fault is still unclear so far.

[0084] As Figure 3 shown, in order to obtain the geometry of the Wenmao fault, the long-term slip rates along the Wenmao fault were collected from previous studies (Ma et al., 2005; Rongjun et al., 2007; Ran et al., 2013; Shen et al., 2019). A total of 5 measurement points and 7 long-term slip rates were inverted. Although the distribution of the observation points does not cover the entire fault area well, it is effective in constraining the first-order geometric characteristics of the Wenmao fault. Figure 3 Each rectangular box in it contains the long-term slip rates from different literature researches. Dextral refers to the dextral strike-slip rate, Up(NW) represents the vertical slip rate, where the northwestern block is the hanging wall. The solid line represents the fault trace on the ground surface. The fault near the northwest is the Wenmao fault, and the dashed line represents the Lixian fault. The accurate fault trace of the Wenmao fault comes from the field investigation by Xie Xinsheng et al. (2011). The five-pointed stars represent the positions of the measurement points of the long-term slip rate. The circles represent the aftershock distribution of the 2008 earthquake (Yin et al., 2018).

[0085] Due to the uncertainty in estimating the geological age of sediments, the long-term slip rate cannot be accurately measured. Therefore, 100 bootstrap resamplings of the entire data set are used to estimate the uncertainty of the inversion, and it is assumed that the data from each measurement point follows a Gaussian distribution. For each inversion, as Figure 4 (a) shown, the same initial model is set: τ *= 5, h0 = 18 km; Figure 4 The x - coordinate represents the distance along the strike, the y - coordinate represents the distance orthogonal to the fault strike, and the z - direction represents the depth. In addition, some aftershocks are distributed at the southern end of the Wenmao fault, which can be used as prior information to constrain the fault geometry curve at the southernmost end. Through trial and error, setting τ * = 3.5 and h0 = 18 km can better fit the distribution of aftershocks. Therefore, by modifying the weights of the first - order difference operator L, the parameter τ * at the southernmost end (x = 0 km) can be kept with small changes during the inversion process. As Figure 4 (b) shows, the average value of the inversion parameter τ * for 100 times is used to construct the final fault model. Figure 5 In it, the abscissa represents the distance along the strike, and the ordinate represents the long - term slip rate or the normalized traction force. From Figure 5 it can be seen that due to the constraint of the long - term slip rate, the dip angle near the surface changes significantly along the fault strike. As Figure 6 and Figure 7 show, most of the inversion results of τ * are distributed in a narrow area, and all objective functions decrease steadily with the increase of the number of iterations, indicating the effectiveness of this inversion method. Figure 6 In it, the triangles represent the positions of the interpolation curves, the abscissa is the distance along the strike, and the ordinate is the fault geometry parameter; Figure 7 In it, the abscissa is the number of iterations, and the ordinate is the root - mean - square.

[0086] The accurate determination of the fault structure, especially the dip angle near the free surface, provides key information about the mechanical behavior of the fault system and earthquake rupture. The results of this method show that the dip angles in the shallow part of the Wenmao fault are generally larger than those of the Beichuan fault and the Pengguan fault. The dip angle range near the surface of the Wenmao fault is between 50° and 70°. While Wan et al. (2017) used multiple datasets to invert the fault geometry and slip distribution of the Wenchuan earthquake and found that the near - surface dip angle at the southwestern end of the Beichuan fault is 36°, and the northeasternmost end is close to vertical (dip angle 83°). In addition, the dip angles of the shallow part of the Beichuan fault in Hongkou, Qingping, and Beichuan are 43°, 67°, and 51° respectively, and the dip angle of the Pengguan fault is about 35°, and it is connected to the Beichuan fault zone at a depth of 11 km. The conclusions obtained by the present invention are consistent with the previous conclusions (Hubbard et al., 2010; Jia et al., 2010), indicating the correctness of the inversion results obtained by the present invention.

[0087] In summary, the forward model constructed in the present invention obtains the comprehensive traction force at the surface observation point by calculating the weighted average stress on multiple fault units, thereby relating the fault geometric parameters to the comprehensive traction force near the surface. The inversion objective function establishes the connection between the long-term slip rate direction and the shear stress direction, and the loss function is established by minimizing the difference between the ratios. At the same time, a regularization scheme is introduced, which can inversely obtain a smooth fault model, obtain the geometric shape of the fault, with high inversion efficiency, strong practicability, and solves the problem that the existing methods cannot obtain the geometric shape of the fault without co-seismic observation data.

Claims

1. A fault geometry inversion method, characterized in that, It includes the following steps: S1. Construct an initial fault model: Initialize the fault depth and n fault curves, and represent the initialized parameters with a vector m for representation; S2. Obtain the observation data of the target area, including the horizontal maximum principal stress, the horizontal minimum principal stress, the vertical principal stress, and the long-term sliding rate at different positions on the fault; S3. Establish a forward model and an inversion objective function based on the initial fault model and the observation data; S4. Use a gradient-based method to calculate the inversion objective function, and obtain the optimized vector when the inversion objective function reaches the minimum value, m and then obtain the updated fault model to complete the fault geometry inversion. The specific method of step S3 includes the following sub-steps: S3-1. Discretize the initial fault model into a number of fault units; S3-2. Obtain the strike of each fault unit by spline interpolation method and dip angle ; S3-3. According to the formula: and respectively obtain the unit vectors along the strike and the unit vectors of dip angle , as well as the normal vector ; S3-4. According to the formula: Obtain the horizontal traction force of the calculation point on the j th fault unit and the dip traction force of the calculation point on the j th fault unit ; where is the regional stress tensor; is the rotation matrix; represents the transpose of the matrix; ; is the diag function; is the maximum horizontal principal stress; is the minimum horizontal principal stress; is the vertical principal stress; S3-5. Construct a forward model: in is the combined traction along the strike; is the integrated traction along the inclination; For the i The combined traction along the strike direction; For the i The combined traction along the inclination angle; For the j The weight coefficient of the calculation point on each fault unit; N is the number of weight points within the calculation radius; For the n The curvature parameters of the fault curves; is the fault depth; is the cosine function; The observation point and j The distance between fault units; For the j The integral radius of a fault unit; π is the circumference of a circle; S3-6. Construct an inversion objective function: where is the objective function, and the corresponding vector is obtained when it takes the minimum value; m ; is a constant; is the long-term horizontal sliding rate of the observation point; is the long-term vertical sliding rate of the observation point; is the regularization parameter, ; is the first-order difference operator.

2. The fault geometry inversion method according to claim 1, characterized in that The fault layer curve in step S1 includes at least one of a linear curve and an inverse trigonometric function curve; the expression of the linear curve is: ; The expression of the inverse trigonometric function curve is as follows: ; represents the fault depth; is the curvature parameter; y is the horizontal distance from a certain point on the fault to the surface fault trace; π is the pi.

3. The fault geometry inversion method according to claim 1, characterized in that In step S4, the gradient-based methods include the conjugate gradient method and the quasi-Newton method.

4. An electronic device, characterized in that, Include A memory storing executable instructions; and A processor configured to execute the executable instructions in the memory to implement the fault geometry inversion method according to any one of claims 1 to 3.

5. A readable storage medium having executable instructions stored thereon, characterized in that, When the executable instructions are executed by the processor, the fault geometry inversion method according to any one of claims 1 to 3 is implemented.

Images