Fault bending quantification method based on high-resolution surface fracture

Abstract

The invention discloses a fault bending quantification method based on high-resolution surface fracture, and belongs to the technical field of earthquake monitoring and early warning, and the method comprises the steps: precisely extracting a surface fracture track through fusing InSAR, POT and field survey data; constructing a three-dimensional dislocation model in combination with Delaunay triangulation and slip inversion; and calculating the roughness by using the included angle between the adjacent triangular units, identifying the bending center and the bending interval of the fault, and projecting the bending center and the bending interval to the earth surface to realize the spatial representation and quantitative characterization of the bending characteristics of the fault. According to the method, the problems of insufficient deformation information and inaccurate fault geometric parameter extraction in a traditional method are solved, the fault discontinuity identification precision can be effectively improved, key constraints are provided for research of an earthquake fracture mechanism, and the method has relatively high innovativeness and engineering application value.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of earthquake monitoring and early warning, and particularly relates to a fault bending quantification method based on high-resolution surface rupture. BACKGROUND

[0002] In earthquake research, it is of great significance to deeply explore the complex geometric structure of faults. The complex geometric characteristics of faults not only affect the initiation, propagation and termination of earthquakes, but also control the magnitude, frequency and rupture pattern of earthquakes. By studying the complex geometric shape of faults, the earthquake potential and disaster risk can be more accurately assessed, providing a scientific basis for earthquake prediction and disaster prevention. In addition, the study of fault geometric structure helps to reveal the mechanism of earthquake preparation and occurrence, improve the accuracy of earthquake monitoring and early warning, and optimize the earthquake early warning system. At the same time, these studies can also provide important clues for the formation and evolution of geological structures, and provide important geological information for engineering site selection and resource exploration.

[0003] During the development and maturation of faults, the geometric structure of faults will undergo a series of changes, which have important influence on the preparation and occurrence of earthquakes. The geometric structure of faults mainly includes fault bifurcation, step zone, bending, etc. Among them, the fault bending structure is a form of particular concern in earthquake research, which is characterized by the bending of the fault on the original rupture trajectory, forming a curved or zigzag fault plane. Previous studies have shown that this structure is often the starting point (seismic source) of earthquake rupture or an obstacle to rupture propagation, thereby affecting the rupture process and magnitude of earthquakes. In addition, fault bending structure may also control the spatial distribution of aftershocks, as stress adjustment at the bending site can trigger a series of secondary ruptures. Therefore, studying fault bending structure is of great significance for understanding the mechanism of earthquake rupture, predicting earthquake hazards, and evaluating fault activity.

[0004] Due to the limitation of obtaining accurate surface rupture trajectories, previous studies on the determination and quantification of fault bending have the following limitations: First, most studies focus on large-scale features (such as large-scale arc structures in the rupture trajectory), which are difficult to effectively capture the details of small-scale bending. Second, the quantification method usually relies on a single observation parameter (such as bending width or angle), lacking a multi-dimensional comprehensive representation of the bending shape. Finally, although the accurate quantification of small-scale fault bending relies on laboratory simulation or numerical simulation, these methods are difficult to establish an accurate correspondence with actual observation data of natural earthquakes, thereby limiting the universality and practical application value of research results. Therefore, in order to accurately identify and quantify the bending structure in the process of earthquake rupture, it is urgent to develop a method that accurately describes the fault shape and establishes a quantitative index of fault shape, to solve the problems of accurate description of fault shape, comprehensive evaluation of fault bending degree, and the need for different fault bending scales in research. SUMMARY

[0005] The present application aims to solve the technical problems in the background art, and aims to provide a fault bending quantification method based on high-resolution surface rupture, based on high-resolution surface rupture observation data, a triangular dislocation slip model is constructed, and the angle change of the slip model along the rupture trajectory is used to accurately depict the fault geometric shape; by extracting the maximum value of the angle distribution of the slip model, the bending position is determined, and it is spatially mapped to the real rupture trajectory, combined with the measurement result of the bending step height, a quantification index of the bending degree is established. By setting reasonable geometric parameter threshold, the required fault bending scale can be effectively identified and quantified.

