Earthquake deformation and stress distribution calculation method and system based on Green function

Abstract

The invention belongs to the field of earthquake and stress distribution calculation, and particularly relates to an earthquake deformation and stress distribution calculation method and system based on a Green function. The method comprises the following steps: acquiring a triangular mesh and a to-be-measured point of a fault; constructing an earth fixed coordinate system, a triangular dislocation coordinate system and an angular dislocation coordinate system based on the triangular mesh of the fault; determining a transformation matrix to obtain a coordinate system; calculating interior angles of three vertexes of the triangular mesh in the coordinate system and corresponding supplementary angles, selecting equivalent configuration based on the supplementary angles, calculating barycentric coordinates of orthogonal projection of the to-be-measured point on the triangular dislocation plane, and selecting configuration based on a partition inequality of the barycentric coordinates; according to configuration, a Burgers function is obtained through calculation by using a solid angle formula, incomplete displacement contributions of all angular dislocations are superposed, and a complete displacement field containing slip continuity correction is obtained; a stress field is calculated based on a complete displacement field or through a Green function of strain. And the seismic deformation and stress distribution calculation performance is obviously improved.

Description

Technical Field

[0001] This invention belongs to the field of earthquake and stress distribution calculation, specifically relating to a method and system for calculating earthquake deformation and stress distribution based on Green's function. Background Technology

[0002] The calculation of earthquake deformation and stress distribution is an important research area in seismology, geophysics, and engineering geology, used to assess earthquake risk, predict ground motion, analyze precursors of volcanic activity, and guide seismic design and disaster early warning. Traditional methods are mostly based on elastic medium theory, using Green's function to simulate the displacement, stress, and deformation fields caused by earthquake sources. Green's function, as a core tool for elastic wave propagation and deformation calculation, is widely used in the integral solution of point-source dislocations to simulate the response of complex geological structures. For example, based on the Green's function method proposed by Steketee (1958), a point-source dislocation solution for half-space models has been developed and has become the standard framework for seismic motion simulation. Since the 1950s, this method has been widely used in near-field and far-field earthquake simulations after theoretical refinement and numerical optimization, and plays a crucial role in practical earthquake risk assessment.

[0003] In the existing technology, common calculation methods based on Green's function include: (1) Empirical Green's function method, which is used to simulate strong ground motion. It uses small earthquake events as empirical Green's function to correct the amplitude spectrum of the target event. However, it has strict requirements for the selection of small earthquakes (such as similar magnitude and close location). In complex media, it is easily affected by path effects and differences in source mechanisms, resulting in a decrease in accuracy; (2) Numerical Green's function method, which combines the spectral element method or the finite difference method to build a Green's function library in a three-dimensional layered Earth model. It is used for source rupture process inversion or centroid moment tensor (CMT) calculation. For example, its application in the Australian continent model can handle deep earthquakes. Earthquakes with a magnitude of up to 50 km are involved, but they involve a large number of numerical integrations and series summations, resulting in low computational efficiency. In particular, the algorithm is complex and resource-intensive in the self-weight compressible sphere model. (3) The stress drop calculation method based on Lg waves corrects the path effect through a broadband attenuation model and fits the observed and theoretical spectra to estimate the stress distribution. This method shows advantages in areas such as the Qinghai-Tibet Plateau, but it is sensitive to data quality and cannot process large-scale events in real time. (4) Passive source monitoring technology uses noise cross-correlation to reconstruct the empirical Green function to realize the time-varying monitoring of structures, such as volcanic activity identification. However, it relies on high signal-to-noise ratio data and is greatly affected by environmental interference.

[0004] Furthermore, within the framework of dislocation theory, Green's functions are commonly used to simulate finite dislocation sources, such as rectangular dislocation (RD) and triangular dislocation (TD) models. The rectangular dislocation solution proposed by Okada (1992) calculates displacement and stress fields in half-space and has been widely applied to earthquake and volcanic deformation modeling. However, limited by rectangular geometry, it cannot flexibly simulate curved surfaces or complex faults. In contrast, the triangular dislocation model approximates arbitrary surfaces through triangular meshes, providing greater geometric flexibility. Early TD solutions were based on Yoffe's (1960) angular dislocation superposition, solving for displacement and stress in full space. However, they suffer from artifact singularities and numerical instability along the extended edges of the TD, leading to incorrect results when the calculation points are close to these lines. Subsequent improvements include the half-space TD solutions of Jeyakumaran et al. (1992) and Meade (2007), which approximate the elimination of artifact singularities through numerical integration and Taylor series expansion. However, the output is scale-dependent and computationally inefficient. Nikkhoo and Walter (2015) proposed a fully analytical TD solution that avoids pseudo-singularities and numerical instability in both full and half-space. Through two equivalent angular dislocation configurations and a barycentric coordinate partitioning strategy, they ensure computation is independent of dislocation size, location, and slip vector. This method improves Burgers function computation and is suitable for multi-scale applications, but further optimization is still needed to handle complex media responses in real-time emergency scenarios.

[0005] Despite the progress made by these existing technologies in calculating seismic deformation and stress distribution, the following shortcomings remain: First, accuracy is limited by model simplification, such as neglecting the Earth's curvature and layered structure in semi-infinite spaces, which increases errors in global or complex geological environments; second, computational efficiency is low, as numerical methods require the construction of large Green's function libraries or integration, making it difficult to respond to large earthquake emergencies in real time; third, they have poor adaptability to multi-source, non-homogeneous media, with prominent boundary effects and scale problems, especially since pseudo-singularities in the TD model, although partially resolved, still affect reliability due to numerical dependence; fourth, there is a lack of efficient hybrid frameworks, making it impossible to seamlessly integrate analytical and numerical solutions to balance accuracy and time.

