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

Abstract

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

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of earthquake prediction, and particularly relates to a fault finite element analysis method based on a viscoelastic model, and also relates to a computer device and a storage medium, and is suitable for calculating and inverting fault locking degree. BACKGROUND

[0002] Fault locking refers to the phenomenon that the rock masses on both sides of the fault cannot slide relative to each other due to factors such as friction and stress accumulation. Accurate evaluation of fault locking degree helps to reveal the sliding mechanism, seismic activity pattern and source characteristics, and thus provides basic data for earthquake prediction and disaster prevention and reduction. Especially in the earthquake active region (such as the Himalayan region), studying fault locking can reveal the stress accumulation and release process, and provide key support for regional seismic risk assessment.

[0003] The existing method mainly relies on an elastic model and a simplified fault geometric assumption, and cannot accurately reflect the actual situation. The elastic model fails to consider the viscous properties of the fault, and the simplified geometric model cannot accurately describe the complex morphology of the fault, resulting in large deviation of the evaluation result. SUMMARY

[0004] The present application aims to solve the above problems existing in the prior art, and provides a fault finite element analysis method based on a viscoelastic model, and also provides a computer device and a storage medium.

[0005] The above purpose of the present application is achieved by the following technical means:

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

[0007] Step 1, obtaining fault geometric initial data and pre-processing the fault geometric initial data to obtain corrected fault geometric data;

[0008] Step 2, constructing a fault geometric model according to the corrected fault geometric data;

[0009] Step 3, configuring a viscoelastic material model for the fault geometric model in a finite element solver to obtain a three-dimensional viscoelastic finite element model of a target fault region, and setting initial physical boundary conditions of the viscoelastic material model;

[0010] Step 4, obtaining GPS observation data of each GPS station in the target fault region, extracting surface displacement field observation data from each GPS observation data, and constructing observation constraint conditions according to the GPS observation data;

[0011] Step 5, calculate the Green function of the thrust back slip rate and the surface displacement, the strike back slip rate and the surface displacement respectively, construct the objective function of the inversion, take the surface displacement field observation data as the input of the inversion, and obtain the optimal back slip rate distribution of the current three-dimensional viscoelastic finite element model by minimizing the objective function;

[0012] Step 6, search for different viscous parameters by using the grid search method, respectively obtain the back slip rate distribution of different viscous parameters, obtain the corresponding predicted surface displacement field data according to the back slip rate distribution of each viscous parameter, compare the predicted surface displacement field data of each viscous parameter with the surface displacement field observation data respectively, and screen out the optimal viscous parameter, thereby obtaining the optimal three-dimensional viscoelastic finite element model.

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

[0014] Step 1.1, obtain the fault geometry initial data, fault depth initial data, surface relief information, and underground layered model of the target fault region, and pre-process the fault geometry initial data, which comprises the following steps:

[0015] Step 1.2, correct the geometric structure of the target fault region according to the fault geometry initial data and the fault depth initial data of the target fault region by using the seismic wave imaging technology, and obtain the fault geometry correction data;

[0016] Step 1.3, smooth and coordinate transform the fault geometry correction data, and transform the fault geometry correction data;

[0017] Step 1.4, coordinate rotation and resampling processing are performed on the transformed fault geometry correction data, and the fault geometry data is obtained;

[0018] Step 1.5, extract the elevation data and the projection of the fault on the surface from the surface relief information, obtain the geometric shape of the fault at different depths underground according to the underground layered model, and correct the fault geometry data to obtain the corrected fault geometry data.

[0019] As described above, step 2 is specifically: input the corrected fault geometry data into a finite element mesh generator, perform three-dimensional mesh modeling on the target fault region, generate a mesh file and export it; read the mesh file of the exported fault geometry model through a finite element solver, assign material parameters to the fault geometry initial model in the finite element solver, and obtain the fault geometry model.

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

[0021] Step 3.1, configure a viscoelastic material model for the fault geometry model in the finite element solver, and the viscoelastic material model selects a 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 . The 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 boundaries 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 method of fault finite element analysis based on a viscoelastic model, characterized by, The method comprises the following steps: Step 1, obtaining fault geometry initial data and pre-processing the fault geometry initial data to obtain corrected fault geometry data; Step 2, constructing a fault geometry model according to the corrected fault geometry data; Step 3, configuring a viscoelastic material model for the fault geometry model in a finite element solver to obtain a three-dimensional viscoelastic finite element model of a target fault region, and setting initial physical boundary conditions of the viscoelastic material model; Step 4, obtaining GPS observation data of each GPS station in the target fault region, extracting surface displacement field observation data from each GPS observation data, and constructing observation constraint conditions according to the GPS observation data; Step 5, calculating the thrust and strike-slip rate and the surface displacement of the thrust and strike-slip rate, respectively, constructing an inversion target function, taking the surface displacement field observation data as the input of the inversion, and obtaining the optimal strike-slip rate distribution of the current three-dimensional viscoelastic finite element model by minimizing the target function; Step 6, searching for different viscous parameters by using a grid search method to obtain the strike-slip rate distribution of each viscous parameter, obtaining the corresponding predicted surface displacement field data according to the strike-slip rate distribution of each viscous parameter, comparing the predicted surface displacement field data of each viscous parameter with the surface displacement field observation data, and screening the optimal viscous parameter to obtain the optimal three-dimensional viscoelastic finite element model.

