Surface deformation simulation and three-dimensional fault inversion method, system, medium and product

Abstract

The application provides a surface deformation simulation and three-dimensional fault inversion method, system, medium and product, which comprises the following steps: obtaining InSAR data, constructing a preliminary geometric shape of a fault, using relocated aftershock data to constrain the dip angle change of the fault, performing local plane fitting and surface fusion processing to form a continuous and smooth three-dimensional non-planar fault surface; discretizing the three-dimensional non-planar fault surface into triangular dislocation units and performing adaptive encryption, and further subdividing into multiple rectangular subunits; for each rectangular subunit, calculating the displacement contribution by using the Okada analytical kernel to construct a Green function matrix; integrating the InSAR data and GNSS data, combining the Green function matrix, and through regularization constraint and boundary condition control, the non-uniform slip distribution of the fault is obtained by inversion; based on the non-uniform slip distribution result of the fault, combined with the mixed dislocation model, the surface deformation simulation is completed. The application has high precision and stability.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of geophysical and geodetic comprehensive inversion, and particularly relates to a key method and technology for numerical simulation of coseismic deformation and fault slip inversion of an earthquake by using multi-source observation data such as InSAR, GNSS and the like. BACKGROUND

[0002] In the continental active tectonic belt, strong earthquake ruptures often occur on strike-slip or reverse fault systems with significant bending, segmentation and irregular geometry. In order to understand the earthquake rupture process and evaluate the stress transfer effect on the surrounding tectonic units, the academic circles generally rely on global navigation satellite system (GNSS), synthetic aperture radar interferometry (InSAR) and strong earthquake records to invert the fault geometry and slip distribution, and on this basis, analyze the Coulomb stress change and seismic risk. The elastic dislocation theory provides a basic framework for the above problems, and the analytical solution of rectangular plane dislocation given by Okada (1985, 1992) is widely used in surface displacement forward modeling and slip inversion.

[0003] However, in specific applications, in order to simplify the calculation and parameterization, the existing researches usually simplify the real fault geometry into a single-angle plane model, or use several plane rectangular sub-faults for segmentation and splicing, and disperse the fault into a regular rectangular grid for slip inversion. Such models based on rectangular dislocation elements (RDE) assume that the fault plane is planar or linearly spliced by several planes in three-dimensional space in terms of geometry, and have limited ability to describe the curvature changes along the strike and depth directions. Since the simple geometry cannot reflect the bending, step zone and roughness characteristics of the real fault near the surface and in the deep part, artificial complexity is often introduced in the slip distribution during the inversion process to fit the observed deformation: on the one hand, it may cause systematic bias in the estimation of seismic moment, and the position and size of the strong slip area are smoothed or lengthened unreasonably; on the other hand, it may cause numerical artifacts at the segmented connection, and the local strain concentration mode and subsequent Coulomb stress change result are distorted, thereby affecting the physical interpretation of the rupture segmentation, stress transfer path and risk of adjacent faults. The existing researches on several continental strong earthquakes have shown that the single plane or simple multi-segment rectangular fault model has obvious limitations in fitting the coseismic InSAR / GNSS deformation, surface rupture pattern and spatial distribution of aftershocks, especially in the vicinity of the geometric bending or segmented transition zone, and it is difficult to meet the geodetic constraints and seismological constraints at the same time.

[0004] To improve the ability of geometric characterization, some studies introduce the idea of non-planar fault modeling, trying to use multi-segmented polyline faults, layered dip models or irregular meshes constructed by finite elements / boundary elements to approximate real curved surface faults. To some extent, this kind of method can reflect the dip angle changes of the fault along the strike and depth direction, but often constructs three-dimensional geometry in the form of local plane splicing, and the segmented junctions are prone to gaps, overlaps or normal mutations, leading to unstable numerical response of stress and displacement field near the geometric connection zone. In addition, although the finite element / finite difference method can describe complex geometry on volumetric meshes, the cost of meshing and solving is high, and it is not realistic to construct large-scale Green function matrix for inversion at high resolution, so it is more used for forward simulation of a few schemes, and less directly used for large-scale slip inversion driven by multi-source observation data.

[0005] For more detailed description of complex fault geometry, triangular dislocation elements (TDE) have attracted attention in recent years. TDE can naturally fit the curved surface structure of strong bending, multi-step and rough boundary by discretizing the fault surface into non-overlapping and gap-free triangular elements, better avoiding the geometric discontinuity problem at the segmented connection of rectangular sub-faults, and effectively reducing the numerical error caused by segmented boundaries. Existing studies (such as Maerten et al., 2005; Nikkhoo and Walter, 2015; Meade, 2007; Jiang et al., 2013; Zhu and Shan, 2020) have shown that geometric modeling based on triangular dislocation elements is significantly better than simple rectangular plane models in fitting complex deformation data. However, the traditional TDE method often relies on numerical integration of dislocation Green's function, and the computational cost is huge when the number of elements is large, making it difficult to efficiently construct the Green's function matrix in the inversion scenario of jointing a large amount of InSAR / GNSS data. Therefore, most of the existing TDE-related work stays at the level of forward simulation or small-scale parameter scanning, and has not yet formed a unified and efficient framework suitable for large-scale inversion.

[0006] In the joint inversion of multi-source observations, a large number of studies use GNSS and InSAR data to constrain fault slip, but usually still build on simplified geometric models: either assuming a single planar fault with fixed dip angle, or using a few rectangular sub-faults to represent fault segmentation. Even if some work tries to introduce non-planar faults or curved surface approximations, their geometric construction process is often independent of the slip inversion process, lacking a unified curved surface-discrete-Green's function integrated design, making it difficult to achieve systematic optimization between geometric construction, element division, regularization and data weight, resulting in uncertainty in the interpretation of the inversion results in geometrically complex areas.

[0007] From the above domestic and foreign researches, it can be seen that a three-dimensional non-planar slip inversion model capable of simultaneously considering complex curved fault geometry, triangular dislocation discretization and multi-source observation data joint constraint is still blank, and it is urgent to innovate the corresponding numerical simulation method and develop a complex fault simulation model. The key difficulties in this field mainly lie in the following two aspects:

[0008] (1) The continuous construction and high-quality discretization method of three-dimensional complex curved fault geometry is not mature: when constructing a continuous three-dimensional fault surface from surface rupture, earthquake distribution, geological structure and other multi-source constraint information, the existing method can only obtain an approximate multi-segment plane or a local approximate surface, and it is difficult to meet the requirements of geometric smoothness and structural rationality in the entire fault range. How to discretize the triangular dislocation unit for strong bending and multi-segment fault without seams and overlaps, and to realize adaptive encryption in the area with large curvature change, is an important technical bottleneck in the current three-dimensional fault geometry modeling.

[0009] (2) The efficient Green function construction and stable inversion framework based on curved surface triangular dislocation has not been formed: traditional TDE forward usually relies on numerical integration to calculate the dislocation Green function, which has high computational cost and is complex to implement, and it is difficult to be directly used for large-scale InSAR / GNSS joint inversion. How to make full use of the efficiency of Okada analytical kernel while maintaining the flexibility of TDE geometry, construct a large-scale Green function matrix through reasonable local mapping or mixed dislocation kernel technology, and on this basis, introduce geometric adaptive regularization and slip direction constraint to realize physically consistent and computationally efficient curved fault slip inversion, is a core technical problem that needs to be broken through in the current surface deformation numerical simulation and inversion field. SUMMARY

[0010] The purpose of the present application is to provide a surface deformation simulation and three-dimensional fault inversion method, system, medium and product with high precision and stability.

[0011] The technical scheme provided by the present application is as follows:

[0012] In the first aspect, the present application provides a surface deformation simulation and three-dimensional fault inversion method, comprising:

[0013] InSAR data acquisition: acquiring InSAR data of satellite ascending and descending orbits, and obtaining surface rupture traces based on the InSAR data;

[0014] Three-dimensional non-planar fault geometry modeling: constructing a preliminary fault geometry based on the surface rupture traces, constraining the dip angle change of the fault by using relocated aftershock data, and performing local plane fitting and surface fusion processing on the preliminary fault geometry to form a continuous and smooth three-dimensional non-planar fault surface;

[0015] Constructing the hybrid dislocation model: discretizing the three-dimensional non-planar fault surface into triangular dislocation units and adaptively encrypting, further subdividing each triangular dislocation unit into multiple rectangular subunits; for each rectangular subunit, using Okada analytical kernel to calculate the displacement contribution (i.e. the displacement response generated at the observation point) of the unit slip to each observation point; based on the displacement contribution calculation result, constructing the Green function matrix to represent the linear relationship between the observation displacement and the fault slip vector;

[0016] Fault slip distribution inversion: integrating InSAR data and GNSS data, combining the Green function matrix, and through the regularization constraint and boundary condition control, obtaining the fault non-uniform slip distribution; based on the fault non-uniform slip distribution result, combining the hybrid dislocation model, completing the surface deformation simulation.

[0017] In a possible implementation manner, in the three-dimensional non-planar fault geometry modeling, the local plane fitting is achieved by projecting the three-dimensional seismic source position to a local two-dimensional strike-depth profile, and using the least square method; the strike of the local fault plane and the dip angle The calculation formula is , , The strike of the local fault plane is represented by The dip angle of the local fault plane is represented by , The linear fitting coefficients in the local fault plane fitting process are respectively

[0018] In a possible implementation manner, the subdividing each triangular dislocation unit into multiple rectangular subunits comprises: for each triangular dislocation unit, mapping it to an approximate plane in a local range; further subdividing each approximate plane into multiple rectangular subunits.

[0019] In a possible implementation manner, based on the displacement contribution calculation result, constructing the Green function matrix comprises:

[0020] For each triangular dislocation unit, the displacement contribution of the unit slip of all the rectangular subunits after splitting to each observation point is superimposed (vector superposition, considering direction and size) according to the observation point, to obtain the total displacement contribution of the triangular dislocation unit to each observation point, thereby integrating the Green function matrix Each element of the matrix corresponds to the total displacement contribution of a unit slip (strike-slip / dip-slip / crack) of a triangular dislocation unit to an observation point;

[0021] The GPU parallel computing is used to accelerate the construction process of the Green function matrix.

[0022] In a possible implementation, in the fault slip distribution inversion, an InSAR interferogram is obtained based on InSAR data; the InSAR interferogram is subjected to a quadtree downsampling processing, and InSAR line-of-sight displacement observations and GPS three-component displacement observations are uniformly included in an inversion framework; the observation displacement and a fault slip vector satisfy a linear relationship:

[0023] ;

[0024] wherein, is an observation displacement vector, which collects all InSAR line-of-sight displacement observations and GNSS three-component displacement observations; is a Green function matrix calculated by a hybrid dislocation model, is a fault slip vector, is an error vector, which represents a comprehensive effect of observation errors and model errors.

[0025] In a possible implementation, in the fault slip distribution inversion, a Tikhonov regularization is used to construct a weighted least squares objective function:

[0026] ;

[0027] wherein, is a fault slip vector, is a data covariance matrix weighted according to observation uncertainties, is a discrete Laplacian operator acting on the slip components, is a regularization coefficient for controlling a balance between model smoothness and data fitting degree;

[0028] A non-negative least squares (NNLS) is used to constrain the slip direction to be consistent with a focal mechanism and a regional tectonic background, and a zero slip boundary condition is applied to a fault boundary and a weakly geometric or data-constrained area.

[0029] In a possible implementation, based on the fault non-uniform slip distribution result, a surface deformation simulation is completed in combination with the hybrid dislocation model, including:

[0030] a fault slip vector obtained by inversion is distributed to each triangular dislocation element to obtain a corresponding slip amount of each triangular dislocation element;

[0031] ​The slip amount allocated to each triangular dislocation unit is uniformly allocated to each rectangular sub-unit according to the area proportion of the rectangular sub-unit (the larger the area of the sub-unit, the more slip amount is allocated), ensuring the conservation of total slip amount (the sum of the slip amounts of n sub-units = the total slip amount of the corresponding triangular dislocation unit);