[0006] Therefore, existing methods for calculating seismic deformation and stress distribution based on Green's function still need improvement in terms of accuracy, efficiency, and applicability to meet the practical needs of seismic risk assessment, volcano monitoring, and engineering applications. This invention aims to solve these problems by providing a novel calculation method that improves the accuracy and efficiency of the calculations. Summary of the Invention The purpose of this invention is to address the shortcomings of existing methods for calculating seismic deformation and stress distribution, such as limited accuracy due to model simplification, low computational efficiency, poor adaptability to multi-source inhomogeneous media, and lack of an efficient hybrid framework. To achieve the above objectives, the present invention adopts the following technical solution: In a first aspect, the present invention provides a method for calculating seismic deformation and stress distribution based on Green's function, comprising: Obtain the triangular mesh of the fault and the points to be measured; Based on the triangular mesh of the fault, we construct the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system; and determine the transformation matrix between the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system to obtain the coordinate system. Calculate the interior angles of the three vertices of the triangular mesh in the coordinate system and their corresponding supplementary angles. Based on the supplementary angles, select an equivalent configuration. Calculate the centroid coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane. Select a configuration based on the partition inequality of the centroid coordinates. Based on the selected configuration, the Burgers function is calculated using the solid angle formula, and the incomplete displacement contribution of each angular dislocation is superimposed to obtain a complete displacement field including slip continuity correction. The stress field can be calculated based on the complete displacement field or through the Green's function of strain.

[0007] Preferably, the fault-based triangular mesh constructs a fixed Earth coordinate system, a triangular dislocation coordinate system, and an angular dislocation coordinate system, and determines the transformation matrix between these three systems to obtain the coordinate system; specifically: Construct a fixed coordinate system for the Earth with East, North, and Up as the X, Y, and Z axes; With the second vertex of the triangular mesh of the fault as the origin, the x-axis is perpendicular to the triangular dislocation plane, the y-axis is along the strike direction, and the z-axis is along the dip direction to construct a triangular dislocation coordinate system; Using the coordinates of the unit normal vector, strike vector, and dip vector of the fault's triangular mesh in the Earth's fixed coordinate system as columns of the transformation matrix, and combining them with the translation matrix, the transformation matrix from the triangular dislocation coordinate system to the Earth's fixed coordinate system is determined. Based on the triangular dislocation coordinate system, an angular dislocation coordinate system is constructed, and the transformation matrix from the angular dislocation coordinate system to the triangular dislocation coordinate system is determined.

[0008] Preferably, the calculation involves determining the interior angles and corresponding supplementary angles of the three vertices of the triangular mesh within the coordinate system, selecting an equivalent configuration based on the supplementary angles, calculating the barycenter coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane, and selecting a configuration based on the partition inequality of the barycenter coordinates, including: Calculate the interior angles of the three vertices of the triangular mesh: , , The supplementary angle is ; Using the supplementary angle as the angular dislocation input, and employing two equivalent angular dislocation configurations, the barycenter coordinates of the orthogonal projection of the test point onto the triangular dislocation plane are calculated. ,in ; Configuration selection based on partitioning inequalities using centroid coordinates.

[0009] Preferably, the centroid coordinates are calculated using the following formula:

[0010] In the formula, Represents the x-coordinate of the receiving point; Indicates the x-coordinate of the first triangular unit; Indicates the x-coordinate of the second triangle unit; Indicates the x-coordinate of the third triangle unit; Represents the ordinate of the first triangular element; Represents the ordinate of the receiving point; Indicates the ordinate of the second triangular unit; Indicates the ordinate of the third triangle element; Indicates the coordinates of the first centroid; This indicates the coordinates of the second centroid.

[0011] Preferably, the incomplete displacement contributions of each angular dislocation are superimposed to obtain a complete displacement field including slip continuity correction, specifically as follows:

[0012] in, It is a fault plane. Corresponding to the x, y, and z directions. It is the Burgers function; It is an incomplete displacement.

[0013] Preferably, the calculation of the stress field based on the complete displacement field or through the Green's function of strain specifically includes: The strain tensor is calculated from the complete displacement field based on triangular dislocation sources:

[0014] Calculating the stress field using Hooke's law:

[0015] in, It is Lamé's constant; For strain tensor; subscript and Represents the coordinate direction, typically taking the value... or ; It is the displacement vector in Components of direction; representing coordinates The partial derivatives; Let be the stress tensor.

[0016] Preferably, at any point in space, the deformation and stress distribution caused by the fault are a linear superposition of the calculation results of all discrete elements.

[0017] Preferably, when superimposing the incomplete displacement contribution of each angular dislocation, the angle between the angular dislocation and the z-axis is determined. When the angle is less than 1e-2 radians, the angular dislocation component and the stress component are directly set to 0.

[0018] In a second aspect, the present invention provides a calculation system for seismic deformation and stress distribution based on Green's function, comprising: The acquisition unit is used to acquire the triangular mesh of the fault and the points to be measured. Transformation units are constructed to build the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system based on the triangular mesh of the fault; and the transformation matrix between the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system is determined to obtain the coordinate system. The first calculation unit is used to calculate the interior angles of the three vertices of the triangular mesh in the coordinate system and the corresponding supplementary angles. Based on the supplementary angles, an equivalent configuration is selected, and the centroid coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane are calculated. The configuration is selected based on the partition inequality of the centroid coordinates. The second calculation unit, based on the selected configuration, calculates the Burgers function using the solid angle formula, and superimposes the incomplete displacement contributions of each angular dislocation to obtain a complete displacement field including slip continuity correction. The third calculation unit calculates the stress field based on the complete displacement field or through the Green's function of strain.

