Fault finite element analysis method, device and storage medium based on viscoelastic model

Abstract

The application discloses a fault finite element analysis method based on a viscoelastic model, and discloses a computer device and a storage medium. The search mechanism of the viscoelastic model is introduced, the optimal three-dimensional viscoelastic finite element model of a target fault area is obtained, the response of the fault locking degree to different stress field conditions can be comprehensively analyzed, the evolution characteristics of the fault under the viscoelastic model can be understood, and in particular, the important influence of the viscoelastic effect on the fault sliding and locking behavior can be generated in the long-term stress accumulation and release process. The inversion method used in the application has strong robustness, can well explain the locking behavior under the elastic model, and can still maintain high calculation accuracy in the area with large geological condition change.

Description

Technical Field

[0001] This invention belongs to the field of earthquake prediction technology, specifically relating to a fault finite element analysis method based on a viscoelastic model, and also to computer equipment and storage media, which is suitable for calculating and inverting the degree of fault blocking. Background Technology

[0002] Fault blocking refers to the phenomenon where rock blocks on either side of a fault cannot slide relative to each other due to factors such as friction and stress accumulation. Accurate assessment of the degree of fault blocking helps reveal sliding mechanisms, seismic activity patterns, and focal characteristics, thus providing fundamental data for earthquake prediction and disaster prevention and mitigation. Especially in seismically active areas (such as the Himalayas), studying fault blocking can reveal stress accumulation and release processes, providing crucial support for regional seismic hazard assessment.

[0003] Existing methods mainly rely on elastic models and simplified fault geometry assumptions, which are difficult to accurately reflect the actual situation. Elastic models fail to consider the viscous properties of faults, and simplified geometric models cannot accurately describe the complex morphology of faults, resulting in significant biases in the assessment results. Summary of the Invention

[0004] The purpose of this invention is to address the aforementioned problems in the prior art by providing a tomographic finite element analysis method based on a viscoelastic model, as well as providing computer equipment and storage media.

[0005] The above-mentioned objectives of the present invention are achieved by the following technical means:

[0006] The finite element method for fault analysis based on the viscoelastic model includes the following steps:

[0007] Step 1: Obtain the initial fault geometry data and preprocess the initial fault geometry data to obtain the corrected fault geometry data;

[0008] Step 2: Construct a fault geometric model based on the corrected fault geometric data;

[0009] Step 3: Configure a viscoelastic material model for the fault geometry model in the finite element solver to obtain a three-dimensional viscoelastic finite element model of the target fault region, and set the initial physical boundary conditions of the viscoelastic material model.

[0010] Step 4: Obtain GPS observation data from each GPS station in the target fault area, extract surface displacement field observation data from each GPS observation data, and construct observation constraints based on the GPS observation data.

[0011] Step 5: Calculate the Green's function for the reverse slip rate and surface displacement, and the strike slip rate and surface displacement, respectively. Construct the objective function for inversion. Use the surface displacement field observation data as the input for inversion. By minimizing the objective function, obtain the optimal slip rate distribution of the current three-dimensional viscoelastic finite element model.

[0012] Step 6: Use the grid search method to search for different viscosity parameters and obtain the backslip rate distribution of different viscosity parameters. Based on the backslip rate distribution of each viscosity parameter, obtain the corresponding predicted surface displacement field data. Compare the predicted surface displacement field data of each viscosity parameter with the observed surface displacement field data to select the optimal viscosity parameter, thereby obtaining the optimal three-dimensional viscoelastic finite element model.

[0013] As described above, step 1 specifically includes the following steps:

[0014] Step 1.1: Obtain initial fault geometry data, initial fault depth data, surface undulation information, and subsurface stratification model for the target fault region. Preprocess the initial fault geometry data, including the following steps:

[0015] Step 1.2: Based on the initial fault geometry data and initial fault depth data of the target fault region, the geometric structure of the target fault region is corrected using seismic wave imaging technology to obtain fault geometry correction data;

[0016] Step 1.3: Smooth and transform the fault geometric correction data, and then transform the fault geometric correction data;

