Active fault stability evaluation method and system based on in-situ stress data

Abstract

The application discloses a method and system for evaluating active fault stability based on ground stress data, and relates to the technical field of geological engineering, comprising the following steps: collecting multi-source heterogeneous stress data and auxiliary data and preprocessing the same to obtain preprocessed multi-source heterogeneous stress data and auxiliary data; performing space-time gridding interpolation processing on the preprocessed multi-source heterogeneous stress data and auxiliary data to obtain initial stress field data of a target active fault region; performing optimization reconstruction of the initial stress field data under the constraint of physical laws to output a high-confidence reconstructed stress field and a stress data abnormal intensity distribution map; constructing an active fault dynamics model based on geological and geophysical data of the target active fault region; and based on the stress data abnormal intensity distribution map. The application solves the prediction deviation problem caused by the physical incompatibility of the initial stress field and the solidification management of the model parameters in the traditional method.

Description

Technical Field

[0001] This invention relates to the field of geological engineering technology, and in particular to a method and system for evaluating the stability of active faults based on geostress data. Background Technology

[0002] Existing technologies typically obtain discrete stress data through methods such as hydraulic fracturing and borehole collapse, and use interpolation methods to construct an initial stress field. Then, they combine numerical models to calculate the risk of fault slip. Although measured data can be used, the construction of the initial stress field mainly relies on statistical interpolation and lacks physical constraints, resulting in high uncertainty and physical inconsistency in areas without data.

[0003] Due to the complexity of geological conditions and the limitations of observation methods, measured stress data are often sparsely and unevenly distributed in space. This leads to significant uncertainty in the initial stress field obtained by conventional interpolation methods in areas without data control, and it is difficult to guarantee that it strictly satisfies basic physical laws such as stress balance. This physical inconsistency will further be transmitted to the subsequent dynamic simulation process, affecting the accuracy of model predictions. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a method for evaluating the stability of active faults based on geostress data to solve the problem of physical inconsistency in the initial stress field caused by data sparsity and defects in interpolation methods.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0007] In a first aspect, the present invention provides a method for evaluating the stability of active faults based on geostress data, which includes collecting multi-source heterogeneous stress data and auxiliary data and performing preprocessing to obtain preprocessed multi-source heterogeneous stress data and auxiliary data.

[0008] Spatiotemporal gridding interpolation is performed on the preprocessed multi-source heterogeneous stress data and auxiliary data to obtain the initial stress field data of the target active fault region.

[0009] The initial stress field data is optimized and reconstructed under physical constraints, and a high-confidence reconstructed stress field and stress data anomaly intensity distribution map are output. Based on the geological and geophysical data of the target active fault area, an active fault dynamic model is constructed.

[0010] Based on the stress data anomaly intensity distribution map, anomaly sub-regions in the target active fault region whose stress response behavior deviates from the linear elasticity assumption are identified. According to the amplitude of the stress data anomaly intensity distribution map, the equivalent physical property parameters of the corresponding sub-regions in the active fault dynamic model are corrected to obtain the active fault dynamic model with corrected physical property parameters.

[0011] The newly acquired observation data is fused and updated with the predicted state of the active fault dynamic model through a dynamic assimilation algorithm, and the fault slip probability curve is output.

[0012] Based on the fault slip probability curve and the updated active fault dynamics model, an active fault stability early warning report is generated.

[0013] As a preferred embodiment of the active fault stability evaluation method based on geostress data described in this invention, the method includes the following steps: collecting multi-source heterogeneous stress data and auxiliary data and preprocessing them to obtain preprocessed multi-source heterogeneous stress data and auxiliary data.

[0014] Collect multi-source heterogeneous stress data and auxiliary data in the target active fault region;

[0015] The multi-source heterogeneous stress data and auxiliary data of the target active fault area are processed by coordinate system one and time alignment to obtain preprocessed multi-source heterogeneous stress data and auxiliary data.

[0016] As a preferred embodiment of the active fault stability evaluation method based on geostress data described in this invention, the method involves: performing spatiotemporal gridded interpolation processing on the preprocessed multi-source heterogeneous stress data and auxiliary data to obtain the initial stress field data of the target active fault region, including the following steps:

[0017] Based on the spatial distribution range and time span of the preprocessed multi-source heterogeneous stress data and the stress data in the auxiliary data, a three-dimensional spatial grid and time series covering the target active fault area are established.

[0018] Using the preprocessed multi-source heterogeneous stress data and the stress data in the auxiliary data, the magnitude of the maximum horizontal principal stress, the magnitude of the minimum horizontal principal stress, the direction of the maximum horizontal principal stress, and the direction of the minimum horizontal principal stress are obtained at each grid node of the three-dimensional spatial grid and time series covering the target active fault area through spatial interpolation.

[0019] Based on the spatial clustering characteristics of stress data and prior knowledge of geological structure, an uncertainty measure of the interpolation results at each grid node is obtained.

[0020] The magnitude of the maximum horizontal principal stress, the magnitude of the minimum horizontal principal stress, the direction of the maximum horizontal principal stress, the direction of the minimum horizontal principal stress, and the variance of the interpolation results at each grid node are combined according to the grid structure and time sequence of the three-dimensional spatial grid covering the target active fault region to form the initial stress field data of the target active fault region.

[0021] As a preferred embodiment of the active fault stability evaluation method based on geostress data described in this invention, the method includes the following steps: optimizing and reconstructing the initial stress field data under physical constraints to output a high-confidence reconstructed stress field and stress data anomaly intensity distribution map:

[0022] Based on the physical equilibrium criterion that the stress tensor divergence is zero, the degree to which the stress tensor components in the initial stress field data of the target active fault region conform to the physical laws of stress equilibrium is analyzed.

[0023] Based on the physical equilibrium criterion that the stress tensor divergence is zero, the stress tensor components in the initial stress field data of the target active fault region are evaluated to assess the satisfaction of the physical law of stress equilibrium, and grid nodes that have poor satisfaction of the physical law of stress equilibrium are identified as potential interference terms.

[0024] Based on the physical laws of stress equilibrium, potential interference terms identified in the initial stress field data of the target active fault region are smoothly adjusted so that the adjusted stress tensor components satisfy the stress equilibrium condition within the target active fault region.

[0025] The smoothed stress field data is compared point by point with the initial stress field data of the target active fault area to generate a stress data anomaly intensity distribution map reflecting the adjustment range at each grid node. The smoothed stress field data is then used as a high-confidence reconstructed stress field.

[0026] As a preferred embodiment of the active fault stability evaluation method based on geostress data described in this invention, the active fault dynamic model is constructed based on the geological and geophysical data of the target active fault area, including the following steps:

[0027] Based on the geological and geophysical data of the target active fault area, information on fault geometry, lithological distribution and boundary conditions of the target active fault area is extracted;

[0028] By utilizing the fault geometry, lithological distribution, and boundary condition information of the target active fault region, a parameterized constitutive relation framework describing the stress accumulation and release process in the target active fault region is constructed.

[0029] The parameterized constitutive framework describing the stress accumulation and release process in the target active fault region is combined with the geometric mesh of the target active fault region to form an active fault dynamic model.

[0030] As a preferred embodiment of the active fault stability evaluation method based on geostress data described in this invention, the method includes the following steps: Identifying anomalous sub-regions in the target active fault region whose stress response behavior deviates from the linear elastic assumption based on anomaly intensity distribution maps of stress data.

[0031] Based on the stress data anomaly intensity distribution map, an intensity threshold is set to distinguish between normal and abnormal responses. Areas in the stress data anomaly intensity distribution map with values ​​greater than twice the average standard deviation are marked as candidate anomaly areas.