[0019] In a third aspect, the present invention provides an electronic device including a processor and a memory, the processor being configured to execute a computer program stored in the memory to implement the steps of the calculation method for seismic deformation and stress distribution based on Green's function as described in any one of the preceding claims.

[0020] In a fourth aspect, the present invention provides a computer-readable storage medium storing at least one instruction that, when executed by a processor, implements the steps of the method for calculating seismic deformation and stress distribution based on Green's function as described in any one of the preceding claims.

[0021] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention significantly improves the performance of seismic deformation and stress distribution calculations by improving the Green's function-based triangular dislocation (TD) model and combining it with a fully analytical solution method. The specific effects are as follows: High-precision simulation, free from pseudo-singularity interference. Traditional methods often suffer from computational errors due to pseudo-singularities along the extension lines of dislocation edges and numerical instability, especially when the calculation point is close to these regions. This invention employs two equivalent angular dislocation configurations and partitions the calculation points by the barycentric coordinates on the TD plane through orthogonal projection (based on matrix determinant calculation). ,in By selecting a configuration far from the pseudo-singularity line, displacement and stress field calculations are ensured to be independent of dislocation scale and location. Compared to existing methods, the accuracy is significantly improved, and it is applicable to complex fault geometries (such as curved faults or magma chambers).

[0022] Highly efficient computation, supporting real-time applications. Traditional numerical methods, requiring the construction of large Green's function libraries and extensive integration, suffer from low computational efficiency, making them unsuitable for emergency response to major earthquakes. This invention utilizes a fully analytical solution (requiring no numerical approximation or series expansion), combined with an improved Burgers function (…). ,in (Calculated via vector triple product), reducing computational complexity and significantly shortening computation time. The MATLAB code implementation facilitates rapid deployment and supports real-time monitoring of earthquake sequences or precursors of volcanic activity, with a response speed several times faster than existing methods.

[0023] Enhanced applicability and adaptability to complex geological environments: Existing methods are mostly based on simplified models (such as semi-infinite spaces), which have limited performance in non-homogeneous media or multi-source scenarios. This invention simulates arbitrary curved surfaces using triangular meshes, extending to half-space (using the mirror principle to superimpose the contributions of the principal dislocation (TD) and the mirror TD), and integrates analytical solutions with numerical methods (such as the boundary element method), effectively handling complex terrain and multi-fault interactions. It is applicable to seismic activity at depths up to 50 kilometers or complex regions such as the Qinghai-Tibet Plateau, and its applicability is superior to traditional rectangular dislocation models.

[0024] Numerical Stability and Verifiability: The singularity problem in traditional Volterra dislocation theory leads to computational instability. This invention eliminates numerical instability near the 1 / r singularity through geometric methods and an improved Burgers function. Displacement Calculation Formula The results are consistent with existing analytical solutions (such as Okada, 1992), ensuring their reliability and applicability to engineering design and academic research.

[0025] Multi-scale applications and engineering value: This invention supports deformation simulation from local fault slip to global scales, through coordinate transformation matrices and slip vector components. Its flexible input adapts to the diverse needs of modern measurements such as InSAR and GPS. In engineering, it can guide seismic structural design and disaster early warning, resulting in significant economic and social benefits.

[0026] In summary, this invention surpasses existing technologies in terms of accuracy, efficiency, applicability, and stability, filling the gap in real-time, efficient deformation calculation under complex geological conditions and possessing broad scientific and engineering application prospects. Attached Figure Description

[0027] The accompanying drawings, which form part of this application, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings: Figure 1 This is a schematic diagram of the method flow according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the partition configuration according to an embodiment of the present invention; Figure 3 This is a diagram showing the calculated surface deformation results caused by the rupture of the Eastern Anatolian Fault in an embodiment of the present invention. Figure 4 This is a comparison diagram of the surface deformation observed by GPS measurements (Qi et al., 2011) and simulated data obtained in an embodiment of the present invention; wherein, (a) is the horizontal diagram of GPS; and (b) is the vertical diagram of GPS. Figure 5 The stress variation diagrams on the Wenchuan-Maoxian fault caused by the Mw7.9 Wenchuan earthquake, calculated using Green's function in this embodiment of the invention, are as follows: (a) is the normal stress variation diagram; (b) is the shear stress variation diagram; (c) is... The Coulomb stress variation diagram; (d) is The Coulomb stress variation diagram; Figure 6 This is a system structure block diagram according to an embodiment of the present invention; Figure 7 This is a structural block diagram of an electronic device according to an embodiment of the present invention. Detailed Implementation

[0028] The present invention will now be described in detail with reference to the accompanying drawings and embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other.

[0029] The following detailed description is exemplary and intended to provide further detailed explanation of the invention. Unless otherwise specified, all technical terms used in this invention have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains. The terminology used in this invention is for the purpose of describing particular embodiments only and is not intended to limit the scope of exemplary embodiments according to the invention.

[0030] See Figure 1 This application discloses a method for calculating seismic deformation and stress distribution based on Green's function, including: S1: Obtain the triangular mesh of the fault and the points to be measured; S2: Based on the triangular mesh of the fault, construct the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system; and determine the transformation matrix between the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system to obtain the coordinate system; S3: Calculate the interior angles of the three vertices of the triangular mesh in the coordinate system and the corresponding supplementary angles. Based on the supplementary angles, select the equivalent configuration, calculate the centroid coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane, and select the configuration based on the partition inequality of the centroid coordinates. S4: Based on the selected configuration, the Burgers function is calculated using the solid angle formula, and the incomplete displacement contribution of each angular dislocation is superimposed to obtain a complete displacement field including slip continuity correction. S5: Calculate the stress field based on the complete displacement field or through the Green's function of strain.