[0032] For each rectangular sub-unit, the displacement contribution of the actual slip amount allocated to it to each observation point is calculated using the Okada analytical kernel;

[0033] For each triangular dislocation unit, the displacement contribution of the actual slip amount of all rectangular sub-units after splitting to each observation point is superimposed (vector superposition, considering direction and size), to obtain the actual displacement contribution of the triangular dislocation unit to each observation point.

[0034] For each observation point, the actual displacement contribution of all triangular dislocation units to it is superimposed to obtain the total displacement of the observation point; the total displacement of all observation points constitutes the final surface deformation simulation result.

[0035] In a second aspect, the present application provides a surface deformation simulation and three-dimensional fault inversion system, comprising a memory and a processor;

[0036] The memory is configured to store a computer program.

[0037] The processor is configured to call the computer program to execute the method as described above.

[0038] In a third aspect, the present application provides a computer readable storage medium, wherein the computer readable storage medium stores a computer program, and the computer program is configured to enable an electronic device to implement the method as described above when the computer program is run on the electronic device.

[0039] In a fourth aspect, the present application provides a computer program product comprising a computer program, wherein the computer program is configured to enable an electronic device to implement the method as described above when the computer program is run on the electronic device.

[0040] The specific implementation manners of the second to fourth aspects of the present application can refer to the implementation manners of the first aspect, which will not be described here.

[0041] Advantages:

[0042] This application proposes a method, system, medium, and product for surface deformation simulation and 3D fault inversion. By integrating the geometric flexibility of triangular dislocation units with the high computational efficiency of the Okada analytical kernel within a unified framework, it achieves accurate description of complex 3D fault geometry and high-resolution dislocation field solutions. This method can adaptively approximate fault surfaces with arbitrary curvature, bending, and segmentation characteristics while maintaining the stability of analytical solutions, effectively overcoming problems such as geometric distortion, stress discontinuity, and systematic biases caused by traditional rectangular dislocation models under curved surface conditions. By combining fault geometry smoothing, triangular dislocation unit discretization, rectangular kernel function mapping, and GPU parallel acceleration, this application constructs a 3D non-planar fault simulation method with high accuracy, high stability, and high efficiency. It can accurately reproduce the spatial geometry of complex faults and their control effect on surface deformation, providing a reliable tool for seismic tectonic analysis, fault geometry modeling, seismic hazard assessment, and crustal dynamics research. It has significant scientific value and application prospects for improving the seismic simulation and inversion capabilities of complex tectonic regions. Attached Figure Description

[0043] Figure 1 : Flowchart of an embodiment of this application.

[0044] Figure 2 : Fault segment dip angle fitting diagram in the embodiments of this application; wherein Figure 2 (a) to 2 (d) respectively show the spatial distribution characteristics of aftershocks (yellow dots) relative to the fault rupture trace (red line) under different scenarios. The horizontal axis in the figure represents the horizontal distance along the fault strike, and the vertical axis represents the depth perpendicular to the fault strike. The unit is km.

[0045] Figure 3 The embodiments of this application show two rectangular approximation schemes; wherein 3(a) shows the rectangular sub-unit shape of the triangular dislocation unit (△ABE region) with a non-uniform / extreme aspect ratio of 1:1, and 3(b) shows the rectangular sub-unit shape of the triangular dislocation unit (△ABE region) divided with a medium aspect ratio.

[0046] Figure 4 Comparison of errors between two rectangular approximation schemes (aspect ratios of 1:1 and 1:0000, respectively) in the embodiments of this application;

[0047] Figure 5 : Comparison of GPU acceleration and original model computation speed in the embodiments of this application;

[0048] Figure 6 : A comparison diagram of the simulation results of this application and the Okada analytical reference solution in the embodiments of this application; wherein Figure 6(a) Comparison of the distribution of the direction displacement (u) of the present application model and the Okada model at different depths (-4.0 km, -2.0 km, 0.0 km depth) Figure 6 (b) Comparison of the distribution of the direction displacement (u) of the present application model and the Okada model at different depths Figure 6 (c) Comparison of the distribution of the z direction displacement (w) of the present application model and the Okada model at different depths Figure 6 (d) Comparison of the distribution of the derivative of displacement (u) with respect to (x) of the present application model and the Okada model at different depths Figure 6 (e) Comparison of the distribution of the derivative of displacement (u) with respect to (y) of the present application model and the Okada model at different depths Figure 6 (f) Comparison of the distribution of the derivative of displacement (u) with respect to (z) of the present application model and the Okada model at different depths

[0049] Figure 7 : The three-dimensional geometric graph and X-Z side view of the fault surface using the precise surface discretization and multi-segment rectangular approximation in the embodiments of the present application; wherein Figure 7 (a) is a three-dimensional comparison graph of the fault surface using the precise surface discretization and multi-segment rectangular approximation; Figure 7 (b) is an X-Z cross section comparison graph of the fault surface using the precise surface discretization and multi-segment rectangular approximation; wherein the coordinate axes x, y are horizontal directions (unit: km), and z is a depth direction (unit: km, negative value represents underground); red / blue / green respectively corresponds to different rectangular / rectangular projection cross sections.

[0050] Figure 8 : Comparison graphs of the simulation results (curved surface-the present application model) of the present application, the Okada analytical reference solution (rectangular-Okada model) and the tdcalc analytical reference solution (curved surface-tdcalc model) in the embodiments of the present application; wherein Figure 8 (a) Comparison of the distribution of the direction displacement (u) of the present application model and the Okada model at different depths (-4.0 km, -2.0 km, 0.0 km depth) Figure 8 ​​​​​​​​​​​​​​​​(b) Comparison of the distributions of the horizontal displacement (u) at different depths for the curved-dislocation model, the rectangular-Okada model, and the curved-tdcalc model. Figure 8 (c) Comparison of the distributions of the z-direction displacement (w) at different depths for the curved-dislocation model, the rectangular-Okada model, and the curved-tdcalc model. Figure 8 (d) Comparison of the distributions of the derivative of the displacement (u) with respect to the depth at different depths for the curved-dislocation model, the rectangular-Okada model, and the curved-tdcalc model. Figure 8 (e) Comparison of the distributions of the derivative of the displacement (w) with respect to the depth at different depths for the curved-dislocation model, the rectangular-Okada model, and the curved-tdcalc model. Figure 8 (f) Comparison of the distributions of the derivative of the displacement (u) with respect to the depth at different depths for the curved-dislocation model, the rectangular-Okada model, and the curved-tdcalc model.