[0017] Step 1.4: Perform coordinate rotation and resampling on the converted fault geometric correction data to obtain fault geometric data;

[0018] Step 1.5: Extract elevation data and the projection of the fault on the ground surface from the surface undulation information. Obtain the geometric shape of the fault at different depths underground based on the underground stratification model. Correct the fault geometric data to obtain the corrected fault geometric data.

[0019] As described above, step 2 specifically involves: inputting the corrected fault geometry data into the finite element mesh generator, performing three-dimensional mesh modeling on the target fault region, generating a mesh file and exporting it; reading the exported fault geometry model mesh file through the finite element solver, assigning material parameters to the initial fault geometry model in the finite element solver, and obtaining the fault geometry model.

[0020] As described above, step 3 specifically includes the following steps:

[0021] Step 3.1: In the finite element solver, configure a viscoelastic material model for the fault geometry model. The viscoelastic material model selected is the Maxwell body material model.

[0022] Step 3.2: Set the initial physical boundary conditions in the finite element solver. The initial physical boundary conditions include the initial stress distribution of the fault and the fixed no-slip boundary.

[0023] As described above, step 4 specifically includes the following steps:

[0024] Step 4.1: Obtain GPS observation data from multiple GPS observation stations in the target fault area and the uncertainty of the corresponding GPS observation stations;

[0025] Step 4.2: Extract the surface displacement field observation data from each GPS observation data, set the ratio of the reverse slip rate to the strike slip rate, and obtain the surface displacement field sub-data of reverse slip and strike slip from the surface displacement field observation data.

[0026] Step 4.3: Obtain the uncertainty of backflip and slipback and the uncertainty of run-slipback based on the uncertainty of each GPS observation station, and construct the weight matrix of backflip and slipback and the weight matrix of run-slipback based on the uncertainty.

[0027] As described above, step 5 specifically includes the following steps:

[0028] Step 5.1: Based on the established three-dimensional viscoelastic finite element model, calculate the Green's function of the backflip rate and the surface displacement through numerical simulation in the finite element solver.

[0029] Step 5.2: Based on the established three-dimensional viscoelastic finite element model, calculate the Green's function of the slip-back rate and the surface displacement through numerical simulation in the finite element solver.

[0030] Step 5.3: Using the surface displacement field data of thrust-back and strike-slip as input for inversion, the thrust-back rate distribution and strike-slip rate distribution are obtained respectively based on the Green's function of thrust-back rate and surface displacement and the Green's function of strike-slip rate and surface displacement.

[0031] Step 5.4: Construct the objective function. By minimizing the objective function, obtain the optimal backslip rate distribution and run-slip rate distribution of the current three-dimensional viscoelastic finite element model, and then obtain the optimal backslip rate distribution of the current three-dimensional viscoelastic finite element model. The backslip rate is the vector sum of the backslip rate and the run-slip rate.

[0032] As mentioned above, the objective function is based on the following formula:

[0033] ;

[0034] In the formula, Let be the objective function. The Green's function is given by the reverse slip velocity and the surface displacement. Let Green's function be the ratio of strike-slip velocity to surface displacement. The backflip velocity distribution is shown. Let m be the travel-back velocity distribution, and m be the return velocity distribution. For the surface displacement field data of reverse slip, For strike-slip and return-slip surface displacement field data. The weight matrix for the backflip and slippage is... This is the weight matrix for the take-off and return movements. For regularization operators, The regularization coefficient is . Norm symbol, It is a dot product.

[0035] As described above, step 6 specifically includes the following steps:

[0036] Step 6.1: Use the mesh search method to search for different viscosity parameters. Repeat steps 5.1 to 5.4 to obtain the optimal slip rate distribution of the three-dimensional viscoelastic finite element model under different viscosity parameters.

[0037] Step 6.2: Based on the optimal backslip rate distribution of the three-dimensional viscoelastic finite element model under different viscosity parameters obtained in Step 6.1, calculate the predicted surface displacement field data of the three-dimensional viscoelastic finite element model under different viscosity parameters according to the two Green's functions in Steps 5.1 and 5.2 respectively.