2. The fault finite element analysis method based on a viscoelastic model according to claim 1, characterized by, The step 1 specifically comprises the following steps: Step 1.1, obtaining fault geometry initial data, fault depth initial data, surface relief information, and a subsurface layered model of a target fault region, and pre-processing the fault geometry initial data, which comprises the following steps: Step 1.2, correcting the geometric structure of the target fault region by a seismic wave imaging technology according to the fault geometry initial data and the fault depth initial data of the target fault region to obtain fault geometry correction data; Step 1.3, smoothing and coordinate conversion are performed on the fault geometry correction data, and the converted fault geometry correction data; Step 1.4, performing coordinate rotation and resampling processing on the converted fault geometry correction data to obtain fault geometry data; Step 1.5, extracting elevation data and the projection of the fault on the ground from the surface relief information, obtaining the geometric shape of the fault at different depths underground according to the subsurface layered model, and correcting the fault geometry data to obtain corrected fault geometry data.

3. The fault finite element analysis method based on a viscoelastic model according to claim 1, characterized by, The step 2 specifically comprises: inputting the corrected fault geometry data into a finite element grid generator, performing three-dimensional grid modeling on the target fault region, generating a grid file and exporting it; reading the grid file of the exported fault geometry model through the finite element solver, and assigning material parameters to the fault geometry initial model in the finite element solver to obtain the fault geometry model.

4. The fault finite element analysis method based on a viscoelastic model according to claim 1, characterized by, The step 3 specifically comprises the following steps: Step 3.1, configuring a viscoelastic material model for the fault geometry model in the finite element solver, and selecting a Maxwell body material model for the viscoelastic material model; Step 3.2, setting initial physical boundary conditions in the finite element solver, and the initial physical boundary conditions include initial stress distribution of the fault and fixed non-slip boundary.

5. The fault finite element analysis method based on a viscoelastic model according to claim 4, characterized by, The step 4 specifically comprises the following steps: Step 4.1, obtaining GPS observation data of a plurality of GPS observation sites in a target fault region and corresponding uncertainties of the GPS observation sites; Step 4.2, extracting surface displacement field observation data from each GPS observation data respectively, setting a ratio of thrust and strike slip rates, and obtaining thrust and strike slip surface displacement field data from the surface displacement field observation data respectively; Step 4.3, obtaining thrust and strike slip uncertainties according to the uncertainties of the GPS observation sites respectively, and constructing a thrust weight matrix and a strike slip weight matrix according to the uncertainties respectively.

6. The fault finite element analysis method based on a viscoelastic model according to claim 5, characterized by, The step 5 specifically comprises the following steps: Step 5.1, calculating a Green function of thrust and surface displacement by numerical simulation in a finite element solver according to the established three-dimensional viscoelastic finite element model; Step 5.2, calculating a Green function of strike and surface displacement by numerical simulation in the finite element solver according to the established three-dimensional viscoelastic finite element model; Step 5.3, taking the thrust and strike slip surface displacement field data as input for inversion, and obtaining thrust and strike slip rate distributions based on the Green functions of thrust and surface displacement and the Green functions of strike and surface displacement respectively; Step 5.4, constructing an objective function, minimizing the objective function, and obtaining optimal thrust and strike slip rate distributions of the current three-dimensional viscoelastic finite element model, and further obtaining optimal slip rate distribution of the current three-dimensional viscoelastic finite element model, wherein the slip rate is a vector sum of the thrust and strike slip rates.

7. The fault finite element analysis method based on a viscoelastic model according to claim 6, characterized by, The objective function is based on the following formula: ； wherein, is the objective function, is the Green's function of reverse-slip rate and surface displacement, is the Green's function of strike-slip rate and surface displacement, is the reverse-slip rate distribution, is the strike-slip rate distribution, m is the slip rate distribution, is the surface displacement field data of reverse-slip, is the surface displacement field data of strike-slip, is the weight matrix of reverse-slip, is the weight matrix of strike-slip, is the regularization operator, is the regularization coefficient, is the norm symbol, is the dot product.

8. The fault finite element analysis method based on a viscoelastic model according to claim 7, characterized by, The step 6 specifically comprises the following steps: Step 6.1, searching for different viscous parameters by using a grid search method, repeating steps 5.1-5.4, and sequentially obtaining optimal slip rate distributions of the three-dimensional viscoelastic finite element model under different viscous parameters; Step 6.2, based on the optimal slip rate distributions of the three-dimensional viscoelastic finite element model under different viscous parameters obtained in step 6.1, calculating predicted surface displacement field data of the three-dimensional viscoelastic finite element model under different viscous parameters according to the two Green functions of steps 5.1 and 5.2; Step 6.3, comparing each viscous parameter predicted surface displacement field data with the surface displacement field observation data respectively, screening out the optimal viscous parameter, and thus obtaining the optimal three-dimensional viscoelastic finite element model. 9.A computer device, comprising a memory and a processor, wherein the memory stores a computer program, and the computer device is configured to perform the method according to any one of claims 1-8 when the computer program is executed by the processor. The processor executes the computer program to realize the steps of the fault finite element analysis method based on the viscoelastic model in any one of claims 1-8.

10. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to realize the steps of the fault finite element analysis method based on the viscoelastic model in any one of claims 1-8.

Images