[0051] Figure 9 Figure 1: Comparison of the curved fault and its grid discretization schemes in the embodiments of the present application; wherein Figure 9 (a) is the rectangular grid discretization scheme; Figure 9 (b) is the triangular dislocation element grid discretization scheme.

[0052] Figure 10 Figure 2: Comparison of the simulation results of the curved-dislocation model and the analytical reference solution of the tdcalc model in the embodiments of the present application; wherein Figure 10 (a) Comparison of the distributions of the horizontal displacement (u) at different depths (-4.0 km, -2.0 km, 0.0 km depth) for the curved-dislocation model and the tdcalc model, Figure 10 (b) Comparison of the distributions of the horizontal displacement (u) at different depths for the curved-dislocation model and the tdcalc model; Figure 10 (c) Comparison of the distributions of the z-direction displacement (w) at different depths for the curved-dislocation model and the tdcalc model, Figure 10 (d) Comparison of the distributions of the derivative of the displacement (u) with respect to the depth at different depths for the curved-dislocation model and the tdcalc model. ​​​​​​​​​​​​​​​​​​​) distribution; Figure 10 (e) Comparison of displacement of the model of the present application and the tdcalc model at different depths of the present application derivative of the present application ) distribution; Figure 10 (f) Comparison of displacement of the model of the present application and the tdcalc model at different depths of the present application derivative of the present application ) distribution.

[0053] Figure 11 : Comparison of analytical reference solution simulation velocity of the present application and the tdcalc in the embodiment of the present application;

[0054] Figure 12 : Earthquake surface rupture trajectory and aftershock distribution map of a certain region in the embodiment of the present application;

[0055] Figure 13 : Fault geometry three-dimensional diagram of a certain region in the embodiment of the present application;

[0056] Figure 14 : Multi-perspective diagram of fault slip distribution of a certain region in the embodiment of the present application; wherein Figure 14 (a) is a back view of the fault slip distribution; Figure 14 (b) is a front view of the fault slip distribution. DETAILED DESCRIPTION

[0057] In order to make the person in the art better understand the present application scheme, the technical scheme of the present application will be further described in detail below in combination with the embodiments of the present application and the drawings.

[0058] The actual fault plane of the continental active tectonic belt has the characteristics of bending, segmentation, multi-scale roughness, etc. The variation of the dip angle in the strike and depth directions will form strong three-dimensional geometric heterogeneity. However, the existing seismic deformation simulation and slip inversion mostly use a single dip angle or a multi-segment spliced planar rectangular dislocation model to approximate the real fault. This geometric simplification method easily introduces systematic geometric distortion in the strong curvature or segmented transition zone of the fault: it will not only cause the overlap or gap of adjacent rectangular sub-faults, resulting in non-physical high values or sudden changes of the local strain and stress field, but also cannot restore the continuous curvature of the fault surface, and needs to introduce "artificial complexity" to compensate for the geometric error, thereby causing the estimation deviation of the seismic moment, the distortion of the position and size of the strong slip area, and the misjudgment of the Coulomb stress change field. At the same time, with the development of high-resolution geodetic measurement technologies such as InSAR and GNSS, the planar simplified model is difficult to fully absorb the spatial information of the observation data, and is also easy to produce inversion residuals and non-physical slip concentration artifacts at the bending and segmented connection of the fault geometry, affecting the quantitative understanding of the segmentation of the earthquake rupture, the strain concentration zone and the stress transfer of adjacent faults. In addition, the traditional triangular dislocation method has low calculation efficiency and cannot support large-scale inversion constrained by multi-source observation, which makes it difficult to reliably predict the displacement field, slip distribution and stress change of complex fault systems.

[0059] In view of the above-mentioned key problems of the traditional method, such as insufficient fault geometry description, low-efficiency dislocation calculation, unstable inversion results and easy-to-produce non-physical artifacts, the present application constructs an integrated multi-source geodetic data, high-precision curved surface construction, discrete and slip inversion method system that adapts to the characteristics of complex three-dimensional faults. Under the constraint of multi-source geology and geophysics, the system generates a continuous and smooth three-dimensional non-planar fault surface, uses triangular dislocation units to discretize the surface with high quality, and develops a hybrid dislocation simulation technology combining triangular units and Okada analytical kernel to realize efficient solving of the surface deformation Green function on the complex curved surface. On this basis, the system introduces geometric self-adaptive regularization and slip direction constraint to establish a stable and physically consistent curved surface fault slip inversion framework, and accurately restores the spatial distribution characteristics and strain release mode of the real rupture.

[0060] The application establishes a complete technical process from InSAR and GNSS data processing, fault geometry reconstruction, curved surface grid discretization, mixed dislocation core construction to geometric self-adaptive constraint inversion, realizes high-resolution slip solving of complex three-dimensional fault structure and accurate simulation of ground surface displacement. The method is not only suitable for coseismic deformation inversion, fault geometry reconstruction, but also can be extended to stress change analysis, tectonic deformation research, earthquake risk assessment, ground surface deformation monitoring and other three-dimensional dislocation geophysical inversion problems, can provide high-precision tools for continental tectonic seismic dynamics research and earthquake monitoring engineering, significantly improve the reliability of coseismic slip model and stress change analysis, support earthquake risk assessment and early warning and tectonic deformation research, and has wide prospect in complex fault system related research and engineering application in continental tectonic region, and has important scientific significance and engineering application value.

[0061] The specific embodiments according to the present application will be described below with reference to the accompanying drawings.

[0062] Embodiment one:

[0063] As shown in the figure, the embodiment of the application provides a ground surface deformation simulation and three-dimensional fault inversion method, comprising: Figure 1

[0064] S1, InSAR data acquisition: acquiring InSAR data of satellite ascending and descending orbits, and obtaining ground surface rupture traces based on the InSAR data.