[0038] Step 6.3: Compare the predicted surface displacement field data of each viscosity parameter with the observed surface displacement field data, select the optimal viscosity parameter, and thus obtain the optimal three-dimensional viscoelastic finite element model.

[0039] A computer device includes a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps of the fault finite element analysis method based on the viscoelastic model as described above.

[0040] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the fault finite element analysis method based on a viscoelastic model as described above.

[0041] Compared with the prior art, the present invention has the following advantages:

[0042] (1) The inversion method used in this invention has strong robustness and can better explain the locking behavior under the elastic model. It can still maintain high calculation accuracy in areas with large geological conditions (such as the Qinghai-Tibet Plateau).

[0043] (2) In the process of model selection, the present invention introduces a search mechanism of viscous model to obtain the optimal three-dimensional viscoelastic finite element model of the target fault region, which can more comprehensively analyze the response of the fault locking degree to different stress field conditions. This search mechanism helps to understand the evolution characteristics of the fault under the viscous model, especially the important influence that the viscous effect may have on the fault sliding and locking behavior during the long-term stress accumulation and release process.

[0044] (3) Accuracy and efficiency of inversion calculation: The fault-locking inversion method of the present invention can well explain the fault behavior under different three-dimensional viscoelastic finite element models, and shows high advantages in both accuracy and calculation efficiency. By combining elastic and viscous models, the fault-locking mechanism can be understood more comprehensively. The difference between the inversion results and the actual observation data is small, indicating that the method has high adaptability and reliability under different geological conditions. Attached Figure Description

[0045] Figure 1 This is a flowchart of the method of the present invention;

[0046] Figure 2 This is a three-dimensional viscoelastic finite element model of the main faults in the Himalayas, as described in Embodiment 1 of the present invention.

[0047] Figure 3 This is a spatial distribution and GPS fitting map of the degree of closure of major faults in the Himalayas in Embodiment 1 of the present invention;

[0048] Figure 4 This is a GPS prediction fitting map of the Himalayan region in Embodiment 1 of the present invention. Detailed Implementation

[0049] To facilitate understanding and implementation of the present invention by those skilled in the art, the present invention will be further described in detail below with reference to embodiments. The embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0050] Example 1:

[0051] Fault finite element analysis method based on viscoelastic model, such as Figure 1 As shown, it includes the following steps:

[0052] Step 1: Obtain initial fault geometry data and preprocess the initial fault geometry data to obtain preprocessed fault geometry data. This includes the following steps:

[0053] Step 1.1: Obtain initial fault geometry data, initial fault depth data, surface relief information, and subsurface layering model of the target fault region using geological survey data, active tectonic maps, earthquake catalogs (such as focal depth and mechanism solutions), and focal inversion results. Preprocess the initial fault geometry data, including the following steps:

[0054] Step 1.2: Although these preliminary data can reflect the approximate location and geometric configuration of the fault, they often contain certain uncertainties. To improve the accuracy of the fault model during modeling, it is necessary to correct the geometric structure of the target fault region using seismic wave imaging techniques (such as receiver function imaging, reflection / refraction tomography, and teleseismic tomography) based on the initial fault geometry and initial fault depth data of the target fault region. This includes correcting the fault strike distribution, dip distribution, and depth distribution in the geometric structure of the target fault region to obtain more realistic three-dimensional fault geometry correction data, providing high-precision geometric boundary conditions for subsequent numerical simulations.

[0055] Step 1.3: Smooth and transform the fault geometric correction data, and then transform the fault geometric correction data;

[0056] Step 1.4: To improve the spatial accuracy and consistency of the modeling data, coordinate rotation and resampling are performed on the converted fault geometry correction data to align it with the model mesh and reduce errors caused by coordinate system or uneven sampling, thus obtaining fault geometry data.

[0057] Step 1.5: Simultaneously, the converted fault geometry data is processed based on the surface undulation information and the underground stratification model: the elevation data and the projection of the fault on the surface are extracted from the surface undulation information, and the geometric shape of the fault at different underground depths is obtained based on the underground stratification model. This corrects the converted fault geometry data, ensuring that the distribution of the fault on the surface conforms to the actual terrain, and obtaining high-precision corrected fault geometry data.