[0032] Spatial clustering analysis is performed on the marked candidate anomaly regions to merge spatially adjacent candidate anomaly regions with similar anomaly intensities, forming anomaly sub-regions with continuous spatial distribution.

[0033] By comparing the linear elastic stress response theoretical distribution of the anomalous sub-region with the target active fault region, the extent to which the stress response behavior in the anomalous sub-region deviates from the linear elastic assumption is confirmed.

[0034] As a preferred embodiment of the active fault stability evaluation method based on geostress data described in this invention, the method includes the following steps: Based on the amplitude of the anomaly intensity distribution map of the stress data, the equivalent physical property parameters of the corresponding sub-region in the active fault dynamic model are corrected to obtain the active fault dynamic model with corrected physical property parameters.

[0035] Based on the amplitude of the anomaly intensity distribution map of stress data, the influence mode of the anomaly intensity amplitude on the equivalent physical property parameters in the dynamic model of active faults is defined;

[0036] Based on the influence model, the adjustment direction and magnitude of the equivalent physical property parameters of the corresponding anomalous sub-regions in the dynamic model of active faults are determined;

[0037] Based on the adjustment direction and magnitude of the equivalent physical property parameters, the equivalent physical property parameters of the corresponding anomalous sub-regions in the active fault dynamic model are adjusted to obtain the active fault dynamic model with corrected physical property parameters.

[0038] As a preferred embodiment of the active fault stability evaluation method based on geostress data described in this invention, the method includes the following steps: a dynamic assimilation algorithm is used to fuse and update newly acquired observation data with the predicted state of the active fault dynamic model, outputting a fault slip probability curve.

[0039] New observational data were acquired by collecting data from stress monitoring instruments deployed in the target active fault region.

[0040] The newly acquired observation data is transformed into an observation vector with the same dimension as the state variables of the active fault dynamic model by using an ensemble Kalman filter dynamic assimilation algorithm.

[0041] The transformed observation vectors are fused with the predicted state set of the active fault dynamics model to update the posterior probability distribution of the state variables of the active fault dynamics model;

[0042] Based on the posterior probability distribution of the state variables of the updated active fault dynamics model, the probability that the fault slip volume will exceed the preset slip volume threshold within a specified future time period is calculated, and a fault slip probability curve is generated.

[0043] As a preferred embodiment of the active fault stability evaluation method based on geostress data described in this invention, the method includes the following steps: generating an active fault stability early warning report based on the fault slip probability curve and the updated active fault dynamic model state.

[0044] Based on the fault slip probability curve, the stability warning level of active faults is divided. Combining the active fault stability warning level with the spatial distribution of the updated active fault dynamic model state, key warning areas and recommended monitoring frequencies are determined.

[0045] The active fault stability warning level, key warning areas, and recommended monitoring frequency are formatted into a structured active fault stability warning report.

[0046] Secondly, the present invention provides an active fault stability evaluation system based on geostress data, including a preprocessing module for collecting multi-source heterogeneous stress data and auxiliary data and performing preprocessing to obtain preprocessed multi-source heterogeneous stress data and auxiliary data.

[0047] The interpolation module performs spatiotemporal gridded interpolation on the preprocessed multi-source heterogeneous stress data and auxiliary data to obtain the initial stress field data of the target active fault region.

[0048] The module constructs an optimized reconstruction of the initial stress field data under physical constraints, outputs a high-confidence reconstructed stress field and stress data anomaly intensity distribution map, and constructs an active fault dynamic model based on the geological and geophysical data of the target active fault area.

[0049] The correction module identifies anomalous sub-regions in the target active fault region whose stress response behavior deviates from the linear elastic assumption based on the anomaly intensity distribution map of the stress data. According to the amplitude of the anomaly intensity distribution map of the stress data, it corrects the equivalent physical property parameters of the corresponding sub-regions in the active fault dynamic model, and obtains the active fault dynamic model with corrected physical property parameters.

[0050] The fusion module uses a dynamic assimilation algorithm to fuse and update newly acquired observation data with the predicted state of the active fault dynamics model, and outputs the fault slip probability curve.

[0051] The early warning module generates an early warning report on the stability of active faults based on the fault slip probability curve and the updated state of the active fault dynamic model.

[0052] The beneficial effects of this invention are as follows: by collecting and preprocessing multi-source heterogeneous stress data to construct an initial stress field, and then using physical law constraint optimization reconstruction technology to obtain a high-confidence reconstructed stress field and anomaly intensity distribution map, the physical parameters of the active fault dynamic model are adaptively corrected based on anomaly information, and the real-time prediction and early warning report generation of fault slip probability are realized by fusing new observation data through dynamic assimilation algorithm. This solves the prediction deviation problem caused by the physical inconsistency of the initial stress field and the fixed management of model parameters in traditional methods. Attached Figure Description

[0053] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0054] Figure 1 This is a flowchart of a method for evaluating the stability of active faults based on geostress data.

[0055] Figure 2 This is a schematic diagram of an active fault stability evaluation system based on geostress data. Detailed Implementation

[0056] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0057] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0058] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0059] Reference Figures 1-2 As one embodiment of the present invention, this embodiment provides a method for evaluating the stability of active faults based on geostress data, comprising the following steps:

[0060] S1. Collect multi-source heterogeneous stress data and auxiliary data and preprocess them to obtain preprocessed multi-source heterogeneous stress data and auxiliary data.

[0061] S1.1 Collect multi-source heterogeneous stress data and auxiliary data in the target active fault area.

[0062] Furthermore, various geophysical observation methods were utilized, including obtaining absolute stress values ​​using hydraulic fracturing testing devices, obtaining borehole collapse direction using borehole television imaging equipment, and recording seismic waveforms using broadband seismographs to invert focal mechanism solutions. Simultaneously, auxiliary data within the corresponding spatiotemporal range were recorded, including surface deformation information obtained from global navigation satellite system receivers, atmospheric pressure changes recorded by barometers, and theoretical calculations of solid tides, forming a comprehensive dataset containing different spatial scales, different observation principles, and different temporal resolutions.

[0063] Specifically, the limitations of a single data source can be overcome by using multiple complementary sources. For example, hydraulic fracturing data provides point-like absolute stress but is costly, while focal mechanism solution data has a wide coverage but provides relative stress information. Combining the two can enhance the spatial representativeness of the data. Simultaneous acquisition of stress data and auxiliary data reflecting environmental changes provides a foundation for subsequent analysis of environmental driving factors of stress changes, demonstrating a qualitative extension from single stress observation to multi-physics field coupled observation.

[0064] S1.2. The multi-source heterogeneous stress data and auxiliary data of the target active fault area are processed by coordinate system one and time alignment to obtain preprocessed multi-source heterogeneous stress data and auxiliary data.

[0065] Furthermore, the spatial coordinates of all stress data and auxiliary data are uniformly transformed to the same geodetic coordinate system or local engineering coordinate system to eliminate spatial position deviations caused by the use of different coordinate benchmarks; secondly, the timestamps of all data are uniformly aligned to the standard time system, and data with different sampling frequencies are resampled or interpolated to ensure that all data have consistent sampling times in the time series; obvious outliers or missing record segments in the data are marked.

[0066] Specifically, it solves the difficulties in fusion of multi-source data caused by differences in coordinate systems and time bases, unifies spatial coordinates, and emphasizes the importance of time alignment, because stress changes in active faults may have short-period characteristics. Precise time synchronization is the key to capturing transient stress responses, and spatiotemporal dual standardization processing ensures the physical consistency of data fusion.