[0065] In order to depict the coseismic deformation field caused by an earthquake, the application processes InSAR data of Sentinel-1 satellite ascending and descending orbits, which covers the ground surface rupture zone. Both images are obtained by using the interferometric wide (IW) mode, and reasonable azimuth resolution and range resolution are set. By selecting image pairs with reasonable time and vertical baseline, the interference coherence in the study area is ensured.

[0066] S2, three-dimensional non-planar fault geometry modeling: constructing a preliminary fault geometry based on the ground surface rupture traces, using relocated aftershock data to constrain the dip angle change of the fault, and performing local plane fitting and curved surface fusion processing on the preliminary fault geometry to form a continuous and smooth three-dimensional non-planar fault surface.

[0067] ​Natural faults often exhibit significant curvature, kinks, segmentation, and other geometric irregularities in three-dimensional space, and the geometry of a fault has a significant control on the complexity of the rupture behavior. If a complex fault plane is directly discretized by rectangular elements, it often leads to unrealistic numerical artifacts such as concentrated slip, discontinuous slip, over-smoothed stress-strain field, and biased stress-strain field estimation due to the overlap or gap between adjacent elements. Therefore, it is necessary to construct a reliable three-dimensional fault surface. When the complex fault geometry constrained by aftershock distribution is introduced as prior information into inversion, the fault slip model constrained by InSAR observations is usually more robust, and the model uncertainty caused by geometric simplification can be effectively reduced. Compared with geodetic data (GNSS, InSAR, etc.) and focal mechanism solutions, relocated high-precision aftershock data can better characterize the deep fault structure with higher spatial resolution, and have stronger constraints on three-dimensional features such as dip angle variation and geometric bending.

[0068] To solve the above problems and accurately capture the geometric characteristics of the three-dimensional curved fault surface in the source area, the application first constructs the preliminary geometry of the fault based on the surface rupture trace. Then, the aftershock is used to constrain the dip angle variation of the fault, so as to construct a three-dimensional non-planar fault geometry closer to the true underground morphology. Assuming that aftershocks mainly occur along the main fault plane, the aftershock dense area can be regarded as the spatial projection of the underground fault. The fault length is determined by the surface rupture and the aftershock segmentation characteristics. According to the spatial aggregation of aftershocks along the strike, the aftershocks are divided into several sub-areas; in each sub-area, the three-dimensional source position is projected onto the local two-dimensional strike-depth profile, and the least squares method is used to fit the local fault plane to extract the dip angle and its variation along the strike Figure 2 . The strike and dip angle of the fault segment are:

[0069] ;

[0070] wherein, represents the strike of the local fault plane; represents the dip angle of the local fault plane; , are the linear fitting coefficients in the local fault plane fitting process.

[0071] Each fitted plane represents a local fault geometry.

[0072] After completing the local plane fitting, redundant nodes are deleted and a multi-segment cubic spline interpolation method is used to fuse all local fault planes to make the entire fault surface have geometric continuity and smoothness in space (first and second order continuity), and avoid the problems of boundary discontinuity, normal jump and non-physical changes in curvature.

[0073] S3, constructing a mixed dislocation model: discretizing the three-dimensional non-planar fault surface into triangular dislocation elements and adaptively encrypting, further subdividing each triangular dislocation element into multiple rectangular sub-elements; for each rectangular sub-element, using Okada analytical kernel to calculate the displacement contribution (i.e. the displacement response generated at the observation point) of the unit slip to each observation point; based on the displacement contribution calculation result, constructing a Green function matrix to represent the linear relationship between the observed displacement and the fault slip vector.

[0074] Since Steketee (1958) first introduced dislocation theory into geodetic deformation research, it has been developed for nearly seven decades and widely applied to calculate the surface displacement, strain, gravity change and surface topography change caused by earthquakes and other problems. Early research includes the analytical expression of point dislocation in Poisson solid derived by Steketee, and the near-field solution of vertical rectangular strike-slip fault given by Chinnery (1961). These works laid the foundation for the theoretical relationship between fault slip and surface deformation in a homogeneous elastic half-space. With the deepening understanding of the source mechanism and the complexity of fault geometry, dislocation theory has gradually expanded to adapt to more general three-dimensional fault configurations. Among them, Okada (1985, 1992) proposed a complete and widely used family of analytical formulas based on the summary of previous research, which can be used to calculate the displacement, strain and tilt produced by rectangular faults with arbitrary inclination in a homogeneous elastic half-space. For any observation point , the displacement produced by a rectangular dislocation surface can be written as:

[0075] ;

[0076] where, represents the displacement amount generated by fault slip at any observation point in a homogeneous elastic half-space; , , represent the strike-slip, dip-slip and opening amount respectively; , and are all Okada Green functions, which explicitly depend on the geometry of the fault, the medium properties and the relative position of the observation point and the dislocation surface.

[0077] After the establishment of continuous three-dimensional fault geometry, the curved surface is divided into high-quality triangular dislocation elements (TDEs) to meet the flexible expression ability of non-planar faults in geometry. However, traditional triangular dislocation elements rely on numerical integration for displacement and stress calculation, which is extremely costly and difficult to meet the large-scale and efficient operation required by InSAR / GNSS joint inversion. To solve this problem, this application proposes a hybrid calculation strategy, that is, through "curved surface triangular dislocation element → local plane mapping → multiple rectangular sub-element approximation → Okada analytical kernel" to realize the efficient coupling between TDE and analytical solution. This method maintains the high-fidelity description ability of triangular dislocation elements to complex curved surface fault geometry, and fully inherits the high precision and high computational efficiency of the Okada model in half-space dislocation calculation. After completing the triangular dislocation element division, each triangular dislocation element is further subdivided into multiple rectangular sub-elements to efficiently couple with the Okada (1985, 1992) rectangular dislocation analytical kernel.

[0078] Wherein, the aspect ratio of the rectangular sub-element can be verified by multiple sets of comparative numerical experiments, and the optimal aspect ratio interval is selected as the core standard for the minimum numerical oscillation and the highest stability of the mixed dislocation model calculation.