[0031] This method utilizes the fully analytical solution of the triangular dislocation (TD) model, avoiding pseudo-singularities and numerical instability, and achieves accurate calculation of displacement and stress fields. Specifically, this invention aims to: Improved computational accuracy: By using analytical solutions of TD in full space and half space, pseudo-singularities along the extension line of dislocation edges are eliminated, ensuring that the calculation results are independent of dislocation scale, location and slip vector, which is suitable for deformation simulation of complex geological structures.

[0032] Improve computational efficiency: Employ angular dislocation superposition and barycentric coordinate partitioning strategies to reduce the need for numerical integration and series summation, supporting real-time processing of large-scale earthquake events and emergency responses.

[0033] Enhanced applicability: Integrating analytical and numerical methods to form a hybrid framework, it can flexibly simulate the response of arbitrary curved faults, multi-source dislocations and non-uniform media, and is suitable for earthquake risk assessment, volcano monitoring and seismic design of engineering projects.

[0034] By achieving the above objectives, the present invention can better meet the practical needs of the fields of seismology and geophysics and promote the progress of related technologies.

[0035] In some embodiments, the fault-based triangular mesh constructs a fixed Earth coordinate system, a triangular dislocation coordinate system, and an angular dislocation coordinate system, and determines the transformation matrix between these three systems to obtain the coordinate system; specifically: Construct a fixed coordinate system for the Earth with East, North, and Up as the X, Y, and Z axes; With the second vertex of the triangular mesh of the fault as the origin, the x-axis is perpendicular to the triangular dislocation plane, the y-axis is along the strike direction, and the z-axis is along the dip direction to construct a triangular dislocation coordinate system; Using the coordinates of the unit normal vector, strike vector, and dip vector of the fault's triangular mesh in the Earth's fixed coordinate system as columns of the transformation matrix, and combining them with the translation matrix, the transformation matrix from the triangular dislocation coordinate system to the Earth's fixed coordinate system is determined. Based on the triangular dislocation coordinate system, an angular dislocation coordinate system is constructed, and the transformation matrix from the angular dislocation coordinate system to the triangular dislocation coordinate system is determined.

[0036] In some embodiments, the calculation of the interior angles and corresponding supplementary angles of the three vertices of the triangular mesh within the coordinate system, selecting an equivalent configuration based on the supplementary angles, calculating the barycentric coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane, and selecting a configuration based on the partition inequality of the barycentric coordinates, includes: Calculate the interior angles of the three vertices of the triangular mesh: , , The supplementary angle is ;in, , and These represent the interior angles of the triangle at vertices 1, 2, and 3, respectively. Indicates from vertex Pointing to the vertex The unit direction vector, ; Using the supplementary angle as the angular dislocation input, and employing two equivalent angular dislocation configurations, the barycenter coordinates of the orthogonal projection of the test point onto the triangular dislocation plane are calculated. ,in ; Configuration selection based on partitioning inequalities using centroid coordinates.

[0037] In some embodiments, the centroid coordinates are calculated using the following formula:

[0038] In the formula, Represents the x-coordinate of the receiving point; Indicates the x-coordinate of the first triangular unit; Indicates the x-coordinate of the second triangle unit; Indicates the x-coordinate of the third triangle unit; Represents the ordinate of the first triangular element; Represents the ordinate of the receiving point; Indicates the ordinate of the second triangular unit; Indicates the ordinate of the third triangle element; Indicates the coordinates of the first centroid; Indicates the coordinates of the second barycenter. This indicates the coordinates of the third centroid.

[0039] In some embodiments, the superposition of the incomplete displacement contributions of each angular dislocation to obtain a complete displacement field including slip continuity correction specifically involves:

[0040] in, It is the Burgers function; These are the components of the complete displacement field in the x, y, and z directions; , These are the Burgers vector components in the x, y, and z directions; It is the first The contribution of incomplete displacement in the x-direction of each angular dislocation; It is the first The contribution of incomplete displacement in the y-direction of each angular dislocation; It is the first The contribution of incomplete displacement in the z-direction of an angular dislocation.

[0041] In some embodiments, the calculation of the stress field based on the complete displacement field or through the Green's function of strain specifically includes: The strain tensor is calculated from the complete displacement field based on triangular dislocation sources:

[0042] Calculating the stress field using Hooke's law:

[0043] in, For strain tensor; subscript and Represents the coordinate direction; It is the displacement vector in Component of direction; Represents coordinates The partial derivatives; It is the displacement vector in Component of direction; Represents coordinates The partial derivatives; For stress tensor; It is Lamé's constant; For indicator functions, This is the volumetric strain term.

[0044] In some embodiments, at any point in space, the deformation and stress distribution caused by the fault are a linear superposition of the calculation results of all discrete elements.

[0045] In some embodiments, when superimposing the incomplete displacement contributions of each angular dislocation, the angle between the angular dislocation and the z-axis is determined. When the angle is less than 1e-2 radians, the angular dislocation component and stress component are directly set to 0. This invention discloses a method for calculating seismic deformation and stress distribution based on Green's function, including: obtaining the triangular mesh of the fault and the points to be measured; constructing a fixed Earth coordinate system, a triangular dislocation coordinate system, and an angular dislocation coordinate system based on the triangular mesh, and determining the transformation matrix between them; then calculating the interior angles and supplementary angles of the vertices of the triangular mesh, selecting an equivalent configuration, and further selecting a configuration based on the centroid coordinate partitioning inequality; then, according to the selected configuration, calculating the Burgers function using the solid angle formula, superimposing the incomplete displacement contributions of each angular dislocation, and performing special processing when the angle between the angular dislocation vector and the Z-axis is close to 0° to obtain the displacement field of a single element; based on the displacement field and stress field results, strain and stress fields can be calculated using constitutive equations. The fault deformation, strain, and stress at any point in space should be a linear superposition of the calculation results of all elements. This method improves computational accuracy and efficiency, and is suitable for simulating displacement and stress fields caused by earthquakes and for disaster assessment.