[0067] S2. Spatiotemporal gridding interpolation is performed on the preprocessed multi-source heterogeneous stress data and auxiliary data to obtain the initial stress field data of the target active fault region.

[0068] S2.1 Based on the spatial distribution range and time span of the preprocessed multi-source heterogeneous stress data and the stress data in the auxiliary data, a three-dimensional spatial grid and time series covering the target active fault area are established.

[0069] Furthermore, the smallest outer envelope cube of all stress data points in space is determined, and the range of this cube is used as the coverage area of ​​the three-dimensional spatial grid. The grid resolution is set according to the research accuracy requirements, for example, the grid is densified near the fault zone and appropriately sparsed in the far field region. In the time dimension, the start and end times of the time series are set based on the time span of all stress data, and the time step is divided according to a fixed time interval, thereby constructing a four-dimensional framework that is regularly discrete in space and sampled at equal intervals in time.

[0070] Specifically, the division of spatial grids is associated with the geological structural characteristics of faults. For example, anisotropic grids are used near known fault traces, with the grid size along the fault strike direction being larger than that perpendicular to the strike direction, in order to adapt to the continuity of the stress field along the fault direction, thereby achieving a balance between efficiency and accuracy.

[0071] S2.2 Using the preprocessed multi-source heterogeneous stress data and the stress data in the auxiliary data, on each grid node of the three-dimensional spatial grid and time series covering the target active fault area, the magnitude of the maximum horizontal principal stress, the magnitude of the minimum horizontal principal stress, the direction of the maximum horizontal principal stress, and the direction of the minimum horizontal principal stress at each grid node are obtained by spatial interpolation.

[0072] Furthermore, for each grid node in the three-dimensional spatial grid and time series covering the target active fault region, all available stress data points are searched within its spatial neighborhood. These stress data points may originate from the absolute stress values ​​provided by hydraulic fracturing tests, the principal stress directions given by borehole collapse statistics, and the regional stress field orientation obtained from source mechanism inversion. For the scalar stress magnitude component, a Kriging interpolation algorithm based on a variogram model is used, considering the spatial autocorrelation of stress data, to provide the optimal linear unbiased estimate for each grid node. For the vector stress direction component, a direction statistical interpolation method on a unit sphere is used, treating the direction data as a unit vector. The interpolation result is obtained by weighted averaging of the unit vector, avoiding the physical meaning errors that may result from directly averaging angle values. For example, when two directions differ by nearly 180 degrees, the arithmetic average will yield an incorrect result perpendicular to the true direction, while the unit vector average can correctly maintain directional continuity. The maximum horizontal principal stress magnitude estimate, minimum horizontal principal stress magnitude estimate, maximum horizontal principal stress direction unit vector, and minimum horizontal principal stress direction unit vector are output for each grid node at the corresponding time step.

[0073] Sparse, multi-source, and heterogeneous stress observation data are transformed into a continuous distribution on a regular spatiotemporal grid, providing a unified input format for subsequent analysis. For the special physical quantity of stress direction, a mathematical processing method different from that of scalar magnitude is adopted, which fully respects the spherical statistical characteristics of directional data. Thus, physically reasonable interpolation results can be obtained even in sparse data regions. For example, when only a few borehole collapse data points indicate an approximately parallel direction, the interpolation results can maintain the trend of that direction, rather than converging to an irrelevant direction.

[0074] S2.3. Based on the spatial clustering characteristics of stress data and prior knowledge of geological structure, the uncertainty measure of the interpolation results at each grid node is obtained.

[0075] Furthermore, spatial clustering analysis was performed on the preprocessed multi-source heterogeneous stress data and the stress data in the auxiliary data. A density-based clustering algorithm was used to identify regional clusters with similar stress characteristics, such as high stress gradient clusters, low stress variability clusters, and stress direction consistency clusters. Each cluster represents a different stress response mode. Combined with the geological structure map of the target active fault area, the tectonic environment category of each grid node was obtained, such as primary fault zone, secondary fault zone, intact rock mass area, and lithological contact zone. For each grid node, its stress clustering characteristics were matched with the geological structure category. The standard kriging interpolation variance is corrected based on the matching results. For example, for grid nodes located in the main fault zone and belonging to the high stress gradient cluster, the uncertainty measure should be multiplied by a structural complexity factor greater than one based on the statistical variance to reflect the additional unpredictability caused by the drastic changes in the stress field within the fault zone. For grid nodes located in the intact rock mass area and belonging to the low variability cluster, the standard variance can be used or even multiplied by a confidence factor less than one to reflect the high reliability of the structurally stable area. The uncertainty measure of the interpolation result on each grid node is output after dual correction by geological structure and data characteristics.

[0076] Specifically, it provides a more geologically accurate uncertainty quantification method that integrates data-driven statistical variance with knowledge-driven structural constraints. For example, in areas with sparse data but simple structures, uncertainty estimation can be reduced based on structural homogeneity, while at fault intersections with dense data but complex structures, a higher level of uncertainty warning is still required. This fusion method significantly improves the geological rationality of uncertainty measurement, enabling subsequent optimization and reconstruction processes to address high-uncertainty areas more effectively.

[0077] S2.4. The magnitude of the maximum horizontal principal stress, the magnitude of the minimum horizontal principal stress, the direction of the maximum horizontal principal stress, the direction of the minimum horizontal principal stress, and the variance of the interpolation results at each grid node are combined according to the grid structure and time sequence of the three-dimensional spatial grid covering the target active fault region to form the initial stress field data of the target active fault region.

[0078] Furthermore, an independent multidimensional array is created for each physical quantity. The array dimension is completely consistent with the number of grid nodes and time steps of the three-dimensional spatial grid covering the target active fault region and the time series, ensuring that each grid node has a corresponding data storage location at each time step. This includes the estimated value of the maximum horizontal principal stress, the estimated value of the minimum horizontal principal stress, the two horizontal components of the unit vector of the maximum horizontal principal stress direction, the two horizontal components of the unit vector of the minimum horizontal principal stress direction, and the corrected interpolation variances of the maximum and minimum principal stresses, which are then filled into the corresponding array elements. These arrays are encapsulated according to a unified spatiotemporal index structure to form a structured data object containing complete stress tensor information and uncertainty information. This data object not only records the best estimate of the stress state but also fully retains the variance information characterizing the reliability of the estimate.

[0079] Specifically, the data structure design fully considers the needs of subsequent physical law constraint optimization and reconstruction. For example, the optimization algorithm can dynamically adjust the strength of the smoothing constraint according to the size of the variance array, apply stronger physical law constraints in regions with large variance to suppress noise, and trust the interpolation results more in regions with small variance.

[0080] S3. Optimize and reconstruct the initial stress field data with physical constraints, and output a high-confidence reconstructed stress field and stress data anomaly intensity distribution map.

[0081] S3.1 Based on the physical equilibrium criterion that the stress tensor divergence is zero, analyze the degree to which the stress tensor components in the initial stress field data of the target active fault region conform to the physical laws of stress equilibrium.

[0082] Furthermore, stress tensor components, including normal stress components and shear stress components, are extracted from the initial stress field data of the target active fault region for each grid node. For each internal grid node, the divergence value of its stress tensor is calculated, which is the sum of the partial derivatives of each component of the stress tensor in the spatial coordinate direction. This calculation can be approximated by the finite difference method. The spatial derivative is estimated by using the stress values ​​of adjacent grid nodes. The stress tensor divergence value is compared with zero to obtain the stress imbalance at each grid node. This imbalance reflects the degree to which the initial stress field data at that location deviates from the static equilibrium condition. The stress imbalance is statistically analyzed for the entire region to generate a spatial distribution map reflecting the degree of stress equilibrium compliance.