[0006] In order to solve the technical problems, the technical scheme of the present application is:

[0007] A fault bending quantification method based on high-resolution surface rupture, the method comprises:

[0008] S1: DInSAR or TSInSAR technology is used on SAR image to extract regional scale surface deformation field; for the area where InSAR is out of phase or the deformation gradient is too large, POT technology is further used to process optical image or SAR image to obtain high-precision displacement information; finally, combined with remote sensing image interpretation result and field investigation data, the surface rupture trajectory line is accurately extracted;

[0009] S2: based on the surface rupture trajectory geometry obtained in S1, a three-dimensional dislocation fault initial model is constructed by using Delaunay triangulation algorithm; taking the surface deformation data obtained by InSAR or GNSS as a constraint, the fault slip distribution is inverted based on the least square principle, and the fitting degree is evaluated by residual analysis; the fault geometric parameters are adjusted by grid search method, and the optimal slip model is finally determined, and the fault geometric structure, slip distribution, residual distribution and fitting degree index are output;

[0010] S3: using the triangular dislocation unit set in the slip model, the average angle between each unit and its adjacent unit is calculated as the roughness of the unit; the roughness in the direction of inclination is counted, and the roughness distribution of the model in the direction of strike is obtained; the roughness-strike position curve is drawn, the extreme point is identified as the bending center, the roughness threshold is set to delimit the bending range, and the roughness value of each unit, the bending center and its spatial position and the bending interval range are output;

[0011] S4: based on the roughness distribution and bending center position obtained in S3, combined with the surface rupture trajectory extracted in S1 and the three-dimensional slip model constructed in S2, the bending center in the fault model is spatially projected to the surface; the vertical distance between the bending center and the bottom of the fault is measured to quantify the bending strength of the fault; the corresponding relationship graph of the bending center and the surface rupture, the bending degree value of each bending center and the spatial distribution graph of the fault bending are output.

[0012] Further, the step S2 discretizes the fault surface into a triangular element grid according to the ground trace and the fault depth, meets the Delaunay criterion, and determines the grid size according to the complexity of the fault rupture, and inverses the fault slip distribution based on the least square principle with InSAR or GNSS ground deformation data as constraints, and specifically includes:

[0013] ①establishing an observation equation: D=G×S, wherein D is an InSAR / GNSS deformation observation vector, G is a Green function matrix, and S is a slip vector;

[0014] ②introducing a smoothing constraint matrix, constructing a slip model to inverse the slip distribution (i.e. solving S) based on the least square criterion, and the observation equation after introducing the smoothing factor is:

[0015]

[0016] In the formula, γ is a smoothing constraint matrix;

[0017] ③solving the parameter D in ② according to the fault geometric parameters and the slip distribution, i.e. forward deformation field, calculating the model residual, and evaluating the model fitting degree. The residual calculation method is:

[0018]

[0019] In the formula, D obs is the deformation value obtained by GNSS / InSAR, D model is the deformation value obtained by model forward, and N is the number of deformation points;

[0020] ④combining the grid search method to adjust the fault geometric parameters, calculating the model residual again, and determining the fault geometric structure and slip distribution corresponding to the minimum model residual as the optimal model by comparing the model residuals under different parameters.

[0021] Further, the step S3 includes:

[0022] Firstly, the roughness of each triangular dislocation is calculated, assuming that the angles of patch2, patch3 and patch4 with patch1 are θ1, θ2 and θ3 (0-90°) respectively, and the roughness θ calculation method follows the formula:

[0023] (n=2 when the triangular dislocation block is located at the edge of the model, otherwise n=3);

[0024] Then, the cumulative value of the triangular dislocation roughness in the dip direction is calculated as the model roughness in the strike direction;

[0025] Finally, the roughness curve of the model along the model trend is drawn, the extreme point of the angle is set as the bending center, the bending width threshold is set according to the actual research needs, and the bending position of the model is determined.

[0026] Compared with the prior art, the present application has the advantages that:

[0027] Based on high-resolution ground observation data and a triangular dislocation model capable of accurately depicting the complex geometry of a fault, the present application realizes accurate quantitative description of the fault morphology.

[0028] Based on the slip model, the bending position is determined, the surface rupture is considered, the deep fault structure is considered, and the morphology determination considering the three-dimensional structure of the fault is realized.

[0029] In the actual research process, according to the specific research needs, by setting appropriate threshold values, different scale bending structures can be accurately extracted, thereby meeting the diversified analysis requirements of the fault bending. BRIEF DESCRIPTION OF DRAWINGS

[0030] Figure 1 The method technical roadmap of the present application (a is the fault bending degree calculation process, b is the determination of the fault bending schematic diagram in the model dimension, and c is the quantitative fault bending degree schematic diagram according to the surface rupture trajectory) ;

[0031] Figure 2 The surface rupture trajectory of the source earthquake interpreted from multiple sources;

[0032] Figure 3 The source earthquake fault slip model and residual of the source earthquake;

[0033] Figure 4 The roughness curve of the model roughness along the fault trend;

[0034] Figure 5 The quantitative value of the fault bending calculated according to the surface rupture trajectory. DETAILED DESCRIPTION

[0035] The specific embodiments of the present application will be described below in conjunction with the embodiments:

[0036] It should be noted that the structures, proportions, sizes, etc. shown in the present specification are only used to cooperate with the content disclosed in the specification, so that those skilled in the art can understand and read, and are not used to limit the implementation conditions of the present application. Any modification of the structure, change of the proportion relationship or adjustment of the size, which does not affect the effect and purpose that can be achieved by the present application, should still fall within the scope of the technical content disclosed by the present application.

[0037] Meanwhile, the terms such as "upper", "lower", "left", "right", "middle", and "one" cited in the specification are only for the convenience of clear description, and are not intended to limit the scope of the application, and the change or adjustment of the relative relationship is also considered as the implementation of the application without substantial changes in the technical content.

[0038] Key terms and abbreviations:

[0039] SAR: Synthetic Aperture Radar

[0040] DInSAR: Differential Interferometric Synthetic Aperture Radar

[0041] TSInSAR: Time Series Interferometric Synthetic Aperture Radar

[0042] POT: Pixel Offset Tracking

[0043] Example 1:

[0044] As shown in Figure 1 , a precise quantification method of fault curvature based on high-resolution surface rupture trace and slip model, specifically comprising:

[0045] The first step is to obtain high spatial resolution surface rupture trace. ① First, process the SAR image through DInSAR or TSInSAR technology to obtain the surface deformation field. ② Then, in the case of large surface deformation, process the optical image or SAR image through POT technology to obtain the surface deformation and solve the problem of out-of-phase coherence caused by large deformation gradient in InSAR technology. ③ Finally, interpret the surface rupture trace through the unmanned aerial image obtained by field investigation and field measurement data.

[0046] The second step is to establish a triangular dislocation slip model. ① First, based on the geometric features of the surface rupture trace, the Delaunay triangulation algorithm is used to construct the initial model of the three-dimensional dislocation fault, which contains the basic elements of fault surface geometric parameters (strike, dip angle, buried depth) and dislocation distribution. ② Take the surface deformation observation data obtained by InSAR or GNSS as the constraint condition, and based on the least square principle, the fault slip distribution is inverted, and the model fitting degree is evaluated by residual analysis. ③ Finally, adjust the fault geometric parameters through grid search method, and when the fitting error no longer significantly reduces, determine the optimal geometric structure of the fault.

[0047] The third step is to calculate the model roughness. ① First, calculate the roughness of each triangular dislocation. Take patch1 in the following figure b2 as an example, and assume that the angles between patch2, patch3 and patch4 and patch1 are θ1, θ2 and θ3 (0-90°) respectively. The roughness θ is calculated according to the following formula.

[0048] (n=2 when the triangular dislocation block is located at the edge of the model, otherwise n=3)