[0046] In some embodiments, this invention discloses a method for calculating seismic deformation and stress distribution based on Green's function. This method utilizes the fully analytical solution of the triangular dislocation (TD) model in both full and half-space, employing angular dislocation superposition, a barycentric coordinate partitioning strategy, and an improved Burgers function to avoid pseudo-singularities and numerical instability, achieving accurate simulation of displacement, stress, and deformation fields. This method is applicable to seismic risk assessment of complex geological structures, volcano monitoring, aftershock research, tsunami modeling, and seismic design of engineering projects, supporting multi-scale applications and real-time calculations. Specifically, this invention models the seismic source or deformation source (such as faults or magma chambers) as a triangular dislocation surface, using Volterra dislocation theory and the Green's function method to solve for the elastic medium response, avoiding the limitations of traditional rectangular dislocation models in simulating curved surface geometry.

[0047] The core of this method lies in modeling the earthquake source as a triangular dislocation, solving for the response in the elastic medium through Green's function integration, applying Yoffe's (1960) principle of triangular dislocation superposition, and eliminating pseudo-singularities along the extension of the TD edge through two equivalent configurations. The specific steps are as follows: Establish coordinate systems and transformation matrices: Define the Earth Fixed Coordinate System (EFCS), the Triangular Dislocation Coordinate System (TDCS), and the Angular Dislocation Coordinate System (ADCS). The EFCS uses East, North, and Up as the X, Y, and Z axes, respectively; the TDCS uses the second vertex P2 of TD as its origin, with the x-axis perpendicular to the TD plane (normal), the y-axis along the azimuth direction, and the z-axis along the dip direction. Use transformation matrices to perform coordinate and vector transformations. For example, use the A_TE matrix for a direct transformation from TDCS to EFCS, where its columns are the unit normal vectors of TD. , directional vector and tilt vector Coordinates in EFCS, combined with translation matrix Implement position vector transformation. Among them, , , Similarly, the A_AT matrix is ​​used for the transformation from TDCS to ADCS, for example, for the first ADCS:

[0048] This matrix is ​​defined and normalized based on the TD vertex coordinates, ensuring accurate conversion of calculation points and slip vectors between different systems. The inverse transformation uses matrix transpose and corresponding translation to achieve coordinate transformation of displacement, stress, and strain tensors.