[0083] Specifically, the physical rationality of the initial stress field data is quantitatively assessed by discretizing the continuous physical equilibrium differential equations and applying them to the gridded stress field data. This transforms abstract physical laws into calculable numerical indicators. For example, near fault zones, due to sparse data or complex structures, the stress tensor divergence value is often large, indicating that there may be physical inconsistencies in the stress field in this region, providing a clear target area for subsequent optimization and reconstruction.

[0084] S3.2. Based on the physical equilibrium criterion that the stress tensor divergence is zero, evaluate the satisfaction of the stress tensor components in the initial stress field data of the target active fault region with the physical law of stress equilibrium, and identify grid nodes that have poor satisfaction with the physical law of stress equilibrium as potential interference terms.

[0085] Furthermore, based on the stress tensor divergence value at each grid node, a threshold reflecting the deviation of stress balance from an acceptable range is set. This threshold can be determined based on the statistical distribution characteristics of the stress tensor divergence values ​​across the entire region. For example, grid nodes with an absolute divergence value greater than the average plus twice the standard deviation are initially marked as candidate anomalies. Spatial clustering analysis is performed on these candidate anomalies, merging spatially adjacent candidate points with similar divergence values ​​to form anomaly region clusters. Combining the geological structural information of the target active fault area, anomalies located near known structural boundaries, such as the main fault trace, are filtered out because structural boundaries themselves may cause stress concentration, resulting in non-zero divergence. The remaining grid nodes that are not located at the main structural boundary but whose stress tensor divergence values ​​still significantly deviate from zero are formally identified as potential interference terms that poorly satisfy the physical laws of stress balance.

[0086] Specifically, by screening out potential noise or error concentration areas in the data from the perspective of physical laws, and combining pure numerical anomaly detection with geological structural knowledge, we can avoid misjudging reasonable stress concentrations caused by tectonics as data interference. This allows us to more accurately locate the physical inconsistencies that are actually caused by measurement errors or interpolation defects, providing reliable target indicators for subsequent optimization and adjustment.

[0087] S3.3 Based on the physical laws of stress balance, the potential interference terms identified in the initial stress field data of the target active fault region are smoothly adjusted so that the adjusted stress tensor components satisfy the stress balance condition within the target active fault region.

[0088] Furthermore, using the initial stress field data of the target active fault region as initial values, an optimization function is constructed with the goal of minimizing the stress tensor divergence. This function contains two terms: the first term is the sum of squares of the stress tensor divergence, reflecting the degree of satisfaction with the physical equilibrium condition; the second term is the sum of squares of the differences between the adjusted stress field and the initial stress field, reflecting the fidelity of the original data. An iterative optimization algorithm, such as the conjugate gradient method, is used to solve this optimization problem. In each iteration, stronger smoothing constraints are applied to the stress values ​​of the grid nodes corresponding to the identified potential disturbance terms, while maintaining a high fidelity weight for non-disturbance term regions. When the optimization function converges or reaches the maximum number of iterations, the adjusted stress field data is output. This data eliminates physical imbalances in the potential disturbance term regions through smoothing, while preserving the original observation information as much as possible in the normal regions.

[0089] Specifically, while correcting the physical inconsistency problem, the system retains the effective data features to the maximum extent, transforms the physical constraints into an optimizable mathematical objective function, and achieves targeted correction of abnormal regions through differentiated weight settings. For example, in regions with sparse data and potential interference terms, the stress field value is allowed to be adjusted more significantly to prioritize meeting the equilibrium condition, while in normal regions with dense data, the adjustment range is strictly limited to maintain data authenticity. The adaptive adjustment strategy significantly improves the physical rationality and data consistency of the reconstructed stress field.

[0090] S3.4 Compare the smoothed stress field data with the initial stress field data of the target active fault area point by point to generate a stress data anomaly intensity distribution map reflecting the adjustment range at each grid node, and use the smoothed stress field data as a high-confidence reconstructed stress field.

[0091] Furthermore, for each grid node in the three-dimensional spatial grid and time series covering the target active fault region, the differences in stress tensor components between the smoothed stress field data and the initial stress field data of the target active fault region at corresponding locations are obtained, including the differences in the magnitudes of the maximum and minimum horizontal principal stresses, the angle of change of the direction of the maximum and minimum horizontal principal stresses, and the angle of change of the direction of the direction of the minimum horizontal principal stress. These differences are normalized to obtain a scalar value reflecting the overall adjustment range, which is used as the stress data anomaly intensity value at that grid node. The stress data anomaly intensity values ​​of all grid nodes are organized into a two-dimensional or three-dimensional distribution map according to the spatial grid structure, i.e., a stress data anomaly intensity distribution map. High-value areas in the map indicate locations in the initial data that require significant correction to satisfy physical laws. The smoothed stress field data optimized by physical law constraints is officially output as a high-confidence reconstructed stress field. This data ensures that the stress tensor satisfies the equilibrium condition within the region while maintaining the main characteristics of the original observation.

[0092] Specifically, it provides an optimized stress field, generates diagnostic information characterizing the optimization process, and transforms the abstract optimization process into a visualized spatial distribution map. For example, high abnormal intensity values ​​may appear near fault zones, indicating that there is physical inconsistency in the initial data of this area. The high-confidence reconstructed stress field obtained after optimization eliminates unreasonable fluctuations.

[0093] S4. Based on the geological and geophysical data of the target active fault area, construct a dynamic model of the active fault.

[0094] S4.1 Based on the geological and geophysical data of the target active fault area, extract the fault geometry, lithological distribution and boundary condition information of the target active fault area.

[0095] Furthermore, geological maps, structural outline maps, borehole columnar sections, and geophysical exploration profiles of the target active fault area are collected and interpreted. A combination of manual interpretation and computer-aided identification is used to extract geometric parameters such as fault traces, fault attitude, fault segmentation, and fault fracture zone width from these data. Based on regional geological maps, lithological distribution maps, and well logging data, the distribution range and mechanical property differences of different lithological units are identified. The regional tectonic loading direction and rate are obtained from regional tectonic stress field research results, GPS observation data, and far-field plate movement models. These extracted geometric morphology, lithological distribution, and boundary condition information are digitized and gridded to match the three-dimensional spatial grid and time series grid structure covering the target active fault area.

[0096] Specifically, it provides accurate geometric and physical parameter inputs for the dynamic model of active faults, and systematically integrates and quantifies multi-source, multi-scale geological and geophysical information. For example, it transforms qualitative geological descriptions into computable geometric parameters and discrete well logging data into continuous lithological distribution fields, thereby laying a solid foundation for the subsequent construction of a physically meaningful dynamic model and avoiding simulation distortion caused by oversimplification of model parameters in traditional methods.

[0097] S4.2. Using the fault geometry, lithological distribution and boundary condition information of the target active fault region, construct a parameterized constitutive relation framework to describe the stress accumulation and release process in the target active fault region.

[0098] Furthermore, based on the extracted fault geometry information, the target active fault region is divided into fault zone units and surrounding rock units. For fault zone units, a constitutive relation considering the rate-state dependent friction law is adopted, which can describe the intensity evolution, sliding strengthening, and weakening effects during fault slip. For surrounding rock units, a transversely isotropic elastic or elastoplastic constitutive relation is adopted to reflect the mechanical properties of layered rock masses. According to the lithological distribution information, different mechanical parameters are assigned to different lithological units, such as elastic modulus, Poisson's ratio, cohesion, and internal friction angle. Boundary condition information is transformed into external loading conditions for the model, including far-field tectonic stress loading rate and gravity field effect. These constitutive relations, material parameters, and boundary conditions are integrated into a unified mathematical framework to form a parameterized set of dynamic equations. This set of equations can simulate the stress accumulation process in the fault region under external loading and the stress release process when slip occurs.