[0049] ② Then, calculate the cumulative value of the triangular dislocation roughness of the model in the dip direction as the model roughness in the strike direction. ③ Finally, draw the roughness curve of the model along the strike direction of the model, and set the extreme point of the angle as the bending center. In addition, according to the actual research needs, set the bending width threshold to determine the bending position of the model (see the following figure b3).

[0050] The fourth step is to quantify the degree of fault bending on the high-resolution surface rupture. ① First, based on the calculation results of the model roughness, the fault bending feature parameters (including the bending geometric center and position) are spatially mapped to realize the projection of deep fault structure to the surface rupture zone; ② Finally, by measuring the distance between the bending center point and the bottom of the bending, the fault bending degree is quantitatively represented, which can effectively represent the strain accumulation characteristics of the geometric discontinuity in the strike-slip fault system, and provide key constraints for subsequent deformation and rupture mechanism analysis.

[0051] Example 2:

[0052] Take the bending structure of the bidirectional rupture seismogenic fault of the 2022 Menyuan earthquake as an example; Figure 2 The coseismic surface rupture trajectory of the Menyuan earthquake based on the comprehensive interpretation of Sentinel-1, GF-7 and unmanned aerial vehicle images. Figure 3 The triangular dislocation slip model of the Menyuan earthquake based on the inversion of surface deformation data. Figure 4 The model roughness curve along the strike direction of the seismogenic fault of the Menyuan earthquake and the determined bending center and position. Through the quantitative analysis of the roughness parameters of the fault slip distribution model, Figure 4 ), it is found that the roughness value along the strike direction of the seismogenic fault mainly distributes in the range of 0-2°. In spatial distribution, the model roughness curve presents four groups of peak values: 3.3°, 3.9°, 4.5° and 3.1° appear in the intervals of-9-7km, -4-2km, 5-7km and 11-13km along the strike direction of the fault respectively Figure 4The bending angle and bending width were used to constrain the bending characteristics of the fault. In this study, the bending center was first determined based on the peak screening (threshold > 2°), and then the bending width (the width of the region where 90% of the peaks are located) was introduced, and finally four regions were determined as the bending of the fault.

[0053] Figure 5 To map the fault bending space determined by the model to the actual surface rupture, and thus measure the fault bending degree result. From the results, this fault bending quantification method shows unique advantages in quantifying the fault morphology with complex geometric structure. Specifically, it can accurately identify small-scale bending characteristics, providing more detailed and accurate data support for related research. In addition, this method avoids the identification of large-scale arc structures with slow curvature changes, effectively capturing the morphological characteristics of the fault rupture trajectory.

[0054] Those skilled in the art will appreciate that embodiments of the present application can be provided as methods, systems, or computer program products. Accordingly, the present application can take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present application can take the form of a computer program product implemented on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROMs, optical storage devices, etc.) containing computer usable program code.

[0055] The present application is described with reference to flowcharts and / or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the present application. It should be understood that each flow and / or block in the flowcharts and / or block diagrams, as well as combinations of flows and / or blocks in the flowcharts and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing apparatus to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing apparatus produce a means for implementing the functions specified in the flowcharts and / or block diagrams. Figure 1 The functions specified in a flow or multiple flows and / or blocks Figure 1 The means for implementing the functions specified in a flow or multiple flows and / or blocks.

[0056] These computer program instructions can also be stored in a computer-readable memory that can direct the computer or other programmable data processing apparatus to work in a specific manner, so that the instructions stored in the computer-readable memory produce a product including instruction means, which implements the functions specified in the flowcharts and / or block diagrams. Figure 1 The functions specified in a flow or multiple flows and / or blocks Figure 1 The means for implementing the functions specified in a flow or multiple flows and / or blocks.

[0057] These computer program instructions can also be loaded into a computer or other programmable data processing devices, so that a series of operational steps are generated to realize the computer-implemented processes, and the instructions executed on the computer or other programmable devices provide a process for implementing the functions specified in the flowchart Figure 1 one flow or multiple flows and / or the functions specified in the block Figure 1 one flow or multiple flows and / or the functions specified in the block