[0049] Calculate the internal angles of the TD and select the configuration: Calculate the internal angles of the TD , , ,in Equal to the unit vector of the TD side, and using the supplementary angle ( The angle dislocation is used as the input. Two equivalent angle dislocation configurations (Config-1 and Config-2) are adopted. The barycenter coordinates of the orthogonal projection of the calculation point onto the TD plane are calculated (normalized coordinates are calculated using matrix determinant). ,in The configuration is selected based on the partitioning inequality to avoid calculation points being close to pseudo-singularity lines. For example, the barycenter coordinates are calculated as follows:

[0050] The TD plane is divided into six zones, using the centerline ( (etc.) serves as a boundary to ensure that the calculation point is far from the semi-infinite dislocation line when selecting the configuration, thereby eliminating numerical instability. See details for specific configuration. Figure 2 The triangular dislocation plane has seven partitions and two configurations. The blue and red dashed lines represent the singular lines corresponding to the two configurations, respectively. Each singular line is located within only one configuration partition. To avoid artificial singularities and numerical instability, this application uses different configurations for different receiving points. For points within the TD, both configurations can be used.

[0051] Calculating the Burgers function: An improvement on Yoffe's (1960) Burgers function, calculated using the solid angle formula. ( ,in Let P be the solid angle subtended by the dislocation surface at the calculation point P, calculated using the following formula (ensuring single-valuedness and numerical stability):

[0052] Used to correct discontinuities in displacement. The solid angle corresponding to the triangle. One is the Burgers function, and the rest are vectors, to improve and eliminate numerical instability in the Yoffe solution.

[0053] Calculate displacement and stress field: superimpose the incomplete displacements and stresses of three angular dislocations (incomplete displacements) Through the glide vector component The contribution calculation involves , , (An algebraic expression with equal parameters), combined with the slip vector and the Burgers function, yields the complete displacement:

[0054] Calculation of stress tensor using Hooke's law etc., among which pass (Partial derivatives and algebraic terms are calculated), supporting full-space models. Displacement and stress fields are calculated in TDCS and then converted to EFCS.

[0055] To improve the numerical stability of angular dislocation calculations, we need to introduce an angle threshold. When the included angle of the angular dislocation model is close to 0° (i.e., almost parallel to the z-axis), it may lead to singularities or numerical instability in the Burgers function or displacement / stress calculations (such as accuracy loss when the input to the acos function is close to 1, or the denominator of the solid angle formula approaching zero at small angles). Therefore, when the calculated included angle beta is less than the threshold 1e-2 radians (approximately 0.0057°, a very small threshold that ensures triggering only when approximately parallel), the components of the angular dislocation (e.g., the Burgers vector contribution) and the associated stress components are directly set to 0. This avoids unnecessary computational overhead and improves the overall simulation accuracy.

[0056] Extending to half-space and hybrid frameworks: For half-space, the mirror principle is used to extend the full-space solution, superimposing the contributions of the main TD and the mirror TD, and calculating the harmonic function contribution (avoiding surface pseudo-singularities through equivalent configuration). Integrating analytical solutions with numerical methods (such as the boundary element method) forms a hybrid framework, suitable for non-homogeneous media, multi-source simulation, and realistic terrain simulation, improving adaptability to complex geological environments.

[0057] Actual fault geometry is often complex (bending, branching, segmentation, etc.), so it is usually discretized into a large number of small rectangular or triangular dislocation elements. The size of each element is much smaller than the scale of the study area, and the slip (displacement vector) is assumed to be uniform within the element. Therefore, the effects of multiple dislocation sources at the same point can be directly linearly superimposed.

[0058]

[0059] Where N is the total number of discrete elements, and superscript(k) represents the field generated at point x when the k-th element acts alone.

[0060] This method is implemented using Python code, and the input parameters include the TD vertex coordinates, the slip vector, and the Poisson's ratio. This invention calculates the coordinates of the points and outputs the displacement and stress field. Compared to existing methods, this invention improves accuracy (independent of scale, position, and slip vector), efficiency (reducing numerical integration and series summation, supporting real-time emergency response), and applicability (flexibly simulating arbitrary curved surface faults and multi-source dislocations). Verification shows consistency with solutions such as the Okada rectangular dislocation solution, making it suitable for multi-scale geophysical problems.

[0061] Example 1 In seismology, coseismic slip refers to the slip distribution along a fault during an earthquake, typically obtained through inversion from observational data such as InSAR and GPS. Calculating surface deformation (usually displacement) and stress using Green's function is a standard method based on elastic dislocation theory. Green's function describes the response of a unit slip source to an observation point; the total response can be obtained through integration or superposition. This invention provides a fully analytical solution for the triangular dislocation (TD) model to simulate displacement and stress caused by finite sources, avoiding the problem of pseudo-singularities. The detailed steps and mathematical expressions are as follows.

[0062] 1. Prepare data and model Obtaining coseismic slip distribution: The slip vector distribution on the fault is obtained by inversion from observation data (such as InSAR or GPS displacement fields), which is usually represented as slip vector. (Including directional sliding, tilt sliding, and opening components), among which It's the source point location. The sliding vector is used... These represent the directions corresponding to the normal, strike, and dip angle, respectively.

[0063] Choose an Earth model: Use a uniform half-space or full-space model; Building the Green function library: Green functions Represents the source point The unit j-direction sliding relative to the observation point The displacement in the i-direction response.

[0064] 2. Calculate surface deformation (displacement field) Mathematical expression: Earth's surface displacement Calculated by integration using the Green's function:

[0065] in It is a fault plane. Corresponding to the x, y, and z directions.

[0066] For the TD model, the displacement is corrected by superimposing the incomplete displacements of the three angular dislocations and adding the Burgers function:

[0067] in It is the Burgers function (based on solid angles) calculate), Equal to incomplete displacement.

[0068] step: The fault is discretized into triangular meshes (TDs), each TD having its own vertex coordinates and sliding vector.

[0069] Establish a coordinate system: using EFCS (Earth-fixed coordinates), TDCS (triangular dislocation coordinates), and ADCS (angular dislocation coordinates), through transformation matrices (such as... and Transform coordinates and vectors.

[0070] Calculate internal angles and supplementary angle Choose the equivalent configuration (Config-1 or Config-2) based on the centroid coordinates. Partitioning helps avoid false singularities.

[0071] For each observation point, the displacement contribution is calculated and summed. Figure 3 and Figure 4 The results of surface deformation calculations for the 2023 Türkiye earthquake and the 2008 Wenchuan earthquake are provided. Figure 3 The calculated surface deformation results for the rupture of the East Anatolian Fault show the Green's function calculations for the East Anatolian Fault and the Nurda Fault. The map shows the surface displacement distribution in the Fault area. The amount of surface displacement is represented by a color gradient, ranging from blue (displacement below -3 meters) to red (displacement above +3 meters), passing through green (displacement close to 0 meters). Arrows and lines indicate the strike and direction of slip of the fault, with accompanying color bands showing the gradual change in displacement. The legend provides the logarithm of the slip rate (log10(s)) and a reference range for surface displacement (-3 meters to +3 meters). Geographic coordinates from 36.5° to 38.5° latitude and 36° to 38° longitude cover the area. Figure 4 This paper compares observed GPS measurements (Qi et al., 2011) with simulated surface deformation. Figures show the horizontal (a) and vertical (b) ground displacements from GPS. Surface deformation data derived from the Global Positioning System (GPS) provides valuable constraints for the Wenchuan earthquake model. The Green's function calculations are compared with the GPS deformation data (Qi et al., 2011). In most areas, the static surface deformation in the model is largely consistent with the observed results, demonstrating the effectiveness of the Green's function in simulating surface deformation. However, some significant horizontal errors exist, particularly in the central and northern parts of the Beichuan Fault.

[0072] 3. Calculate the stress field Stress Tensor It can be obtained through the Green's function of strain or calculated from displacement. First, calculate the strain. Then apply Hooke's Law: in It is Lamé's constant.

[0073] step: Strain is calculated from the displacement field (via partial derivatives and algebraic terms calculation wait).

[0074] Stress is obtained by applying Hooke's law. The TD model can directly calculate stress, avoiding numerical integration.

[0075] For Coulomb stress changes (used to assess aftershock risk): ,in It is a change in normal stress. It is a change in shear stress. It is the effective coefficient of friction. Figure 5 The calculation results of Coulomb stress changes on the Wenmao Fault caused by the Wenchuan earthquake are provided, with the surface trajectory of the Beichuan Fault marked by red lines. WC - Wenchuan Town; MX - Maoxian Town. Based on Green's function, the stress changes on other faults caused by coseismic slip during the Wenchuan earthquake are calculated, characterizing the degree of seismic hazard. Its mathematical expression is (Harris, 1998):

[0076] in and The difference between the normal stress and shear stress and the initial state is denoted as . This indicates an increase in compressive stress. This is the static friction coefficient. Generally speaking, positive... A negative ΔCFS characterizes fault instability, corresponding to increased seismic activity; while a negative ΔCFS characterizes the stress shadow zone, corresponding to decreased seismic activity. The calculated static Coulomb stress changes on the Wenmao Fault after the Wenchuan earthquake are shown below. Figure 5 As shown, the main shock may have significantly reduced the normal and shear stresses on the Wenmao Fault, but the stress distribution on this fault plane is uneven. The largest stress changes are mainly located in the southern part of the Wenmao Fault, near the Beichuan Fault. Because the reduction in yield stress is less than the reduction in shear stress in this area, ΔCFS is significantly reduced. In the middle and upper parts of the Wenmao Fault, because the reduction in yield stress is greater than that in shear stress, ΔCFS is significantly enhanced in this area, gradually receding northward, and mainly distributed near the surface in the northern segment of the Wenmao Fault.

[0077] Example 2 like Figure 6 As shown, based on the same inventive concept as the above embodiments, the present invention also provides a calculation system for seismic deformation and stress distribution based on Green's function, comprising: The acquisition unit is used to acquire the triangular mesh of the fault and the points to be measured. Transformation units are constructed to build the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system based on the triangular mesh of the fault; and the transformation matrix between the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system is determined to obtain the coordinate system. The first calculation unit is used to calculate the interior angles of the three vertices of the triangular mesh in the coordinate system and the corresponding supplementary angles. Based on the supplementary angles, an equivalent configuration is selected, and the centroid coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane are calculated. The configuration is selected based on the partition inequality of the centroid coordinates. The second calculation unit, based on the selected configuration, calculates the Burgers function using the solid angle formula, and superimposes the incomplete displacement contributions of each angular dislocation to obtain a complete displacement field including slip continuity correction. The third calculation unit calculates the stress field based on the complete displacement field or through the Green's function of strain.

[0078] Example 3 like Figure 1 As shown, the present invention also provides an electronic device 100 for implementing the steps of a method for calculating seismic deformation and stress distribution based on Green's function; The electronic device 100 includes a memory 101, at least one processor 102, a computer program 103 stored in the memory 101 and executable on at least one processor 102, and at least one communication bus 104.

[0079] The memory 101 can be used to store the computer program 103. The processor 102 implements the steps of the calculation method of seismic deformation and stress distribution based on Green's function by running or executing the computer program stored in the memory 101 and calling the data stored in the memory 101.

[0080] The memory 101 may primarily include a program storage area and a data storage area. The program storage area may store the operating system, application programs required for at least one function (such as sound playback function, image playback function, etc.), etc.; the data storage area may store data created based on the use of the electronic device 100 (such as audio data), etc. In addition, the memory 101 may include non-volatile memory, such as hard disk, RAM, plug-in hard disk, smart media card (SMC), secure digital (SD) card, flash card, at least one disk storage device, flash memory device, or other non-volatile solid-state storage device.

[0081] At least one processor 102 may be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. Processor 102 may be a microprocessor or any conventional processor. Processor 102 is the control center of electronic device 100, connecting various parts of electronic device 100 via various interfaces and lines.

[0082] The memory 101 in the electronic device 100 stores multiple instructions to implement a method for calculating seismic deformation and stress distribution based on Green's function, and the processor 102 can execute multiple instructions to achieve the following: S1: Obtain the triangular mesh of the fault and the points to be measured; S2: Based on the triangular mesh of the fault, construct the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system; and determine the transformation matrix between the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system to obtain the coordinate system; S3: Calculate the interior angles of the three vertices of the triangular mesh in the coordinate system and the corresponding supplementary angles. Based on the supplementary angles, select the equivalent configuration, calculate the centroid coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane, and select the configuration based on the partition inequality of the centroid coordinates. S4: Based on the selected configuration, the Burgers function is calculated using the solid angle formula, and the incomplete displacement contribution of each angular dislocation is superimposed to obtain a complete displacement field including slip continuity correction. S5: Calculate the stress field based on the complete displacement field or through the Green's function of strain.

[0083] Example 4 If the modules / units integrated in the electronic device 100 are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, all or part of the processes in the methods of the above embodiments of the present invention can also be implemented by a computer program instructing related hardware. The computer program can be stored in a computer-readable storage medium, and when executed by a processor, it can implement the steps of the various method embodiments described above. The computer program includes computer program code, which can be in the form of source code, object code, executable files, or certain intermediate forms. The computer-readable medium can include: any entity or device capable of carrying computer program code, recording media, USB flash drives, portable hard drives, magnetic disks, optical disks, computer memory, and read-only memory (ROM).

[0084] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0085] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations 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, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0086] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0087] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0088] In the description of this specification, references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0089] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A method for calculating seismic deformation and stress distribution based on Green's function, characterized in that, include: Obtain the triangular mesh of the fault and the points to be measured; Based on the triangular mesh of faults, construct the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system; The transformation matrices between the Earth-fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system are determined to obtain the coordinate system. Calculate the interior angles of the three vertices of the triangular mesh in the coordinate system and their corresponding supplementary angles. Based on the supplementary angles, select an equivalent configuration. Calculate the centroid coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane. Select a configuration based on the partition inequality of the centroid coordinates. Based on the selected configuration, the Burgers function is calculated using the solid angle formula, and the incomplete displacement contribution of each angular dislocation is superimposed to obtain a complete displacement field including slip continuity correction. The stress field can be calculated based on the complete displacement field or through the Green's function of strain.