[0079] This application carries out numerical experiment verification on different rectangular approximation schemes (aspect ratios are 1:1 and 1:10000) Figure 3 By comparing the two rectangular approximation schemes, it is found that the configuration with moderate aspect ratio can effectively reduce numerical oscillation and improve calculation stability Figure 4 Therefore, each triangular dislocation element is further subdivided into multiple rectangular sub-elements with moderate aspect ratio (such as 1:10000-1:10). Further, in order to more accurately capture the local changes of slip, the triangular grid is adaptively refined by local curvature, so that the high curvature area or the area with high slip gradient can obtain higher geometric resolution, thereby improving the accuracy of subsequent slip inversion. However, grid refinement will significantly increase the calculation amount. In order to further improve the calculation efficiency, this application combines the GPU parallel computing of the PyTorch platform in the construction of the mixed dislocation model, which accelerates the construction process of large-scale Green function matrix from CPU environment to about 30 times or more Figure 5 ), making it possible to include tens of thousands to hundreds of thousands of data points in InSAR and GNSS joint inversion. This high-performance computing framework makes this application applicable to regional-scale complex three-dimensional fault slip inversion, significantly improving the efficiency bottleneck of traditional methods.

[0080] S4, fault slip distribution inversion: integrate InSAR data and GNSS data, combine with Green's function matrix, through regularization constraint and boundary condition control, the non-uniform slip distribution of the fault is obtained by inversion; based on the non-uniform slip distribution result of the fault, combined with the mixed dislocation model, the surface deformation simulation is completed.

[0081] The application constructs a three-dimensional slip inversion framework combined with InSAR and GNSS. The InSAR interferogram adopts quadtree downsampling to reduce the data amount and computing cost as much as possible while retaining the long-wavelength deformation characteristics, and unifies the InSAR line-of-sight displacement observation and GPS three-component displacement observation into the same inversion framework. The Green's function matrix is generated by the mixed dislocation model , so that the observation displacement and the fault slip vector satisfy the linear relationship:

[0082] ;

[0083] wherein, is the observation displacement vector, which collects all InSAR line-of-sight displacement observations and GNSS three-component displacement observations; is the Green's function matrix calculated by the mixed dislocation model, is the fault slip vector (including the strike-slip and dip-slip components on each triangular dislocation element, and the cracking component can be flexibly added according to the actual fault rupture type), is the error vector, which represents the comprehensive effect of observation error and model error.

[0084] The application adopts Tikhonov regularization to construct a weighted least squares objective function, and the formula is:

[0085] ;

[0086] wherein, is the data covariance matrix weighted according to the observation uncertainty, is the discrete Laplacian operator acting on the slip component, is the regularization coefficient for controlling the balance between model smoothness and data fitting degree.

[0087] wherein, The regularization term is introduced. The present application introduces a geometric adaptive smoothing method, and the smoothing strength of the regularization term is not globally uniform, but is dynamically adjusted according to the local geometric characteristics (local curvature and segmented structure) of the fault surface - in the area with small curvature and gentle geometry, the smoothing strength is appropriately enhanced to avoid the occurrence of high-frequency fluctuations in the slip distribution which is not physically meaningful; in the area with large curvature such as the bending area of the fault and the transition zone of the segment, the smoothing strength is automatically weakened to retain the local details (such as slip gradient and peak concentration characteristics) of the slip distribution, and to avoid excessive smoothing or non-physical slip rupture in the bending area of the fault.

[0088] To ensure the geophysical consistency of the inversion results, the present application further uses the non-negative least squares method (NNLS) to constrain the slip direction, so that it is consistent with the focal mechanism and regional tectonic background. At the same time, zero slip boundary conditions are applied at the fault boundary and in areas with weak geometric or data constraints to suppress the false concentration of slip caused by boundary effects. Among them, the area with weak geometric or data constraints refers to the area within the range of the fault surface, where the constraint ability for fault geometry construction and slip distribution inversion is significantly reduced due to the complexity of its own geometric characteristics or insufficient support of InSAR / GNSS multi-source observation data (insufficient coverage, accuracy and effectiveness of observation data).

[0089] The optimal regularization parameter is determined by the L-curve method to achieve the best trade-off between data fitting ability and model smoothing degree. Overall, the inversion system composed of geometric adaptive regularization-NNLS effectively improves the stability, resolution and physical reasonableness of slip inversion.

[0090] The above surface deformation simulation and three-dimensional fault inversion method proposed by the present application reconstructs a continuous and smooth three-dimensional curved surface fault through multi-source constraints, and realizes high-precision dislocation calculation for any curvature by using a hybrid strategy combining triangular dislocation elements and Okada analytical kernel. Combined with adaptive grid and GPU parallel acceleration, the present application can stably invert the non-uniform slip distribution of complex faults, avoid the geometric distortion and stress discontinuity problems of traditional rectangular models, and has the advantages of high efficiency, high precision and wide applicability, which can provide reliable technical tools for seismic deformation simulation, three-dimensional fault structure reconstruction and seismic risk assessment.

[0091] The accuracy, stability and geometric adaptability of the three-dimensional non-planar fault geometric modeling and hybrid dislocation modeling method proposed by the present application in multiple typical tectonic scenarios are verified below.

[0092] a) Numerical consistency verification under the condition of planar fault

[0093] To verify the accuracy of the proposed method under simple geometric conditions, a idealized planar fault model with uniform slip was first constructed, and the traditional Okada analytical model and the proposed three-dimensional hybrid dislocation method were used for simulation respectively. Since the planar fault is the strict application condition of Okada model, it can be used to verify the numerical consistency of the proposed method under the condition of no curvature. Figure 6 The experimental results show that the displacement components and strain gradient curves of the two models on different depth profiles are almost completely consistent (

[0094] b) Model verification under simple curved fault conditions

[0095] The proposed hybrid three-dimensional fault model was further compared with tdcalc (a boundary element implementation of triangular dislocation elements; Maerten et al., 2014) as a benchmark. tdcalc is a boundary element calculation program that can be used to solve the displacement, strain and stress caused by planar triangular dislocation elements in a full or half-space elastic medium. By combining multiple triangular dislocation elements, this model can simulate fault surfaces with complex curvature.