[0058] Preprocessing the initial fault geometry data laid the foundation for establishing a realistic three-dimensional fault geometry model, thus improving the quality and reliability of the fault geometry model.

[0059] Step 2: Input the preprocessed fault geometry data obtained in Step 1 into Trelis software (finite element mesh generator). Use the mesh generation function of Trelis software to perform three-dimensional mesh modeling of the target fault region, construct a fine mesh structure, generate a mesh file of the fault geometry model suitable for viscoelastic model analysis, and export it to ensure the accuracy and stability of the simulation and ensure that the model can realistically reproduce the morphology and distribution of the fault. Use a finite element solver (in this embodiment, Pylith software is used as the finite element solver) to read the mesh file of the fault geometry model exported from Trelis software. In Pylith software, assign material parameters to the initial fault geometry model. In this embodiment, the material parameters include the elastic modulus, Poisson's ratio, shear modulus, and density of the rock to obtain the fault geometry model.

[0060] Step 3: In Pylith software, based on the fault geometry model established in Step 2, configure a viscoelastic material model for the fault geometry model to obtain a three-dimensional viscoelastic finite element model of the target fault region, and set the initial physical boundary conditions of the viscoelastic material model to provide physical and boundary constraints for the numerical simulation. This specifically includes the following steps:

[0061] Step 3.1: In Pylith software, configure a viscoelastic material model for the fault geometry model generated in Step 2. The viscoelastic material model selected is the Maxwell body material model.

[0062] The Maxwell body is a viscoelastic material model consisting of a spring and a damper connected in series. Its viscous component follows Newtonian fluid properties. In this model, viscous stress is proportional to the strain rate, and the proportionality constant is called the viscosity coefficient. The overall behavior of the material is determined by both elastic response and viscous flow: when the material is subjected to external force, the spring undergoes instantaneous elastic deformation, while the damper induces viscous flow that develops over time. Under constant strain conditions, the stress decays exponentially over time, with the decay rate controlled by the ratio of the viscosity coefficient to the elastic modulus (i.e., relaxation time). If constant stress is maintained, the material will first undergo instantaneous elastic deformation, followed by a viscous flow stage where strain increases linearly with time. The larger the viscosity coefficient, the stronger the material's resistance to flow, the slower the stress relaxation process, and the more likely it is to exhibit viscous-dominated behavior under long-term loads. This model is often used to describe the short-term viscoelastic phenomena of materials such as polymers and biological soft tissues. However, because it only includes a single relaxation time and cannot reflect the steady-state elastic response, it is often necessary to enhance its applicability by extending the model (such as the generalized Maxwell model) in practical applications. Therefore, the Maxwell body material model can reflect the mechanical behavior of faults at different time scales.

[0063] Step 3.2: Set the initial physical boundary conditions in Pylith. The initial physical boundary conditions include the initial stress distribution of the fault and the fixed no-slip boundary to ensure the physical rationality and realism of the simulation results.

[0064] The initial stress distribution of the fault is determined by constructing a stress field model of the target fault region and combining it with the earthquake focal mechanism solution to set the initial normal stress distribution and shear stress distribution of the fault plane.

[0065] Alternatively, by establishing a steady-state slip model of the target fault region, the stress accumulation state of the fault during the inter-seismic phase can be estimated based on the steady-state slip model, thereby obtaining and setting the initial normal stress distribution and shear stress distribution of the fault plane;

[0066] Fixed no-slip boundaries are set in the following way: the bottom of the three-dimensional viscoelastic finite element model is usually set as a fully constrained surface, and the far side (i.e. the boundary away from the target fault area) can be set as a fixed boundary or a semi-fixed boundary according to GPS deformation observation data to approximate the stress boundary of the fault in the real tectonic environment.

[0067] A fixed boundary is defined as a boundary node where the displacement in a specified direction is zero, while other directions allow free movement (e.g., fixing only the horizontal displacement and allowing vertical movement; or fixing only the vertical displacement and allowing horizontal movement).