[0099] Specifically, a physical model that reflects geological reality and is computationally feasible was constructed. A constitutive relation strategy with regional differentiation was adopted. For example, a complex friction law was used in the fault zone to capture the seismic cycle behavior, while a relatively simple constitutive relation was used in the surrounding rock zone to improve computational efficiency. At the same time, through parametric design, the model can adapt to the characteristics of different fault segments by adjusting a few key parameters.

[0100] S4.3. Combine the parameterized constitutive framework describing the stress accumulation and release process in the target active fault region with the geometric mesh of the target active fault region to form an active fault dynamic model.

[0101] Furthermore, the three-dimensional spatial grid covering the target active fault region and the time-series geometric grid are used as the discrete computational grid for the active fault dynamics model, ensuring that the spatial resolution of the active fault dynamics model is consistent with the stress field data. Each element in the grid is assigned a corresponding constitutive relation type and material parameters. Elements located within the fault zone are assigned a rate-state dependent friction constitutive relation, while surrounding rock elements are assigned an elastic or elastoplastic constitutive relation. The material parameter values ​​are determined based on lithological distribution information. Boundary condition information, including displacement boundary conditions and stress boundary conditions, is applied to the boundary nodes of the active fault dynamics model grid. Numerical discretization methods such as the finite element method or finite difference method are used to transform the continuous differential equations in the parameterized constitutive relation framework into a set of algebraic equations on the discrete grid nodes, constructing a complete active fault dynamics model capable of simulating the evolution of stress, strain, and displacement over time at grid nodes under given boundary loading.

[0102] Specifically, it realizes the transformation from a geological conceptual model to a computable numerical model, maintains the consistency between the model grid and the observation data grid, and allows high-confidence reconstructed stress field data to be directly used as the initial conditions of the model. At the same time, the partitioned constitutive relation design enables the model to simultaneously simulate the nonlinear sliding behavior of fault zones and the elastic response of surrounding rocks.

[0103] S5. Based on the abnormal intensity distribution map of stress data, identify abnormal sub-regions in the target active fault region where the stress response behavior deviates from the linear elastic assumption.

[0104] S5.1 Based on the stress data anomaly intensity distribution map, set an intensity threshold to distinguish between normal and abnormal responses, and mark the areas in the stress data anomaly intensity distribution map where the values ​​are greater than twice the average standard deviation as candidate anomaly areas.

[0105] Furthermore, based on the arithmetic mean and standard deviation of the abnormal intensity values ​​of all grid nodes in the stress data anomaly intensity distribution map, the overall level and fluctuation range of the anomaly intensity in the entire region are characterized. The mean plus twice the standard deviation is used as the intensity threshold to distinguish between normal and abnormal responses. This threshold is selected based on statistical principles and can effectively separate outliers that significantly deviate from the background level. Each grid node in the stress data anomaly intensity distribution map is traversed, and nodes with abnormal intensity values ​​greater than the threshold are marked as candidate anomaly points. A preliminary connectivity analysis is performed on spatially adjacent candidate anomaly points to form candidate anomaly regions.

[0106] Specifically, it provides an objective and quantitative standard for identifying abnormal regions. It adaptively determines the threshold by utilizing the statistical characteristics of the data itself, avoiding the incompatibility that may be caused by regional differences due to manually set fixed thresholds. For example, in a high-abnormality background area, the mean and standard deviation are both high, so the threshold is increased accordingly, thus identifying only the most significant anomalies. In a low-background area, it can sensitively detect weak anomalies. The adaptive threshold strategy ensures the reliability of anomaly identification and regional comparability.

[0107] S5.2 Perform spatial clustering analysis on the marked candidate anomaly regions, merge spatially adjacent candidate anomaly regions with similar anomaly intensities, and form anomaly sub-regions with continuous spatial distribution.

[0108] Furthermore, a density-based spatial clustering algorithm is employed. This algorithm can identify densely connected point sets in spatial data. By using labeled candidate anomalies as input data and setting a neighborhood search radius and a minimum point threshold, candidate anomalies that are spatially adjacent are clustered together. Statistical analysis is performed on the anomaly intensity values ​​within each cluster to obtain the average and standard deviation of the anomaly intensity within the cluster. Adjacent clusters with anomaly intensity differences within an acceptable range are merged; for example, two adjacent clusters can be considered for merging if the difference in the average anomaly intensity of the two adjacent clusters is less than the sum of their standard deviations. The algorithm outputs a series of anomaly sub-regions with spatial continuity and similar internal anomaly intensity characteristics.

[0109] This method transforms discrete anomalies into geologically significant continuous anomalies. It combines spatial proximity and physical property similarity as dual criteria for clustering and merging, avoiding the erroneous merging of anomalies with different origins that may occur when clustering is based solely on spatial distance. For example, although high anomaly areas caused by fault zones and high anomaly areas caused by dikes are spatially close, their intensity characteristics may be different. By using intensity similarity constraints, they can be distinguished, thus more accurately reflecting the actual distribution pattern of underground anomaly responses.

[0110] S5.3 Compare the linear elastic stress response theoretical distribution of the anomalous sub-region with the target active fault region to confirm the degree to which the stress response behavior in the anomalous sub-region deviates from the linear elastic assumption.

[0111] Furthermore, based on the geological model and boundary conditions of the target active fault region, the theoretical distribution of stress response in the region under the same external loading is calculated using linear elasticity theory. This yields the spatial distribution of the theoretical maximum and minimum horizontal principal stress magnitudes, directions of the maximum and minimum horizontal principal stresses. The measured stress intensity values ​​within the anomalous sub-regions with continuous spatial distribution are compared with the corresponding theoretical stress response values ​​to obtain the relative deviations. Examples include the ratio of the measured stress gradient to the theoretical stress gradient within the anomalous sub-region, and the difference between the measured principal stress direction deflection angle and the theoretical direction. These deviation indicators are then combined to generate a linear elasticity deviation quantification value for each anomalous sub-region. This value reflects the degree to which the stress response behavior in that region deviates from the ideal elastic behavior.

[0112] Specifically, it provides a means to quantitatively assess the physical nature of abnormal responses by comparing observed anomalies with theoretical expectations, thereby distinguishing deviations of different natures. For example, a high degree of deviation may indicate strong nonlinear behavior such as plastic deformation or damage in the region, while a low degree of deviation may only reflect local lithological changes. Comparative analysis provides a physical basis for subsequent targeted correction of model parameters, avoiding blind adjustments.

[0113] S6. Based on the amplitude of the stress data anomaly intensity distribution map, correct the equivalent physical property parameters of the corresponding sub-region in the active fault dynamic model to obtain the active fault dynamic model after physical property parameter correction.

[0114] S6.1 Based on the amplitude of the abnormal intensity distribution map of stress data, define the influence mode of abnormal intensity amplitude on the equivalent physical property parameters in the dynamic model of active fault.