[0096] To evaluate the geometric expression ability and numerical accuracy of the proposed method under curved fault conditions, a three-dimensional curved fault surface with slight curvature variation was further constructed. The surface was discretized by two different modeling strategies: (1) accurate surface discretization (TDE model): directly use high-quality triangular dislocation elements (tdcalc thought framework) to strictly geometrically discretize the curved surface as a high-precision reference solution; (2) multi-planar rectangular approximation (Okada simplified model): simplify the curved fault into three adjacent planar rectangular fault pieces to simulate their displacement field with typical Okada kernel. Both models were subjected to a uniform slip of 0.5 m, and the sampling grid was kept consistent to allow for comparative testing ( Figure 7 ).

[0097] The simulation results show that the proposed hybrid three-dimensional dislocation method is highly consistent with the benchmark TDE model (tdcalc) in terms of displacement distribution, gradient variation and spatial form, with minimal differences. The traditional segmented rectangular planar approximation model shows significant deformation deviation in the fault bending area, especially in the curvature variation area, resulting in amplitude error or waveform distortion ( Figure 8). These errors are caused by the fact that the planar model cannot accurately represent the continuous curvature of the fault in the strike and depth directions, thus destroying the geometric continuity of the fault surface and the simulation results. Comprehensive analysis shows that the traditional rectangular planar model cannot accurately approximate the curved fault surface, which will introduce a non-negligible systematic error in the deformation simulation; in contrast, the method of the present application can accurately capture the continuity and curvature variation of the curved fault surface in three-dimensional space, and has a significant advantage in the simulation of complex geometric faults.

[0098] c) Non-uniform slip simulation under real curved fault conditions and high-precision solution verification

[0099] To further verify the applicability of the method of the present application under real complex fault geometry and non-uniform slip conditions, a three-dimensional non-planar fault model with significant spatial variation is constructed. The fault surface is reconstructed by spline smoothing technology, and high-quality discretization is performed using adaptive triangular elements to maintain the geometric continuity and local curvature characteristics of the curved surface in the strike and depth directions. On this curved surface, each triangular dislocation element is assigned a different slip magnitude (Sx, Sy, Sz) Figure 9 ) to simulate the non-uniform slip distribution in real seismic rupture.

[0100] The comparison results show that the method of the present application is highly consistent with the high-precision tdcalc model in simulating the surface displacement and displacement gradient field generated by the non-uniform slip fault ( Figure 10 ). It is worth noting that under the same geometric discretization conditions, the computational efficiency of the hybrid method of the present application is significantly better than that of the tdcalc model: as the grid resolution increases, the numerical integration cost of tdcalc increases sharply, while the dislocation solution based on the analytical kernel of the present application can maintain high computational efficiency ( Figure 11 ), showing a high ability to apply to large-scale inversion and regional seismic simulation.

[0101] The above verification results fully prove that the present application not only has high precision in describing the displacement field in geometrically complex regions, but also is much more efficient than the traditional boundary element implementation of the TDE model in terms of computational efficiency, and is particularly suitable for large-scale joint inversion of InSAR / GNSS and high-resolution seismic rupture model construction.

[0102] Simulation experiment:

[0103] The method is used to simulate the slip process of a seismic fault in a certain area in history. Surface investigation and InSAR deformation show that the earthquake produced a surface rupture zone of about 160-170 km. Using the double-difference algorithm of Wang et al. (2022), the aftershocks within nine days after the main shock are relocated, and after spatial filtering, 1,288 high-quality events are obtained, which clearly depict a NW-SE trending fault zone, showing significant segmentation characteristics: dense in the west, sparse in the middle, and gaps in the east, with most events distributed in the 5-15 km depth range Figure 12 ). According to the surface rupture zone obtained from field investigation, the trace of the fault is determined, and the variation of the fault dip is constrained by the three-dimensional distribution of aftershocks, and a three-dimensional non-planar fault geometry closer to the true underground shape is constructed. Further, the fault slip distribution is inverted using the fault slip inversion method.

[0104] The preferred three-dimensional fault model obtained in this application extends from the surface (z=0 km) to about 25 km, which is highly consistent with the depth distribution of aftershock relocation results and the surface rupture shape obtained from InSAR and field investigation Figure 13 ). In the plane, the fault trace shows a slow arc change: the west segment is nearly NW-SE, the middle segment gradually changes to E-W, and the east segment further rotates to ENE. Along the strike direction, the dip also shows systematic changes: the west segment is steeply north-dipping, the middle segment is nearly vertical, and the east segment gradually transitions to south-dipping, forming a bowl-shaped three-dimensional curved surface as a whole, rather than a simple single planar fault.

[0105] The inverted slip model corresponds to a seismic moment magnitude of Mw7.44, which is basically consistent with the teleseismic source inversion result. The coseismic slip is highly non-uniform in space, mainly concentrated in the 0-20 km depth range on the hanging wall of the fault, and rapidly decays below ~22-25 km Figure 14 ). Along the strike, four main slip peak areas can be identified: the largest main rupture zone is located about 5-10 km deep east of the epicenter, with a maximum slip of 5.16 m; the remaining three slip peaks are developed in the west and middle segments, as well as near the shallow southeast end of the fan-shaped bifurcation. Along the strike, the east and west sides of the epicenter show obvious asymmetry: the east segment has deeper slip and more concentrated peaks, while the west segment has more dispersed slip and slightly lower amplitude.

[0106] Embodiment Two

[0107] The embodiment provides a surface deformation simulation and three-dimensional fault inversion system, comprising a storage and a processor;

[0108] The storage is used to store a computer program;

[0109] The processor is configured to invoke the computer program to execute the method in Embodiment One.

[0110] Embodiment Three

[0111] The embodiment provides a computer readable storage medium, and the computer readable storage medium stores a computer program. The computer program is used to enable an electronic device to implement the method in Embodiment One when the computer program is run on the electronic device.

[0112] Embodiment Four

[0113] The embodiment provides a computer program product, and the computer program product includes a computer program. The computer program is used to enable an electronic device to implement the method in Embodiment One when the computer program is run on the electronic device.