[0068] Semi-fixed boundaries are those whose movement is restricted by spring damping or friction conditions, rather than being completely fixed, such as elastic or viscous boundaries.

[0069] Step 4: Obtain GPS observation data from each GPS station in the target fault area, extract surface displacement field observation data from each GPS observation data, and construct observation constraints based on the GPS observation data.

[0070] Step 4.1: Obtain GPS observation data from multiple GPS observation stations in the target fault area and the uncertainty of the corresponding GPS observation stations;

[0071] Step 4.2: Based on Euler's laws of motion, the influence caused by factors such as rigid block motion and internal strain rate is removed from the GPS observation data of each GPS observation station to obtain the surface displacement field observation data caused only by fault blocking. The ratio of thrust slip rate to strike slip rate is set, and the surface displacement field sub-data of thrust slip and strike slip slip are obtained from the surface displacement field observation data. The surface displacement field sub-data is used as the input for inversion.

[0072] Step 4.3: To improve the reliability and physical rationality of the inversion results, the uncertainty of GPS observation data is introduced to assign different weights to the surface displacement field observation data of each GPS observation station. Based on the uncertainty of each GPS observation station, the uncertainty of reverse slip and strike slip are obtained respectively, and the weight matrices of reverse slip and strike slip are constructed respectively. The weight of the surface displacement observation data of each GPS observation station is inversely proportional to the uncertainty of the corresponding GPS observation station. That is, the station with smaller uncertainty and higher data accuracy is in the objective function.

[0073] GPS observation uncertainty typically originates from the fitting results of long-term velocity time series observations, reflecting the accuracy of velocity estimation at each station. Common evaluation methods include least-squares fitting residuals, diagonal elements of the covariance matrix, or upper bounds of error estimated using Bayesian methods. During the inversion process, the weight of each observation point is inversely proportional to its uncertainty; that is, stations with lower uncertainty and higher data accuracy have larger weights in the objective function and a more significant impact on the inversion results, while those with lower uncertainty have smaller weights, thus reducing the adverse effects of high-error data on the overall solution. This uncertainty-weighted inversion strategy helps improve the model's response sensitivity to high-quality data, ensuring that the final inverted fault slip behavior not only conforms to the physical mechanisms of earthquakes but also accurately reflects the spatial variation characteristics of surface observation data.

[0074] Step 5: Calculate the Green's functions for the thrust-back rate and surface displacement, and the strike-slip rate and surface displacement, respectively. Construct the objective function for inversion. Use the observed surface displacement field data as the input for inversion. By minimizing the objective function, obtain the optimal thrust-back rate distribution and strike-slip rate distribution of the current three-dimensional viscoelastic finite element model, and thus obtain the optimal backslip rate distribution of the current three-dimensional viscoelastic finite element model. This allows for further understanding of the stress release during fault slip. Specifically, this includes the following steps:

[0075] Step 5.1: Based on the established three-dimensional viscoelastic finite element model, the Green's function of the thrust-back rate and surface displacement is calculated by numerical simulation in Pylith software. The thrust-back rate reflects the fault thrust-back velocity per unit time, revealing the dynamic process of fault activity. The Green's function describes the surface displacement generated by the unit thrust-back rate of each grid block, and can describe the stress release and surface displacement variation law during the thrust-back process.

[0076] Step 5.2: Based on the established three-dimensional viscoelastic finite element model, the Green's function of strike-slip and return rates and surface displacement is calculated by numerical simulation in Pylith software. The Green's function describes the surface displacement generated by the unit strike-slip and return rate of each grid block. It is applied to simulate the surface displacement under different fault strike-slip rates, providing key data required for inversion analysis and further refining the prediction of surface deformation.

[0077] Step 5.3: Using the surface displacement field data of thrust-back and strike-slip as input for inversion, the thrust-back rate distribution and strike-slip rate distribution are obtained respectively based on the Green's function of thrust-back rate and surface displacement and the Green's function of strike-slip rate and surface displacement.

[0078] Step 5.4: Construct the objective function, which is calculated based on the following formula:

[0079] (1);