[0115] Furthermore, the potential correlation between the abnormal intensity amplitude in the stress data anomaly intensity distribution map and key equivalent physical property parameters in the active fault dynamic model, such as equivalent elastic modulus, equivalent viscosity coefficient, and equivalent friction coefficient, is analyzed. A qualitative relationship is established in which an increase in abnormal intensity amplitude corresponds to an adjustment of equivalent physical property parameters towards weakening. For example, a high abnormal intensity amplitude indicates that the rock mass in this area may have damage or nonlinear behavior, so the equivalent elastic modulus should be appropriately reduced. This influence mode is quantified by defining a functional relationship between the abnormal intensity amplitude and the adjustment amount of the equivalent physical property parameters, such as using a linear proportional relationship or a piecewise functional relationship, so that for every unit increase in abnormal intensity amplitude, the equivalent elastic modulus decreases by a certain proportion.

[0116] Specifically, abstract abnormal intensity data are transformed into specific model parameter adjustment instructions, establishing a physical mapping from observed anomalies to model parameter correction. For example, in the high anomaly zone of the fault zone, the pre-slip behavior caused by stress concentration is simulated by reducing the equivalent friction coefficient, while in the high anomaly zone of the surrounding rock, the rock mass damage is reflected by reducing the equivalent elastic modulus. The differentiated influence mode ensures the physical rationality of the correction of the active fault dynamic model.

[0117] S6.2. Based on the influence model, determine the adjustment direction and magnitude of the equivalent physical property parameters of the corresponding anomalous sub-regions in the dynamic model of active faults.

[0118] Furthermore, for each identified anomalous sub-region, the average amplitude of the anomalous intensity distribution map of its internal stress data is obtained; based on the influence mode function, the adjustment amount of the equivalent physical property parameter corresponding to the average amplitude is obtained, such as the percentage reduction of the equivalent elastic modulus; combined with the structural location and linear elastic deviation of the anomalous sub-region, the adjustment amplitude is fine-tuned, for example, a larger adjustment amplitude is given to anomalous sub-regions located in the main fault zone and with high deviation; the adjustment direction and amplitude value of the equivalent physical property parameter corresponding to each anomalous sub-region are output.

[0119] Specifically, it realizes adaptive parameter adjustment based on anomaly characteristics and introduces linear elastic deviation as a correction factor for adjustment magnitude. This avoids the problem of over-correction of minor anomalies or under-correction of severe anomalies that may be caused by relying solely on the magnitude of anomaly intensity. For example, in two regions with similar anomaly intensity magnitudes but large differences in deviation, the region with higher deviation will receive stronger parameter weakening adjustment, thus more accurately reflecting the actual geological conditions.

[0120] S6.3. Based on the adjustment direction and magnitude of the equivalent physical property parameters, the equivalent physical property parameters of the corresponding abnormal sub-regions in the active fault dynamic model are adjusted to obtain the active fault dynamic model after the physical property parameters are corrected.

[0121] Furthermore, the current values ​​of the equivalent physical property parameters of the corresponding anomalous sub-region grid cells are extracted from the active fault dynamics model. Based on the adjustment direction and magnitude, the adjusted parameter values ​​are obtained, for example, by multiplying the current equivalent elastic modulus by one and subtracting the adjustment percentage to obtain a new value. The adjusted equivalent physical property parameter values ​​are then updated to the corresponding grid cells of the active fault dynamics model. A consistency check is performed on the active fault dynamics model to ensure that the parameter changes are within a physically reasonable range, and the active fault dynamics model with corrected physical property parameters is output.

[0122] Specifically, it achieves dynamic and local optimization of the parameters of the active fault dynamic model, maintains the integrity of the active fault dynamic model mesh structure, and only corrects parameters for anomalous areas rather than reconstructing the entire model. This improves the fitting accuracy of the active fault dynamic model while preserving the physical framework and computational efficiency of the original model to the greatest extent.

[0123] S7. Using a dynamic assimilation algorithm, the newly acquired observation data is fused and updated with the predicted state of the active fault dynamic model, and the fault slip probability curve is output.

[0124] S7.1 Acquire new observation data by using stress monitoring instruments deployed in the target active fault area.

[0125] Furthermore, in-situ stress monitoring equipment such as borehole strain gauges, volumetric strain gauges, and hydraulic fracturing test systems deployed in the target active fault area are used to continuously record stress change data at a preset sampling frequency. At the same time, auxiliary data such as atmospheric pressure and theoretical values ​​of solid tides at the corresponding time are also recorded. Preliminary quality control is performed on the collected raw data to remove obvious outliers and the influence of instrument drift, resulting in a high-quality new observation dataset.

[0126] Specifically, it provides real-time and reliable stress change information by using a variety of stress monitoring instruments for complementary observation. For example, borehole strain gauges provide high-frequency response but may be affected by local interference, while hydraulic fracturing tests provide absolute stress values ​​but have a low sampling rate. The combination of the two can simultaneously capture rapid fluctuations and long-term trends in stress changes, providing a rich source of information for subsequent assimilation.

[0127] S7.2. The newly acquired observation data is transformed into an observation vector with the same dimension as the state variables of the active fault dynamic model by using the ensemble Kalman filter dynamic assimilation algorithm.

[0128] Furthermore, observational components related to the state variables of the active fault dynamics model are extracted from the new observational data. For example, components corresponding to the strain field of the model are extracted from borehole strain data, and components corresponding to the stress field of the active fault dynamics model are extracted from hydraulic fracturing data. An observation operator is established to map the state variables of the active fault dynamics model to the observation space. For example, the stress values ​​on the grid nodes of the active fault dynamics model are converted into stress values ​​at the observation points through spatial interpolation. The difference between the new observational data and the values ​​calculated by the observation operator is used to construct an observation vector. The dimension of this vector is consistent with the number of observational data, and each component represents the deviation between the actual observation and the model prediction at the corresponding observation point.

[0129] Specifically, the standard format conversion between heterogeneous observation data and model states is achieved by designing the observation operator to take into account the spatial representativeness and error characteristics of different observation data. For example, for point-like hydraulic fracturing data, the observation operator uses nearest neighbor interpolation, while for regional strain data, it uses weighted average interpolation, thus more accurately reflecting the difference between observation and model prediction.

[0130] S7.3. Fuse the transformed observation vector with the predicted state set of the active fault dynamics model to update the posterior probability distribution of the state variables of the active fault dynamics model.

[0131] Furthermore, the active fault dynamics model is run to predict the current state set from the posterior state set of the previous time step; the mean and covariance matrix of the predicted state set are obtained, as well as the mean and covariance of the observed prediction values ​​after the predicted state is transformed by the observation operator; the Kalman gain matrix is ​​obtained, which quantifies the relative weights of the model prediction and the observed data; the transformed observation vector and the predicted state set are linearly combined using the Kalman gain to obtain the updated posterior state set, whose statistical properties represent the probability distribution of the model state variables after fusing observation information.

[0132] Specifically, it achieves optimal fusion of active fault dynamic model prediction and real-time observation, and uses an ensemble method to approximate the probability distribution of state variables, avoiding the dependence of traditional Kalman filtering on linear models. It can effectively handle the nonlinear behavior of active fault dynamic models, and automatically estimates the prediction error covariance through ensemble statistics without the need for manual setting, thus improving the adaptability and accuracy of the assimilation process.

[0133] S7.4 Based on the posterior probability distribution of the state variables of the updated active fault dynamic model, the probability that the fault slip volume will exceed the preset slip volume threshold within a specified future time period is calculated, and a fault slip probability curve is generated.

[0134] Furthermore, fault-related state variables, such as fault slip rate and fault shear stress, are extracted from the updated posterior state set. Using the current posterior state as initial conditions, the active fault dynamics model is run to predict the slip evolution process of the fault within a specified future time period, obtaining the future slip time series corresponding to each set member. For each future time point, the proportion of all set members whose slip exceeds a preset slip threshold is counted, and this proportion is the probability that the fault will experience significant slip at that time point. These probability values ​​are connected in chronological order to form a fault slip probability curve.