2. The method for calculating seismic deformation and stress distribution based on Green's function according to claim 1, characterized in that, The fault-based triangular mesh constructs a fixed Earth coordinate system, a triangular dislocation coordinate system, and an angular dislocation coordinate system, and determines the transformation matrices between these three systems to obtain the coordinate system; specifically: Construct a fixed coordinate system for the Earth with East, North, and Up as the X, Y, and Z axes; With the second vertex of the triangular mesh of the fault as the origin, the x-axis is perpendicular to the triangular dislocation plane, the y-axis is along the strike direction, and the z-axis is along the dip direction to construct a triangular dislocation coordinate system; Using the coordinates of the unit normal vector, strike vector, and dip vector of the fault's triangular mesh in the Earth's fixed coordinate system as columns of the transformation matrix, and combining them with the translation matrix, the transformation matrix from the triangular dislocation coordinate system to the Earth's fixed coordinate system is determined. Based on the triangular dislocation coordinate system, an angular dislocation coordinate system is constructed, and the transformation matrix from the angular dislocation coordinate system to the triangular dislocation coordinate system is determined.

3. The method for calculating seismic deformation and stress distribution based on Green's function according to claim 1, characterized in that, The computational coordinate system calculates the interior angles and corresponding supplementary angles of the three vertices of the triangular mesh. Based on the supplementary angles, an equivalent configuration is selected. The centroid coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane are calculated. The configuration is selected based on the partition inequality of the centroid coordinates, including: Calculate the interior angles of the three vertices of the triangular mesh: , , The supplementary angle is ; Using the supplementary angle as the angular dislocation input, and employing two equivalent angular dislocation configurations, the barycenter coordinates of the orthogonal projection of the test point onto the triangular dislocation plane are calculated. Configuration selection based on partitioning inequalities using centroid coordinates; in, , and These represent the interior angles of the triangle at vertices 1, 2, and 3, respectively. Indicates from vertex Pointing to the vertex The unit direction vector, .