[0080] in, Let be the objective function. The Green's function is given by the reverse slip velocity and the surface displacement. Let Green's function be the ratio of strike-slip velocity to surface displacement. and These characterize the responses to retrograde and strike slip, respectively. The backflip velocity distribution is shown. The travel-back velocity distribution is given by m, where m is the return velocity of each grid (the return velocity is the vector sum of the backflip return velocity and the travel-back velocity). For the surface displacement field data of reverse slip, For strike-slip and return-slip surface displacement field data. The weight matrix for the backflip and slippage is... Let L be the weight matrix for the take-slip and return-slip operations, and L be the regularization operator. The regularization coefficient is . The norm symbol, It is a dot product.

[0081] By minimizing the objective function, the optimal backslip rate distribution and run-slip rate distribution of the current three-dimensional viscoelastic finite element model are obtained, and then the optimal backslip rate distribution of the current three-dimensional viscoelastic finite element model is obtained. The backslip rate is the vector sum of the backslip rate and the run-slip rate.

[0082] Step 6: Use a grid search method to search for different viscosity parameters (such as effective viscosity coefficient and shear modulus) to obtain the backslip rate distribution for each viscosity parameter. Based on the backslip rate distribution of each viscosity parameter, obtain the corresponding predicted surface displacement field data. Compare the predicted surface displacement field data for each viscosity parameter with the observed surface displacement field data to select the optimal viscosity parameter, thereby obtaining the optimal three-dimensional viscoelastic finite element model. This specifically includes the following steps:

[0083] Step 6.1: Use the mesh search method to search for different viscosity parameters. Repeat steps 5.1 to 5.4 to obtain the optimal backslip rate distribution of the three-dimensional viscoelastic finite element model under different viscosity parameters in order to accurately quantify the fault-locking degree of the three-dimensional viscoelastic finite element model under different viscosity parameters.

[0084] Step 6.2: Based on the optimal backslip rate distribution of the three-dimensional viscoelastic finite element model under different viscosity parameters obtained in Step 6.1, calculate the predicted surface displacement field data of the three-dimensional viscoelastic finite element model under different viscosity parameters according to the two Green's functions in Steps 5.1 and 5.2 respectively.

[0085] Step 6.3: Compare the predicted surface displacement field data of each viscosity parameter with the observed surface displacement field data, select the optimal viscosity parameter, and thus obtain the optimal three-dimensional viscoelastic finite element model.

[0086] The degree of fault locking in the target fault region is characterized by the backslip rate distribution of the optimal three-dimensional viscoelastic finite element model. The magnitude of the backslip rate reflects the difference between the actual slip rate on the fault plane and the regional tectonic loading rate. A backslip rate close to the tectonic rate indicates strong locking in the region, long-term stress accumulation, and high seismic potential; conversely, if the backslip rate is close to zero, it means that the fault releases strain through stable slip, and the seismic hazard is relatively low.

[0087] Furthermore, a quantitative relationship between fault closure degree and seismic hazard can be established. Seismic hazard can be quantified through various indicators, such as the frequency and magnitude of historical earthquakes, near-field stress accumulation rate, backslip intensity distribution in the source area, and earthquake recurrence interval, providing scientific basis and decision support for earthquake disaster prevention and mitigation.

[0088] Taking the main Himalayan fault in the Qinghai-Tibet Plateau region as an example, this invention is applied to the inversion calculation of the actual fault locking degree. The three-dimensional viscoelastic finite element model of the target fault region is as follows: Figure 2 As shown, the locking degree and fitting results obtained from the three-dimensional viscoelastic model inversion are as follows: Figure 3 As shown. Among them, Figure 2The yellow layers represent subducting plates, the green layers represent overlying plates, the red layers represent the lower crust, and the blue and purple layers represent the mantle on both sides. Figure 3 The pink dots represent historical earthquakes, while the red dots represent earthquakes in the last ten years. The higher the degree of closure, the darker the color. Figure 4 The blue arrows represent actual observed surface displacement field data, and the red arrows represent predicted surface displacement field data.

[0089] The fault-locking inversion method used in this invention has strong robustness and can well explain the locking behavior under the elastic model. It can still maintain high calculation accuracy in areas with large geological conditions (such as the Qinghai-Tibet Plateau).