[0135] Specifically, it provides intuitive and quantitative time-series predictions of fault activity risk, and uses the diversity of the posterior state set to characterize model uncertainty. The resulting slip probability curve not only reflects the most likely development trend, but also includes the risk range caused by the uncertainty of the initial state and active fault dynamics model parameters, providing a more comprehensive information basis for disaster prevention decision-making.

[0136] S8. Generate an active fault stability early warning report based on the fault slip probability curve and the updated active fault dynamic model state.

[0137] S8.1. Based on the fault slip probability curve, classify the active fault stability warning level, and combine the active fault stability warning level with the spatial distribution of the updated active fault dynamic model state to determine the key warning areas and recommended monitoring frequency.

[0138] Furthermore, based on the numerical range of the fault slip probability curve, multiple warning level thresholds are set. For example, a slip probability below 1% is classified as a safe level, between 1% and 5% as a concern level, and above 5% as a warning level. The overall warning level is determined based on the fault slip probability value at the current moment. Combining the spatial distribution of the updated active fault dynamic model state, regions with high stress concentration, high slip rate, or friction coefficient close to the critical value are identified in the model. These regions are superimposed with the warning level to determine key warning areas. Based on the distribution range of the warning level and key warning areas, differentiated recommended monitoring frequencies are formulated. For example, in the fault segment at the warning level and located in the key warning area, hourly high-frequency monitoring is recommended, while in the safe level area, daily monitoring frequency can be maintained.

[0139] Specifically, abstract numerical predictions are transformed into actionable disaster prevention guidelines. Probability predictions in the time dimension are coupled with model state distribution in the spatial dimension for analysis. For example, when the overall slip probability is at the level of concern but the model of a certain fault segment shows abnormally high stress concentration, the area can be upgraded to a key early warning area and the monitoring frequency can be increased, thereby achieving precise targeted monitoring of risk areas.

[0140] S8.2 Format the active fault stability early warning level, key early warning areas and recommended monitoring frequency into a structured active fault stability early warning report.

[0141] Furthermore, a standardized report template is created, including fields such as report title, generation time, data source, warning level summary, key area list, and monitoring recommendations; information such as active fault stability warning level, spatial coordinate range description of key warning areas, and recommended monitoring frequency are filled into the corresponding fields; necessary charts and attachments are added, such as fault slip probability curves, key warning area distribution maps, and model stress field cloud maps; the output is in a structured document format, such as Extensible Markup Language (XML) or Portable Document Format (PDO), to ensure that the report content is machine-readable and easy for human review.

[0142] Specifically, it has achieved standardized and automated output of early warning information. The report structure not only includes conclusive information, but also retains the visualization of key intermediate results. For example, the attached model stress field diagram can help professionals understand the physical basis of the early warning conclusions, improve the scientificity and credibility of the report, and provide comprehensive and clear information support for disaster prevention decision-making.

[0143] This embodiment also provides an active fault stability evaluation system based on geostress data, including: a preprocessing module, which collects multi-source heterogeneous stress data and auxiliary data and performs preprocessing to obtain preprocessed multi-source heterogeneous stress data and auxiliary data;

[0144] The interpolation module performs spatiotemporal gridded interpolation on the preprocessed multi-source heterogeneous stress data and auxiliary data to obtain the initial stress field data of the target active fault region.

[0145] The module constructs an optimized reconstruction of the initial stress field data under physical constraints, outputting a high-confidence reconstructed stress field and anomaly intensity distribution map of stress data. Based on the geological and geophysical data of the target active fault area, it constructs a dynamic model of the active fault.

[0146] The correction module identifies anomalous sub-regions in the target active fault region whose stress response behavior deviates from the linear elastic assumption based on the anomaly intensity distribution map of the stress data. According to the amplitude of the anomaly intensity distribution map of the stress data, it corrects the equivalent physical property parameters of the corresponding sub-regions in the active fault dynamic model, and obtains the active fault dynamic model with corrected physical property parameters.

[0147] The fusion module uses a dynamic assimilation algorithm to fuse and update newly acquired observation data with the predicted state of the active fault dynamics model, and outputs the fault slip probability curve.

[0148] The early warning module generates an early warning report on the stability of active faults based on the fault slip probability curve and the updated state of the active fault dynamic model.

[0149] This embodiment also provides a computer device applicable to the active fault stability evaluation method based on geostress data, comprising: a memory and a processor; the memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions to implement the active fault stability evaluation method based on geostress data as proposed in the above embodiment.

[0150] The computer device can be a terminal, comprising a processor, memory, communication interface, display screen, and input devices connected via a system bus. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, carrier networks, NFC (Near Field Communication), or other technologies. The display screen can be an LCD screen or an e-ink screen. The input devices can be a touch layer covering the display screen, buttons, a trackball, or a touchpad on the computer device's casing, or an external keyboard, touchpad, or mouse.

[0151] This embodiment also provides a storage medium storing a computer program that, when executed by a processor, implements the active fault stability evaluation method based on geostress data as proposed in the above embodiments. The storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read Only Memory (EPROM), Programmable Red-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk.

[0152] In summary, this invention constructs an initial stress field by collecting and preprocessing multi-source heterogeneous stress data, and then uses physical law-constrained optimization and reconstruction technology to obtain a high-confidence reconstructed stress field and anomaly intensity distribution map. Based on anomaly information, it adaptively corrects the physical parameters of the active fault dynamic model, and uses a dynamic assimilation algorithm to fuse new observation data to achieve real-time prediction and early warning report generation of fault slip probability. This solves the prediction bias problem caused by the physical inconsistency of the initial stress field and the fixed management of model parameters in traditional methods.

[0153] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for evaluating the stability of active faults based on geostress data, characterized in that: This includes collecting multi-source heterogeneous stress data and auxiliary data and preprocessing them to obtain preprocessed multi-source heterogeneous stress data and auxiliary data; Spatiotemporal gridding interpolation is performed on the preprocessed multi-source heterogeneous stress data and auxiliary data to obtain the initial stress field data of the target active fault region. The initial stress field data is optimized and reconstructed under physical constraints, and a high-confidence reconstructed stress field and stress data anomaly intensity distribution map are output. Based on the geological and geophysical data of the target active fault area, an active fault dynamic model is constructed. Based on the stress data anomaly intensity distribution map, anomaly sub-regions in the target active fault region whose stress response behavior deviates from the linear elasticity assumption are identified. According to the amplitude of the stress data anomaly intensity distribution map, the equivalent physical property parameters of the corresponding sub-regions in the active fault dynamic model are corrected to obtain the active fault dynamic model with corrected physical property parameters. The newly acquired observation data is fused and updated with the predicted state of the active fault dynamic model through a dynamic assimilation algorithm, and the fault slip probability curve is output. Based on the fault slip probability curve and the updated active fault dynamics model, an active fault stability early warning report is generated.

2. The method for evaluating the stability of active faults based on geostress data as described in claim 1, characterized in that: Collect and preprocess multi-source heterogeneous stress data and auxiliary data to obtain preprocessed multi-source heterogeneous stress data and auxiliary data, including the following steps: Collect multi-source heterogeneous stress data and auxiliary data in the target active fault region; The multi-source heterogeneous stress data and auxiliary data of the target active fault area are processed by coordinate system one and time alignment to obtain preprocessed multi-source heterogeneous stress data and auxiliary data.