4. The method for calculating seismic deformation and stress distribution based on Green's function according to claim 3, characterized in that, The centroid coordinates are calculated using the following formula: In the formula, Represents the x-coordinate of the receiving point; Indicates the x-coordinate of the first triangular unit; Indicates the x-coordinate of the second triangle unit; Indicates the x-coordinate of the third triangle unit; Represents the ordinate of the first triangular element; Represents the ordinate of the receiving point; Indicates the ordinate of the second triangular unit; Indicates the ordinate of the third triangle element; Indicates the coordinates of the first centroid; Indicates the coordinates of the second barycenter. This indicates the coordinates of the third centroid.

5. The method for calculating seismic deformation and stress distribution based on Green's function according to claim 1, characterized in that, The incomplete displacement contributions of each angular dislocation are superimposed to obtain a complete displacement field including slip continuity correction, specifically: in, It is the Burgers function; These are the components of the complete displacement field in the x, y, and z directions; , These are the Burgers vector components in the x, y, and z directions; It is the first The contribution of incomplete displacement in the x-direction of each angular dislocation; It is the first The contribution of incomplete displacement in the y-direction of each angular dislocation; It is the first The contribution of incomplete displacement in the z-direction of an angular dislocation.

6. The method for calculating seismic deformation and stress distribution based on Green's function according to claim 1, characterized in that, The calculation of the stress field based on the complete displacement field or through the Green's function of strain specifically involves: The strain tensor is calculated from the complete displacement field based on triangular dislocation sources: Calculating the stress field using Hooke's law: in, For strain tensor; subscript and Represents the coordinate direction; It is the displacement vector in Component of direction; Represents coordinates The partial derivatives; It is the displacement vector in Component of direction; Represents coordinates The partial derivatives; For stress tensor; It is Lamé's constant; For indicator functions, This is the volumetric strain term.

7. The method for calculating seismic deformation and stress distribution based on Green's function according to claim 1, characterized in that, At any point in space, the deformation and stress distribution caused by the fault are a linear superposition of the calculation results of all discrete elements.

8. The method for calculating seismic deformation and stress distribution based on Green's function according to claim 1, characterized in that, When superimposing the incomplete displacement contribution of each angular dislocation, the angle between the angular dislocation and the z-axis is determined. When the angle is less than 1e-2 radians, the angular dislocation component and stress component are directly set to 0.

9. A calculation system for seismic deformation and stress distribution based on Green's function, characterized in that, include: The acquisition unit is used to acquire the triangular mesh of the fault and the points to be measured. Transformation units are constructed to build the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system based on the triangular mesh of the fault; and the transformation matrix between the Earth fixed coordinate system, the triangular dislocation coordinate system, and the angular dislocation coordinate system is determined to obtain the coordinate system. The first calculation unit is used to calculate the interior angles of the three vertices of the triangular mesh in the coordinate system and the corresponding supplementary angles. Based on the supplementary angles, an equivalent configuration is selected, and the centroid coordinates of the orthogonal projection of the point to be measured onto the triangular dislocation plane are calculated. The configuration is selected based on the partition inequality of the centroid coordinates. The second calculation unit, based on the selected configuration, calculates the Burgers function using the solid angle formula, and superimposes the incomplete displacement contributions of each angular dislocation to obtain a complete displacement field including slip continuity correction. The third calculation unit calculates the stress field based on the complete displacement field or through the Green's function of strain.

10. An electronic device, characterized in that, It includes a processor and a memory, the processor being used to execute a computer program stored in the memory to implement the steps of the method for calculating seismic deformation and stress distribution based on Green's function as described in any one of claims 1 to 8.

Images