[0090] In the model selection process, this invention introduces a viscous model search mechanism. By searching for different viscosity parameters, the backslip rate distribution of the three-dimensional viscoelastic finite element model under different viscosity parameters is obtained. Furthermore, the predicted surface displacement field data of the three-dimensional viscoelastic finite element model under different viscosity parameters is calculated and compared with actual GPS-observed surface displacement field data to select the optimal combination of viscosity parameters, thus obtaining the optimal three-dimensional viscoelastic finite element model. This allows for a more comprehensive analysis of the response of fault locking degree to different stress field conditions. This search mechanism helps us understand the evolution characteristics of faults under viscous models, especially the significant impact that viscosity effects may have on fault slip and locking behavior during long-term stress accumulation and release.

[0091] Accuracy and efficiency of inversion calculations: The fault-locking inversion method of this invention can well explain the fault behavior under different three-dimensional viscoelastic finite element models, showing high advantages in both accuracy and computational efficiency. By combining elastic and viscous models, the fault-locking mechanism can be understood more comprehensively. The difference between the inversion results and actual observation data is small, indicating that the method has high adaptability and reliability under different geological conditions.

[0092] In one embodiment, a computer device is also provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above method embodiments.

[0093] In one embodiment, a computer-readable storage medium is provided having a computer program stored thereon that, when executed by a processor, implements the steps in the above method embodiments.

[0094] In one embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above method embodiments.

[0095] It should be noted that the embodiments described in this invention are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains can make various modifications or additions to the described embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A finite element method for fault analysis based on a viscoelastic model, characterized in that, Includes the following steps: Step 1: Obtain the initial fault geometry data and preprocess the initial fault geometry data to obtain the corrected fault geometry data; Step 2: Construct a fault geometric model based on the corrected fault geometric data; Step 3: Configure a viscoelastic material model for the fault geometry model in the finite element solver to obtain a three-dimensional viscoelastic finite element model of the target fault region, and set the initial physical boundary conditions of the viscoelastic material model. Step 4: Obtain GPS observation data from each GPS station in the target fault area, extract surface displacement field observation data from each GPS observation data, and construct observation constraints based on the GPS observation data. Step 5: Calculate the Green's function for the reverse slip rate and surface displacement, and the strike slip rate and surface displacement, respectively. Construct the objective function for inversion. Use the surface displacement field observation data as the input for inversion. By minimizing the objective function, obtain the optimal slip rate distribution of the current three-dimensional viscoelastic finite element model. Step 6: Use the grid search method to search for different viscosity parameters and obtain the slip rate distribution of different viscosity parameters. Based on the slip rate distribution of each viscosity parameter, obtain the corresponding predicted surface displacement field data. Compare the predicted surface displacement field data of each viscosity parameter with the surface displacement field observation data to select the optimal viscosity parameter, thereby obtaining the optimal three-dimensional viscoelastic finite element model. Step 1 specifically includes the following steps: Step 1.1: Obtain initial fault geometry data, initial fault depth data, surface undulation information, and subsurface stratification model for the target fault region. Preprocess the initial fault geometry data, including the following steps: Step 1.2: Based on the initial fault geometry data and initial fault depth data of the target fault region, the geometric structure of the target fault region is corrected by seismic wave imaging to obtain fault geometry correction data; Step 1.3: Smooth and transform the fault geometric correction data, and then transform the fault geometric correction data; Step 1.4: Perform coordinate rotation and resampling on the converted fault geometric correction data to obtain fault geometric data; Step 1.5: Extract elevation data and the projection of the fault on the ground surface from the surface undulation information; obtain the geometric shape of the fault at different depths underground based on the underground stratification model; correct the fault geometric data to obtain the corrected fault geometric data. The viscoelastic material model used is the Maxwell body material model.

2. The finite element method for fault analysis based on a viscoelastic model according to claim 1, characterized in that, Step 2 specifically involves: inputting the corrected fault geometry data into the finite element mesh generator, performing three-dimensional mesh modeling on the target fault region, generating a mesh file and exporting it; reading the exported fault geometry model mesh file through the finite element solver, assigning material parameters to the initial fault geometry model in the finite element solver, and obtaining the fault geometry model.