3. The method for evaluating the stability of active faults based on geostress data as described in claim 2, characterized in that: The preprocessed multi-source heterogeneous stress data and auxiliary data are subjected to spatiotemporal gridding interpolation to obtain the initial stress field data of the target active fault region, including the following steps: Based on the spatial distribution range and time span of the preprocessed multi-source heterogeneous stress data and the stress data in the auxiliary data, a three-dimensional spatial grid and time series covering the target active fault area are established. Using the preprocessed multi-source heterogeneous stress data and the stress data in the auxiliary data, the magnitude of the maximum horizontal principal stress, the magnitude of the minimum horizontal principal stress, the direction of the maximum horizontal principal stress, and the direction of the minimum horizontal principal stress are obtained at each grid node of the three-dimensional spatial grid and time series covering the target active fault area through spatial interpolation. Based on the spatial clustering characteristics of stress data and prior knowledge of geological structure, an uncertainty measure of the interpolation results at each grid node is obtained. The magnitude of the maximum horizontal principal stress, the magnitude of the minimum horizontal principal stress, the direction of the maximum horizontal principal stress, the direction of the minimum horizontal principal stress, and the variance of the interpolation results at each grid node are combined according to the grid structure and time sequence of the three-dimensional spatial grid covering the target active fault region to form the initial stress field data of the target active fault region.

4. The method for evaluating the stability of active faults based on geostress data as described in claim 3, characterized in that: The initial stress field data is optimized and reconstructed under physical constraints to output a high-confidence reconstructed stress field and anomaly intensity distribution map of stress data, including the following steps: Based on the physical equilibrium criterion that the stress tensor divergence is zero, the degree to which the stress tensor components in the initial stress field data of the target active fault region conform to the physical laws of stress equilibrium is analyzed. Based on the physical equilibrium criterion that the stress tensor divergence is zero, the stress tensor components in the initial stress field data of the target active fault region are evaluated to assess the satisfaction of the physical law of stress equilibrium, and grid nodes that have poor satisfaction of the physical law of stress equilibrium are identified as potential interference terms. Based on the physical laws of stress equilibrium, potential interference terms identified in the initial stress field data of the target active fault region are smoothly adjusted so that the adjusted stress tensor components satisfy the stress equilibrium condition within the target active fault region. The smoothed stress field data is compared point by point with the initial stress field data of the target active fault area to generate a stress data anomaly intensity distribution map reflecting the adjustment range at each grid node. The smoothed stress field data is then used as a high-confidence reconstructed stress field.

5. The method for evaluating the stability of active faults based on geostress data as described in claim 4, characterized in that: Based on geological and geophysical data of the target active fault region, a dynamic model of the active fault is constructed, including the following steps: Based on the geological and geophysical data of the target active fault area, information on fault geometry, lithological distribution and boundary conditions of the target active fault area is extracted; By utilizing the fault geometry, lithological distribution, and boundary condition information of the target active fault region, a parameterized constitutive relation framework describing the stress accumulation and release process in the target active fault region is constructed. The parameterized constitutive framework describing the stress accumulation and release process in the target active fault region is combined with the geometric mesh of the target active fault region to form an active fault dynamic model.

6. The method for evaluating the stability of active faults based on geostress data as described in claim 5, characterized in that: Based on the anomaly intensity distribution map of stress data, identify anomalous sub-regions in the target active fault region whose stress response behavior deviates from the linear elastic assumption, including the following steps: Based on the stress data anomaly intensity distribution map, an intensity threshold is set to distinguish between normal and abnormal responses. Areas in the stress data anomaly intensity distribution map with values ​​greater than twice the average standard deviation are marked as candidate anomaly areas. Spatial clustering analysis is performed on the marked candidate anomaly regions to merge spatially adjacent candidate anomaly regions with similar anomaly intensities, forming anomaly sub-regions with continuous spatial distribution. By comparing the linear elastic stress response theoretical distribution of the anomalous sub-region with the target active fault region, the extent to which the stress response behavior in the anomalous sub-region deviates from the linear elastic assumption is confirmed.

7. The method for evaluating the stability of active faults based on geostress data as described in claim 6, characterized in that: Based on the amplitude of the anomaly intensity distribution map of stress data, the equivalent physical property parameters of the corresponding sub-regions in the active fault dynamic model are corrected to obtain the active fault dynamic model with corrected physical property parameters, including the following steps: Based on the amplitude of the anomaly intensity distribution map of stress data, the influence mode of the anomaly intensity amplitude on the equivalent physical property parameters in the dynamic model of active faults is defined; Based on the influence model, the adjustment direction and magnitude of the equivalent physical property parameters of the corresponding anomalous sub-regions in the dynamic model of active faults are determined; Based on the adjustment direction and magnitude of the equivalent physical property parameters, the equivalent physical property parameters of the corresponding anomalous sub-regions in the active fault dynamic model are adjusted to obtain the active fault dynamic model with corrected physical property parameters.

8. The method for evaluating the stability of active faults based on geostress data as described in claim 7, characterized in that: The newly acquired observational data is fused and updated with the predicted state of the active fault dynamic model using a dynamic assimilation algorithm to output the fault slip probability curve. This process includes the following steps: New observational data were acquired by collecting data from stress monitoring instruments deployed in the target active fault region. The newly acquired observation data is transformed into an observation vector with the same dimension as the state variables of the active fault dynamic model by using an ensemble Kalman filter dynamic assimilation algorithm. The transformed observation vectors are fused with the predicted state set of the active fault dynamics model to update the posterior probability distribution of the state variables of the active fault dynamics model; Based on the posterior probability distribution of the state variables of the updated active fault dynamics model, the probability that the fault slip volume will exceed the preset slip volume threshold within a specified future time period is calculated, and a fault slip probability curve is generated.

9. The method for evaluating the stability of active faults based on geostress data as described in claim 8, characterized in that: Based on the fault slip probability curve and the updated active fault dynamics model, an active fault stability early warning report is generated, including the following steps: Based on the fault slip probability curve, the stability warning level of active faults is divided. Combining the active fault stability warning level with the spatial distribution of the updated active fault dynamic model state, key warning areas and recommended monitoring frequencies are determined. The active fault stability warning level, key warning areas, and recommended monitoring frequency are formatted into a structured active fault stability warning report.

10. An active fault stability evaluation system based on geostress data, based on the active fault stability evaluation method based on geostress data according to any one of claims 1 to 9, characterized in that: This includes a preprocessing module, which collects multi-source heterogeneous stress data and auxiliary data and performs preprocessing to obtain preprocessed multi-source heterogeneous stress data and auxiliary data. The interpolation module performs spatiotemporal gridded interpolation on the preprocessed multi-source heterogeneous stress data and auxiliary data to obtain the initial stress field data of the target active fault region. The module constructs an optimized reconstruction of the initial stress field data under physical constraints, outputs a high-confidence reconstructed stress field and stress data anomaly intensity distribution map, and constructs an active fault dynamic model based on the geological and geophysical data of the target active fault area. The correction module identifies anomalous sub-regions in the target active fault region whose stress response behavior deviates from the linear elastic assumption based on the anomaly intensity distribution map of the stress data. According to the amplitude of the anomaly intensity distribution map of the stress data, it corrects the equivalent physical property parameters of the corresponding sub-regions in the active fault dynamic model, and obtains the active fault dynamic model with corrected physical property parameters. The fusion module uses a dynamic assimilation algorithm to fuse and update newly acquired observation data with the predicted state of the active fault dynamics model, and outputs the fault slip probability curve. The early warning module generates an early warning report on the stability of active faults based on the fault slip probability curve and the updated state of the active fault dynamic model.

Images