[0058] The preferred embodiments of the present application have been described in detail above, but the present application is not limited to the above-described embodiments, and various changes can be made within the knowledge of those skilled in the art without departing from the spirit of the present application.

[0059] Many other changes and modifications can be made within the scope and spirit of the application without departing from the concepts and scope of the application. It should be understood that the application is not limited to a particular embodiment, and the scope of the application is defined by the appended claims.

Claims

1. A method for quantifying fault curvature based on high-resolution surface rupture, characterized in that: The method comprises: S1: Using DInSAR or TSInSAR technology on SAR images, regional-scale surface deformation fields are extracted. For areas where InSAR decoherence or deformation gradients are excessive, optical or SAR images are further processed using POT technology to obtain high-precision displacement information. Finally, combining remote sensing image interpretation results with field survey data, surface rupture trajectories are accurately extracted. S2: Based on the surface rupture trajectory geometry obtained in S1, an initial 3D dislocation fault model is constructed using the Delaunay triangulation algorithm. Using surface deformation data obtained by InSAR or GNSS as constraints, the fault slip distribution is inverted based on the least squares principle, and the goodness of fit is evaluated through residual analysis. The fault geometry parameters are adjusted using a grid search method to ultimately determine the optimal slip model, outputting the fault geometry, slip distribution, residual distribution, and goodness of fit metrics. S3: Using the set of triangular dislocation units in the slip model, calculate the mean angle between each unit and its adjacent units as the roughness of the unit; statistically analyze the roughness along the dip direction and summarize the roughness distribution of the model in the strike direction; draw a roughness-strike position curve, identify the extreme point as the bending center, set the roughness threshold to define the bending range, and output the roughness value, bending center, spatial position and bending range of each unit; S4: Based on the roughness distribution and bending center position obtained by S3, combined with the surface rupture trajectory extracted by S1 and the 3D slip model constructed by S2, the bending center space in the fault model is projected onto the surface; the vertical distance between the bending center and the bottom of the fault is measured to quantify the fault bending strength; and the corresponding relationship diagram between the bending center and the surface rupture, the curvature value of each bending center, and the spatial distribution map of the fault bending are output.

2. The method for quantifying fault curvature based on high-resolution surface rupture according to claim 1, characterized in that: Step S2 discretizes the fault plane into triangular unit grids according to the surface trace and fault depth, satisfying the Delaunay criterion. The grid size is determined according to the complexity of the fault rupture. The fault slip distribution is inverted based on the least squares principle with InSAR or GNSS surface deformation data as a constraint, specifically including: ① Establish the observation equation: D = G × S, where D is the InSAR / GNSS deformation observation vector, G is the Green's function matrix, and S is the slip vector; ②Introduce the smooth constraint matrix, build the slip model based on the least squares criterion to invert the slip distribution (i.e. solve S). The observation equation after introducing the smoothing factor is: Where γ is the smooth constraint matrix; ③ According to the fault geometry and slip distribution, solve the parameter D in ②, that is, the forward deformation field, calculate the model residual, and evaluate the model fit. The residual calculation method is: Where D obs is the deformation value obtained by GNSS / InSAR, D model is the deformation value obtained by forward modeling, and N is the number of deformation points; ④ Combined with the grid search method, the fault geometric parameters are adjusted and the model residuals are calculated again. By comparing the model residuals under different parameters, the fault geometric structure and slip distribution corresponding to the minimum model residual are determined as the optimal model.

3. The method for quantifying fault curvature based on high-resolution surface rupture according to claim 1, characterized in that: The step S3 comprises: First, the roughness of each triangular dislocation is calculated. Assuming that the angles between patch2, patch3, and patch4 and patch1 are θ1, θ2, and θ3 (0 to 90°), the roughness θ is calculated using the following formula: (n = 2 when the triangular dislocation block is at the edge of the model, otherwise n = 3); Then, the cumulative value of the model's triangular dislocation roughness in the dip direction is calculated as the model's roughness in the strike direction; Finally, the roughness curve of the model is drawn along the direction of the model, the extreme point of the angle is set as the bending center, and the bending width threshold is set according to actual research needs to determine the bending position of the model.

Images