[0114] The specific implementation manners of the system, the electronic device, the computer readable storage medium and the computer program product provided in the embodiment of the present application can refer to the specific embodiments of the above method, and will not be described here.

[0115] Obviously, those skilled in the art should understand that each unit or each step of the above-mentioned present application can be realized by using a general computing device, and they can be concentrated on a single computing device, or distributed on a network composed of multiple computing devices, and optionally, they can be realized by using program codes executable by a computing device, so that they can be stored in a storage device and executed by a computing device, or they can be respectively manufactured into each integrated circuit module, or multiple modules or steps among them can be manufactured into a single integrated circuit module to realize. Therefore, the present application is not limited to any specific combination of hardware and software.

[0116] The above only describes the preferred embodiments of the present application and is not used to limit the present application. For those skilled in the art, the present application can have various modifications and changes. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A method of surface deformation simulation and 3D fault inversion, characterized in that, include: InSAR data acquisition: Acquire InSAR data of satellite ascending and descending orbits, and obtain surface rupture traces based on the InSAR data; Three-dimensional non-planar fault geometry modeling: Based on the surface rupture trace, the preliminary geometry of the fault is constructed, the dip angle of the fault is constrained by the relocation aftershock data, and the preliminary geometry of the fault is subjected to local plane fitting and surface fusion processing to form a continuous and smooth three-dimensional non-planar fault surface. Constructing a hybrid dislocation model: The three-dimensional non-planar fault surface is discretized into triangular dislocation units and adaptively refined. Each triangular dislocation unit is further subdivided into multiple rectangular sub-units. For each rectangular sub-unit, the Okada analytical kernel is used to calculate the displacement contribution of unit slip to each observation point. Based on the displacement contribution calculation results, a Green's function matrix is ​​constructed to represent the linear relationship between the observed displacement and the fault slip vector. Fault slip distribution inversion: Integrating InSAR and GNSS data, and combining the Green's function matrix, the non-uniform fault slip distribution is inverted through regularization constraints and boundary condition control. Based on the results of the non-uniform fault slip distribution, and combined with the aforementioned hybrid dislocation model, surface deformation simulation is completed, including: The fault slip vectors obtained by inversion are distributed to each triangular dislocation unit to obtain a corresponding slip amount distributed to each triangular dislocation unit. The fault slip vectors obtained by inversion are distributed to each triangular dislocation unit to obtain a corresponding slip amount distributed to each triangular dislocation unit. The slip amount allocated to each triangular dislocation unit is evenly distributed to each rectangular sub-unit according to the area ratio of the rectangular sub-units to ensure the conservation of total slip amount; For each rectangular sub-element, the Okada analytical kernel is used to calculate the contribution of its assigned actual slip to the displacement of each observation point; For each triangular dislocation element, the actual slip of all its sub-rectangular elements after splitting is added together to obtain the actual displacement contribution of the triangular dislocation element to each observation point. For each observation point, the actual displacement contribution of all triangular dislocation elements is superimposed to obtain the total displacement of that observation point; the total displacement of all observation points constitutes the final surface deformation simulation result.

2. The method according to claim 1, characterized in that, In the three-dimensional non-planar fault geometry modeling, the local plane fitting is achieved by projecting the three-dimensional seismic source position to a local two-dimensional strike-depth profile using the least square method; the strike of the local fault plane and the dip angle The calculation formula is , , The strike of the local fault plane is represented by The dip angle of the local fault plane is represented by , The linear fitting coefficients in the local fault plane fitting process are respectively The surface fusion adopts a multi-segment cubic spline interpolation method to fuse all the local fault planes to obtain a three-dimensional non-planar fault surface, and the three-dimensional non-planar fault surface satisfies the first and second order continuity to ensure continuous smoothness.

3. The method according to claim 1, characterized in that, The step of subdividing each triangular dislocation unit into multiple rectangular subunits includes: mapping each triangular dislocation unit to an approximate plane within a local range; and further subdividing each approximate plane into multiple rectangular subunits.

4. The method according to claim 1, characterized in that, Based on the displacement contribution calculation results, the Green's function matrix is ​​constructed, including: For each triangular dislocation element, the displacement contribution of the unit slip of all its sub-elements to each observation point is summed at each observation point to obtain the total displacement contribution of the triangular dislocation element to each observation point. This summation yields the Green's function matrix. Each element of the matrix corresponds to a unit slip of a triangular dislocation unit and its contribution to the total displacement of an observation point. GPU parallel computing is used to accelerate the construction process of the Green's function matrix.

5. The method according to claim 1, characterized in that, In the fault slip distribution inversion, InSAR interferograms are obtained based on InSAR data; quadtree downsampling is applied to the InSAR interferograms to unify the InSAR line-of-sight displacement observations and GPS three-component displacement observations into the inversion framework; the observed displacements... With fault slip vector Satisfies a linear relationship: ； in, To observe the displacement vector, all InSAR line-of-sight displacement observations and GNSS three-component displacement observations are collected; The Green's function matrix is ​​calculated using the mixed dislocation model. The fault slip vector, This is the error vector, representing the combined effect of observation error and model error.

6. The method according to claim 1, characterized in that, In the fault slip distribution inversion, the regularization constraint uses Tikhonov regularization to construct a weighted least squares objective function: ； in, The fault slip vector, This is the data covariance matrix weighted by observation uncertainty. For the discrete Laplace operator acting on the slip component, Regularization coefficients are used to control the tradeoff between model smoothness and data fit. By using the non-negative least squares (NNLS) method to constrain the slip direction to be consistent with the source mechanism and regional tectonic background, zero slip boundary conditions are applied at fault boundaries and in areas with weak geometric or data constraints.

7. A surface deformation simulation and three-dimensional fault inversion system, characterized in that, include: Memory and processor; The memory is used to store computer programs; The processor is configured to invoke the computer program to perform the method as described in any one of claims 1 to 6.

8. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed on an electronic device, causes the electronic device to perform the method as described in any one of claims 1 to 6.

9. A computer program product, comprising a computer program, characterized in that, When the computer program is run on an electronic device, it causes the electronic device to perform the method as described in any one of claims 1 to 6.

Images