3. The fault finite element analysis method based on the viscoelastic model according to claim 1, characterized in that, Step 3 specifically includes the following steps: Step 3.1: In the finite element solver, configure a viscoelastic material model for the fault geometry model; Step 3.2: Set the initial physical boundary conditions in the finite element solver. The initial physical boundary conditions include the initial stress distribution of the fault and the fixed no-slip boundary.

4. The fault finite element analysis method based on the viscoelastic model according to claim 3, characterized in that, Step 4 specifically includes the following steps: Step 4.1: Obtain GPS observation data from multiple GPS observation stations in the target fault area and the uncertainty of the corresponding GPS observation stations; Step 4.2: Extract the surface displacement field observation data from each GPS observation data, set the ratio of the reverse slip rate to the strike slip rate, and obtain the surface displacement field sub-data of reverse slip and strike slip from the surface displacement field observation data. Step 4.3: Obtain the uncertainty of backflip and slipback and the uncertainty of run-slipback based on the uncertainty of each GPS observation station, and construct the weight matrix of backflip and slipback and the weight matrix of run-slipback based on the uncertainty.

5. The finite element method for fault analysis based on a viscoelastic model according to claim 4, characterized in that, Step 5 specifically includes the following steps: Step 5.1: Based on the established three-dimensional viscoelastic finite element model, calculate the Green's function of the backflip rate and the surface displacement through numerical simulation in the finite element solver. Step 5.2: Based on the established three-dimensional viscoelastic finite element model, calculate the Green's function of the slip-back rate and the surface displacement through numerical simulation in the finite element solver. Step 5.3: Using the surface displacement field data of thrust-back and strike-slip as input for inversion, the thrust-back rate distribution and strike-slip rate distribution are obtained respectively based on the Green's function of thrust-back rate and surface displacement and the Green's function of strike-slip rate and surface displacement. Step 5.4: Construct the objective function. By minimizing the objective function, obtain the optimal backslip rate distribution and run-slip rate distribution of the current three-dimensional viscoelastic finite element model, and then obtain the optimal backslip rate distribution of the current three-dimensional viscoelastic finite element model. The backslip rate is the vector sum of the backslip rate and the run-slip rate.

6. The finite element method for fault analysis based on a viscoelastic model according to claim 5, characterized in that, The objective function is based on the following formula: ； In the formula, Let be the objective function. The Green's function is given by the reverse slip velocity and the surface displacement. Let Green's function be the ratio of strike-slip velocity to surface displacement. The backflip velocity distribution is shown. Let m be the travel-back velocity distribution, and m be the return velocity distribution. This is the surface displacement field data for reverse slip. For strike-slip and return-slip surface displacement field data. The weight matrix for the backflip and slippage is... This is the weight matrix for the take-off and return movements. For regularization operators, The regularization coefficient is . Norm symbol, It is a dot product.

7. The finite element method for fault analysis based on a viscoelastic model according to claim 6, characterized in that, Step 6 specifically includes the following steps: Step 6.1: Use the mesh search method to search for different viscosity parameters. Repeat steps 5.1 to 5.4 to obtain the optimal slip rate distribution of the three-dimensional viscoelastic finite element model under different viscosity parameters. Step 6.2: Based on the optimal backslip rate distribution of the three-dimensional viscoelastic finite element model under different viscosity parameters obtained in Step 6.1, calculate the predicted surface displacement field data of the three-dimensional viscoelastic finite element model under different viscosity parameters according to the two Green's functions in Steps 5.1 and 5.2 respectively. Step 6.3: Compare the predicted surface displacement field data of each viscosity parameter with the observed surface displacement field data, select the optimal viscosity parameter, and thus obtain the optimal three-dimensional viscoelastic finite element model.

8. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the fault finite element analysis method based on the viscoelastic model as described in any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the fault finite element analysis method based on the viscoelastic model as described in any one of claims 1 to 7.

Images