A three-dimensional in-situ stress numerical analysis method for a geological engineering area

Abstract

The present invention provides a three-dimensional in-situ stress numerical analysis method for a geological engineering area, belonging to the technical field of underground stress analysis. The method includes: conducting geological analysis on the engineering area and its surrounding areas to collect geological data; using the hydraulic fracturing stress measurement method and the hollow inclusion stress relief method to conduct measurements at several measurement points in the engineering area; through analyzing the China Crustal Stress Environment Basic Database and the geological structures and stress field characteristics in the engineering area, deeply analyzing the in-situ stress state of the engineering area; using machine learning methods to obtain three-dimensional stress calibration points with spatial interpolation by using the regional boundary in-situ stress data and the three-dimensional in-situ stress measured data points; using the ASPECT geological modeling software to construct a three-dimensional geological finite element model, applying boundary conditions and loading conditions to the model, and calculating the three-dimensional stress tensor distribution of each point in the model through numerical methods; analyzing the three-dimensional stress field state of the engineering area based on the in-situ stress data and characteristics obtained from the above steps.

Description

Technical Field

[0001] The present invention relates to the technical field of underground stress analysis, and particularly to a three-dimensional in-situ stress numerical analysis method for geological engineering areas. Background Art

[0002] The hydraulic fracturing method is one of the internationally recognized and relatively advanced direct stress measurement techniques in recent years. Its core advantage lies in being able to directly measure the current in-situ stress in underground strata without the need to know in detail the rock mechanical parameters, and being able to obtain multiple stress parameters, with the characteristics of simple operation, continuous or repeated testing at any depth, fast measurement speed, reliable measurement values, etc. Therefore, it has been widely used in recent years and a large number of achievements have been obtained.

[0003] However, the traditional hydraulic fracturing method relies on some assumptions that do not fully conform to the actual geological environment, such as the rock being linearly elastic, isotropic and non-permeable, resulting in limited applicability of the method and only being able to provide limited stress parameters, unable to comprehensively describe the stress field characteristics in complex geological environments. These limitations make the traditional method unable to effectively solve the stress measurement problem under complex geological conditions. Summary of the Invention

[0004] Aiming at the problems existing in the prior art, the present invention provides a three-dimensional in-situ stress numerical analysis method for geological engineering areas. By using the measured point data, considering the influence of terrain, etc., a three-dimensional finite element model is established, and through data assimilation, the three-dimensional stress state of the engineering area is further calculated, and the measured data and the model data are fused to improve the accuracy and reliability of the calculation results.

[0005] To achieve the above object, the present invention provides a three-dimensional in-situ stress numerical analysis method for geological engineering areas, which is characterized by including the following steps:

[0006] Step 1: Conduct a natural geography and engineering geology analysis on the engineering area and its surrounding areas to collect geological data, including the natural geography overview, topography and geomorphology, stratigraphic lithology, regional stability and earthquakes, geological structures and in-situ stress fields, and the influence of engineering geological conditions and hydrogeological conditions on the in-situ stress state of the engineering area;

[0007] Step 2: Use the hydraulic fracturing stress measurement method and the hollow inclusion stress relief method to measure at several measurement points in the engineering area to obtain the three-dimensional stress state of the measurement points;

[0008] Step 3: Through analyzing the China Crustal Stress Environment Basic Database and the geological structures and stress field characteristics in the engineering area, deeply analyze the in-situ stress state of the engineering area;

[0009] Step 4: Using the regional boundary geostress data and the three-dimensional geostress measured data points, a machine learning method is used to obtain the three-dimensional stress calibration points of spatial interpolation;

[0010] Step 5: Based on the data collected and analyzed in steps 1-4, a three-dimensional geological finite element model is constructed using ASPECT geological modeling software, boundary conditions and loading conditions are applied to the model, and the three-dimensional stress tensor distribution of each point in the model is calculated by numerical method;

[0011] Step 6: Analyze the three-dimensional stress field state of the project area based on the geostress data and characteristics obtained in the above steps.

[0012] Optionally, the step 2 specifically includes:

[0013] In-situ stress measurements of hydraulic fracturing were conducted in two boreholes with a designed depth of 90-110 m, and stress relief measurements were conducted near one of the boreholes. The fracturing measurements to determine the stress magnitude were arranged in 18 sections, and the impression measurements to determine the direction of the maximum horizontal principal stress were arranged in 6 sections. The number and location of the measurement points were increased to obtain more measurement data.

[0014] The fracture pressure, instantaneous closure pressure and reopening pressure of the rock are extracted through the pressure-time recording curve to calculate the horizontal maximum and minimum principal stresses and the in-situ tensile strength of the rock.

[0015] Optionally, the step 3 specifically includes:

[0016] A rose diagram of stress direction in the project area was drawn based on the basic database of China's crustal stress environment, showing the direction of the maximum principal stress in the project area and the specific stress state data range;

[0017] By analyzing the basic database of China's crustal stress environment and the geological structure and stress field characteristics in the project area, the direction of the stress field in the project area is determined and more calibration points of the stress field at the model boundary are obtained.

[0018] Optionally, step 4 specifically includes:

[0019] Obtain the regional boundary geostress data and three-dimensional geostress measured data points obtained in the previous step;

[0020] The data is stored in a CSV file. Each line represents a measurement point, including longitude, latitude, and three-dimensional geostress data, which corresponds to the three geostress values. The pandas library is used to read the CSV file and convert the data into a numpy array to obtain the initial geostress data.

[0021] Preliminarily check the outliers and noise points in the data set and remove the data that do not meet the standards;

[0022] Standardize the data using the Min-Max or Z-Score method;

[0023] Install the numpy library, matplotlib library, scikit-learn library, and pykrige library for Python. Use the scikit-learn library to train a Kriging interpolation model based on Gaussian process regression: First, import Gaussian process regression and the RBF kernel function. Then, create a Gaussian process regression model that specifies the use of the RBF kernel function and specifies that when optimizing the model parameters, if a non-global optimal solution is found, restart the optimization process multiple times. Finally, use the model.fit method to train the model;

[0024] After the Kriging interpolation model is trained, use the trained model to predict the values at unobserved locations: Use the numpy library to define a set of coordinates for unobserved locations, where lon represents the longitude of the unobserved location and lat represents the latitude of the observed location. After defining the coordinates of the unobserved locations, use the trained Kriging interpolation model and the predict method to predict the values at the unobserved locations;

[0025] Use the matplotlib library to visualize the observed data and the interpolation results: First, import the matplotlib.pyplot module and name it plt. Then, use the plt.scatter function to plot a scatter plot of the observed data coordinates. Next, use the plt.contourf function to plot a contour map of the data at the unobserved locations obtained by Kriging interpolation. Finally, specify the legend and title and use the plt.show function to generate and display a visualization graph that includes the scatter plot of the observed data and the contour map of the interpolation results.

[0026] Optionally, step 5 specifically includes:

[0027] Based on the geological data obtained previously, as well as the setting of boundary conditions and material fields;

[0028] Use the ASPECT geological modeling software to construct a three-dimensional geological finite element model and solve a series of equations, specifically as follows:

[0029] , (1.1)

[0030] (1.2)

[0031] (1.3)

[0032] (1.4)

[0033] Equation (1.1) represents the compressible Stokes equations, where u = u(x, t) is the velocity field, p = p(x, t) is the pressure field, and both fields depend on the spatial position x and time t. The motion of the substance is driven by the gravitational force acting on the object and is proportional to the density and pressure of the object. η is the dynamic viscosity of the fluid, representing the internal frictional resistance of the fluid. is the strain rate tensor, representing the deformation rate of the velocity field. is the divergence of the velocity field, describing the rate of change of the volume of the fluid. p is the internal pressure of the fluid, ρ is the fluid density, and g is the gravitational acceleration vector.

[0034] Equation (1.2) represents the mass conservation equation. is the divergence of the mass flow rate, representing the conservation of mass per unit volume. ρu is the mass flux.

[0035] Equation (1.3) represents the temperature field equation, which includes a heat conduction term and a velocity field with velocity u. The right-hand side term of this equation also includes the generation of internal heat. is the specific heat capacity at constant pressure, and T is the temperature. is the convective term of the temperature, representing the transport effect of the velocity field on the temperature distribution. k is the thermal conductivity coefficient, describing the heat conduction ability. is the heat conduction term, representing the diffusion of heat. H is the heat source term per unit volume. is the viscous dissipation term, representing the heat generated by the deformation of the fluid. α is the coefficient of thermal expansion, describing the volume expansion caused by temperature changes. ΔS is the entropy change rate.

[0036] Equation (1.4) represents the chemical composition conservation equation, which describes the change of chemical composition over time, including convection and diffusion effects. is the concentration of the i-th chemical component. is the convective term, representing the transport effect of the velocity field on the chemical composition. qi is the chemical reaction or external source term.

[0037] Ω represents the computational domain, that is, the modeled region. in Ω means that this set of equations is solved inside the computational domain Ω.

[0038] When establishing the model, an appropriate grid is divided for the model to ensure the balance between the accuracy of the calculation results and the calculation efficiency. Next is the stress field calculation stage. After completing the model initialization, boundary conditions and loading conditions are applied to the model, and the calculation is run through numerical methods. The main principle is as follows:

[0039] Under the action of gravitational force, inertial centrifugal force, Coriolis force, and pressure gradient force, the motion equation is:

[0040] (1.5)

[0041] The equations of motion are as follows:

[0042] (1.6)

[0043] The continuity equation is:

[0044] (1.7)

[0045] where V is the velocity vector, G is the effective gravity, including the resultant force of the earth's gravitational force and the inertial centrifugal force, Ω is the earth's angular velocity vector of rotation, r is the position vector, ρ is the fluid density, is the pressure gradient, is other external forces, ν is the kinematic viscosity coefficient, ΔV is the Laplacian operator of velocity, representing the viscous diffusion term, w 、 u 、 v are the components of velocity in z 、 x 、 y directions, x , y , z are the spatial coordinates, is the component of the Coriolis parameter, g is the gravitational acceleration, is the horizontal turbulent diffusion coefficient;

[0046] Finally, the three-dimensional stress tensor distribution at each point in the model is solved.

[0047] After adopting the above technical solution, the present invention has at least the following beneficial effects:

[0048] By introducing a three-dimensional finite element model and data assimilation technology, the present invention significantly optimizes the application of the traditional hydraulic fracturing method under complex geological conditions, can provide more accurate and reliable underground stress measurement results, can combine multi-dimensional data, and through detailed stress field simulation, overcomes the limitations of traditional methods, especially showing unique advantages in complex environments with multi-level and multi-direction stress fields, overcomes the shortcoming that traditional methods can only provide limited stress parameters, not only improves the measurement accuracy, but also provides more accurate data support for engineering design, and has significant economic benefits and risk control advantages. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

[0050] Figure 1 It is a schematic flowchart of a three-dimensional in-situ stress numerical analysis method for a geological engineering area provided by an embodiment of the present disclosure;

[0051] Figure 2 It is a pressure-time recording curve graph for hydraulic fracturing stress measurement;

[0052] Figure 3 It is a longitudinal section position diagram of three-dimensional visualization of the engineering area;

[0053] Figure 4 It is a three-dimensional in-situ stress distribution diagram of the engineering area. Specific Embodiments

[0054] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, rather than all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts belong to the scope of protection of the present invention.

[0055] The following will detail the solution of the present invention with specific embodiments:

[0056] Reference Figure 1 , a three-dimensional in-situ stress numerical analysis method for a geological engineering area provided by an embodiment of the present disclosure, characterized by including the following steps:

[0057] Step 1: Conduct a natural geography and engineering geology analysis on the engineering area and its surrounding areas to collect geological data, including the natural geography overview, topography, stratigraphic lithology, regional stability and earthquakes, geological structures and in-situ stress fields, and the influence of engineering geological conditions and hydrogeological conditions on the in-situ stress state of the engineering area;

[0058] Step 2: Use the hydraulic fracturing stress measurement method and the hollow inclusion stress relief method to measure at several measurement points in the engineering area to obtain the three-dimensional stress state of the measurement points;

[0059] In this step, in-situ stress measurement by hydraulic fracturing is carried out in two drill holes with a designed hole depth of 90 - 110 m. Stress relief measurement is carried out near one of the drill holes. For the fracturing measurement arrangement to determine the stress magnitude, 18 segments are set, and for the impression measurement to determine the direction of the maximum horizontal principal stress, 6 segments are set. By increasing the number and location of measurement points, more measurement data can be obtained. According to the measurement data, a pressure-time recording curve is plotted. As Figure 2 shown, the fracture pressure of the rock is extracted from the pressure-time recording curve , the instantaneous closure pressure and the reopening pressure . The fracture pressure is obtained by determining the peak pressure of the first circulation round trip. In the hydraulic fracturing test, as the fluid pressure gradually increases, when the rock reaches its tensile strength limit, cracking will occur, and the pressure peak recorded at this time is the fracture pressure. The reopening pressure takes the points with obvious slope changes on the pressure-time curve. Usually, the pressure value of the third circulation round trip or the average value of the second to fourth round trips is used. The instantaneous closure pressure is equal to the minimum horizontal principal stress. The measurement methods include the inflection point method and the single tangent method, etc. The inflection point method determines the instantaneous closure pressure by observing the inflection point on the pressure-time curve, while the single tangent method approximates the position of the instantaneous closure pressure by drawing a tangent. The single tangent method is commonly used to ensure accuracy. The fracture pressure , the instantaneous closure pressure and the reopening pressure are used to calculate the maximum and minimum horizontal principal stresses and the in-situ tensile strength of the rock.

[0060] Step 3: By analyzing the China Crustal Stress Environment Basic Database and the geological structure and stress field characteristics in the project area, deeply analyze the in-situ stress state of the project area;

[0061] In this step, in order to more deeply explore the characteristics and distribution laws of the stress field in the project area, it is necessary to carefully consult and analyze the China Crustal Stress Environment Basic Database, draw a rose diagram of the stress direction, and show the direction of the maximum horizontal principal stress in the region and the specific stress state data range.

[0062] By analyzing the China Crustal Stress Environment Basic Database and the geological structure and stress field characteristics in the project area, the aim is to determine the direction of the regional stress field, obtain more calibration points of the stress field at the model boundary. This comprehensive analysis method helps to reveal the dynamic changes of the regional stress field and provides an important basis for engineering design. The main data includes the following:

[0063] Focal mechanism solution data: also known as fault plane solution, is an important parameter used to study the force at the source and the nature of fault dislocation when an earthquake occurs using seismic observation data. It contains key information such as fault strike, fault dip, and slip angle. These parameters help to understand the type of fault, reveal the specific movement of the fault when an earthquake occurs, and describe the characteristics of the slip plane.

[0064] Fault slip data: refers to data that records the relative sliding of rock blocks on both sides of the fault along the fracture surface. These data usually include information such as the type of fault (such as normal fault, reverse fault, strike-slip fault, etc.), the direction and distance of the sliding, etc.

[0065] Hydraulic fracturing data: data obtained through hydraulic fracturing experiments. Hydraulic fracturing is a method of sealing both ends of a certain section of a borehole, and then injecting high-pressure water into the sealed section to rupture the borehole wall and measure the ground stress. The data recorded during the experiment include the water pressure at the initial cracking, the water pressure changes during the crack expansion process, etc.

[0066] Stress relief data: data obtained through stress relief experiment. Stress relief experiment is a method to estimate the original stress by measuring the deformation or strain of rock or soil during stress relief. The data recorded during the experiment include deformation, strain value, etc.

[0067] Borehole collapse data: data obtained through borehole collapse experiments. Borehole collapse experiments are a method of studying crustal stress by observing the collapse of the borehole wall. The data recorded during the experiment include the depth, direction, and shape of the collapse.

[0068] By obtaining the principal stress square and stress contact data of the region, it can be verified with the hydraulic fracturing in situ stress measurement data obtained in step 2. If the hydraulic fracturing data and the stress relief data are consistent in size and direction, it can be considered that the data has a small fluctuation range in the region; conversely, if the hydraulic fracturing data and the stress relief data are significantly different in size and direction, there may be stimulating fractures in the region.

[0069] Step 4: Using the regional boundary geostress data and the three-dimensional geostress measured data points, a machine learning method is used to obtain the three-dimensional stress calibration points of spatial interpolation;

[0070] In this step, since only a limited number of three-dimensional geostress data points can be obtained in actual measurements, in order to fill these data gaps, site Kriging interpolation (Kriging) can be used. It is a commonly used spatial interpolation method that uses data at measurement points to calculate the covariance function between known point pairs to describe spatial correlation, and estimates the values ​​of unknown points based on these correlations. It can produce interpolation results that are closer to the actual situation and estimate the values ​​at unobserved locations by modeling spatial autocorrelation.

[0071] The specific process is as follows:

[0072] Obtain the regional boundary geostress data and three-dimensional geostress measured data points obtained in the previous step;

[0073] The data is stored in a CSV file. Each line represents a measurement point, including longitude, latitude, and three-dimensional geostress data, which corresponds to the three geostress values. The pandas library is used to read the CSV file and convert the data into a numpy array to obtain the initial geostress data.

[0074] Preliminarily check the outliers and noise points in the data set and remove the data that do not meet the standards;

[0075] Use the Min-Max or Z-Score method to standardize the data to meet the input requirements of the machine learning model and avoid the influence of features of different dimensions on model training;

[0076] Install the numpy library, matplotlib library, scikit-learn library, and pykrige library for Python. The numpy library is used for data processing and calculation, the matplotlib library is used for data visualization, the scikit-learn library is used for the implementation of machine learning methods, and the pykrige library is used for the implementation of Kriging interpolation.

[0077] Use the scikit-learn library to train a kriging interpolation model based on Gaussian process regression: first import the Gaussian process regression and RBF variogram, then create a Gaussian process regression model that specifies the use of RBF variogram, and specify when optimizing model parameters, restart the optimization process multiple times if a non-global optimal solution is found, and finally use the model.fit method to train the model;

[0078] After the Kriging interpolation model training is completed, use the trained model to predict the values at unobserved locations: Define a set of coordinates of unobserved locations using the numpy library, where lon represents the longitude of the unobserved location and lat represents the latitude of the observed location. After completing the definition of the coordinates of the unobserved locations, use the trained Kriging interpolation model and the predict method to predict the values at the unobserved locations;

[0079] Use the matplotlib library to visualize the observed data and the interpolation results: First, import the matplotlib.pyplot module and name it plt. Then use the plt.scatter function to draw a scatter plot of the coordinates of the observed data. Next, use the plt.contourf function to draw a contour map of the data at the unobserved locations obtained by Kriging interpolation. Finally, after specifying the legend and title, use the plt.show function to generate and display a visualization graph that includes the scatter plot of the observed data and the contour map of the interpolation results.

[0080] In this step, combining geographic information systems, geological engineering, and machine learning techniques, an efficient three-dimensional in-situ stress spatial interpolation method is proposed. Using Python's powerful data processing and machine learning libraries, especially the Gaussian process regression module in the scikit-learn library, point Kriging interpolation is realized, effectively filling the spatial gap in in-situ stress data. Through data visualization techniques, the comparison between the interpolation results and the observed data is intuitively demonstrated, verifying the accuracy and practicality of the method. This method is not only applicable to the interpolation of in-situ stress data but also provides new ideas and technical paths for the analysis and processing of other geospatial data.

[0081] Step 5: Based on the data collected and analyzed in Steps 1 - 4, use the ASPECT geological modeling software to construct a three-dimensional geological finite element model, apply boundary conditions and loading conditions to the model, and calculate the three-dimensional stress tensor distribution of each point in the model through numerical methods;

[0082] In this step, obtain the geological data, boundary conditions, and the setting of the material field obtained previously. The material field component is a mathematical model used to simulate the interaction relationships between different materials and is usually used to track the movement process and changes of materials. Calculating the increment of the material field usually needs to be carried out according to the reaction rate q(T, c), which involves finding an equilibrium state between different domains to enable chemical reactions to occur;

[0083] Use the ASPECT geological modeling software to construct a three-dimensional geological finite element model and solve a series of equations as follows:

[0084] , (1.1)

[0085] (1.2)

[0086] (1.3)

[0087] (1.4)

[0088] Equation (1.1) represents the compressible Stokes equation, where u = u(x, t) is the velocity field, p = p(x, t) is the pressure field, both fields depend on the spatial position x and time t, the motion of the substance is driven by the gravity acting on the object and is proportional to the density and pressure of the object, η is the dynamic viscosity of the fluid (unit: Pa·s), representing the internal frictional resistance of the fluid, is the strain rate tensor (unit: ), representing the deformation rate of the velocity field, is the divergence of the velocity field (unit: ), describing the rate of change of the volume of the fluid, p is the pressure (unit: Pa), the internal pressure of the fluid, ρ is the fluid density (unit: kg / m³), g is the gravitational acceleration vector (unit: m / s²), driving the motion of the fluid;

[0089] Equation (1.2) represents the mass conservation equation, is the divergence of the mass flow rate (unit: kg / (m³•s)), representing the conservation of mass per unit volume, ρu is the mass flux (unit: kg / (m²•s)), the product of density and velocity;

[0090] Equation (1.3) represents the temperature field equation, which includes a heat conduction term and a velocity field with velocity u, and the right - hand side term of the equation also includes the generation of internal heat, for example, heat generation by radioactive decay, frictional heating, adiabatic compression of materials, heat generation by phase change, etc., is the specific heat capacity at constant pressure (unit: J / (kg•K)), T is the temperature (unit: K), is the convective term of temperature (unit: K / s), representing the transport effect of the velocity field on the temperature distribution, k is the thermal conductivity (unit: W / (m•K)), describing the heat conduction ability, is the heat conduction term (unit: W / m³), representing the diffusion of heat, H is the heat source term per unit volume (unit: W / m³), is the viscous dissipation term (unit: W / m³), representing the heat generated by the deformation of the fluid, α is the coefficient of thermal expansion (unit: ), describing the volume expansion caused by temperature change, ΔS is the entropy change rate (unit: J / (kg•K•s)), related to phase change or chemical reaction;

[0091] Equation (1.4) represents the conservation equation of chemical composition, describing the change of chemical composition over time, including convection and diffusion effects. is the concentration of the i-th chemical component (unit: mol / m³ or mass fraction). is the convection term (unit: mol / (m³•s)), representing the transport effect of the velocity field on the chemical composition, and qi is the chemical reaction or external source term (unit: mol / (m³•s)), such as phase change, diffusion or external input.

[0092] Ω represents the computational domain, that is, the modeled region, and "in Ω" indicates that this set of equations is solved inside the computational domain Ω.

[0093] When establishing the model, an appropriate grid is divided for the model to ensure the balance between the accuracy of the calculation results and the calculation efficiency. Next is the stress field calculation stage. After completing the model initialization, boundary conditions and loading conditions are applied to the model, and the calculation is run through numerical methods. The main principles are as follows:

[0094] Under the action of gravitational force, inertial centrifugal force, Coriolis force, and pressure gradient force, the motion equation is:

[0095] (1.5)

[0096] The system of motion equations is:

[0097] (1.6)

[0098] The continuity equation is:

[0099] (1.7)

[0100] Among them, V is the velocity vector (unit: m / s), G is the effective gravity (unit: m / s²), including the resultant force of gravitational force and inertial centrifugal force, Ω is the angular velocity vector of the Earth's rotation (unit: rad / s), r is the position vector (unit: m), ρ is the fluid density (unit: kg / m³), is the pressure gradient (unit: Pa / m), is other external forces (such as thermal forcing, unit: N / kg), ν is the kinematic viscosity coefficient (unit: m² / s), ΔV is the Laplacian operator of velocity (unit: m / s³), representing the viscous diffusion term. w 、 u 、 v are the velocities in z 、 x 、 yComponent in the direction (unit: m / s), x , y , z are spatial coordinates (unit: m), is the component of the Coriolis parameter (unit: ), g is the gravitational acceleration (unit: m / s²), is the horizontal turbulent diffusion coefficient (unit: m² / s), characterizing the turbulent mixing effect;

[0101] Finally, the three-dimensional stress tensor distribution at each point in the model is solved.

[0102] Step 6. Analyze the three-dimensional stress field state in the engineering area based on the in-situ stress data and characteristics obtained in the above steps.

[0103] In this step, according to the calculated stress field distribution, analyze the geological stability, evaluate the risks of geological disasters such as landslides and rock bursts, and use the results to guide engineering design, such as key works like tunnel excavation and mine exploitation.

[0104] Use the method described in the above embodiment to conduct three-dimensional in-situ stress analysis on the measured in-situ stress data of an engineering area, and obtain the three-dimensional in-situ stress analysis results of the engineering area. Figure 3 shows the longitudinal section position of the engineering area, which can clearly display the spatial layout of the strata and the overall in-situ stress and important structural features, providing an intuitive reference for subsequent analysis. Figure 4 shows the three-dimensional in-situ stress distribution of the engineering area. The spatial distribution of the in-situ stress data is more effectively displayed, enabling the accuracy of the three-dimensional in-situ stress distribution to be improved to the order of 10 meters, enhancing the credibility and application value of the calculation results. At the same time, it effectively eliminates the deviation between the model and reality, enhances the stress prediction ability in local areas, and provides a more scientific basis for subsequent engineering design, structural optimization, etc.

[0105] Although the present invention has been disclosed as above with embodiments, it is not intended to limit the present invention. Any person skilled in the art within the technical field can make some changes and modifications without departing from the spirit and scope of the present invention. Therefore, the protection scope of the present invention shall be subject to what is defined by the claims.

Claims

1. A three-dimensional geostress numerical analysis method for geological engineering regions, characterized in that: The following steps are involved: Step 1: Conduct physical geography and engineering geology analysis of the project area and its surrounding areas to collect geological data, including physical geography, topography, stratum lithology, regional stability and earthquakes, geological structure and geostress field, engineering geological conditions and hydrogeological conditions, and their impact on the geostress state of the project area; Step 2: Using a hydraulic fracturing stress measurement method and a hollow inclusion stress relief method to measure at several measurement points in the project area, and obtain a three-dimensional stress state of the measurement points; Step 3: By analyzing the basic database of China's crustal stress environment and the geological structure and stress field characteristics in the project area, we can deeply analyze the ground stress state of the project area; Step 4: Using the regional boundary geostress data and the three-dimensional geostress measured data points, a machine learning method is used to obtain the three-dimensional stress calibration points of spatial interpolation; Step 5: Based on the data collected and analyzed in steps 1-4, a three-dimensional geological finite element model is constructed using ASPECT geological modeling software, boundary conditions and loading conditions are applied to the model, and the three-dimensional stress tensor distribution of each point in the model is calculated by numerical method; Step 6: Analyze the three-dimensional stress field state of the project area based on the geostress data and characteristics obtained in the above steps.

2. The method for numerical analysis of three-dimensional geostress in a geological engineering region according to claim 1, characterized in that: The step 2 specifically includes: In-situ stress measurements of hydraulic fracturing were conducted in two boreholes with a designed depth of 90-110 m, and stress relief measurements were conducted near one of the boreholes. The fracturing measurements to determine the stress magnitude were arranged in 18 sections, and the impression measurements to determine the direction of the maximum horizontal principal stress were arranged in 6 sections. The number and location of the measurement points were increased to obtain more measurement data. The fracture pressure, instantaneous closing pressure and reopening pressure of the rock are extracted through the pressure-time recording curve, and the horizontal maximum and minimum principal stresses and the in-situ tensile strength of the rock are calculated.

3. The method for numerical analysis of three-dimensional geostress in a geological engineering region according to claim 1, characterized in that: The step 3 specifically includes: A rose diagram of stress direction in the project area was drawn based on the basic database of China's crustal stress environment, showing the direction of the maximum principal stress in the project area and the specific stress state data range; By analyzing the basic database of China's crustal stress environment and the geological structure and stress field characteristics in the project area, the direction of the stress field in the project area is determined and more calibration points of the stress field at the model boundary are obtained.

4. The method for numerical analysis of three-dimensional geostress in a geological engineering region according to claim 1, characterized in that: The step 4 specifically includes: Obtain the regional boundary geostress data and three-dimensional geostress measured data points obtained in the previous step; The data is stored in a CSV file. Each line represents a measurement point, including longitude, latitude, and three-dimensional geostress data, which corresponds to the three geostress values. The pandas library is used to read the CSV file and convert the data into a numpy array to obtain the initial geostress data. Preliminarily check the outliers and noise points in the data set and remove the data that do not meet the standards; The data were standardized using the Min-Max or Z-Score method; Install numpy, matplotlib, scikit-learn, and pykrige libraries for Python, and use scikit-learn to train a kriging interpolation model based on Gaussian process regression: first import the Gaussian process regression and RBF variogram, then create a Gaussian process regression model that specifies the use of RBF variogram, and specify that when optimizing model parameters, restart the optimization process multiple times if a non-global optimal solution is found, and finally use the model.fit method to train the model; After the Kriging interpolation model training is completed, the trained model is used to predict the values ​​at the unobserved positions: the numpy library is used to define a set of coordinates of the unobserved positions, lon represents the longitude of the unobserved position, and lat represents the latitude of the observed position. After the coordinates of the unobserved positions are defined, the trained Kriging interpolation model is used to predict the values ​​of the unobserved positions using the predict method; Use the matplotlib library to visualize the observed data and interpolation results: first import the matplotlib.pyolot module, name it plt, and use the plt.scatter function to draw a scatter plot of the observed data coordinates. Then use the plt.contourf function to draw a contour map of the data at the unobserved locations obtained by kriging interpolation. Finally, specify the legend and title and use the plt.show function to generate and display a visualization graphic containing a scatter plot of the observed data and a contour map of the interpolation results.

5. The method for numerical analysis of three-dimensional geostress in a geological engineering region according to claim 1, characterized in that: The step 5 specifically includes: Based on the geological data obtained in the early stage and the setting of boundary conditions and material fields; The ASPECT geological modeling software was used to construct a three-dimensional geological finite element model and solve a series of equations, as follows: , （1.1） （1.2） （1.3） （1.4） Equation (1.1) represents the compressible Stokes equation, where u = u (x, t) is the velocity field and p = p (x, t) is the pressure field. Both fields depend on the spatial position x and time t. The movement of matter is driven by the gravity acting on the object and is proportional to the density and pressure of the object. η is the dynamic viscosity of the fluid, which represents the internal friction resistance of the fluid. is the strain rate tensor, which represents the deformation rate of the velocity field, is the divergence of the velocity field, describing the rate of change of the volume of the fluid, p is the pressure inside the fluid, ρ is the fluid density, and g is the gravitational acceleration vector; Equation (1.2) represents the mass conservation equation, is the divergence of mass flow, indicating the conservation of mass per unit volume, and ρu is the mass flux; Equation (1.3) represents the temperature field equation, which contains the heat conduction term and the velocity field with a flow rate of u. The right side of the equation also includes the generation of internal heat. is the specific heat capacity at constant pressure, T is the temperature, is the convection term of temperature, which indicates the transport effect of velocity field on temperature distribution. k is the heat conductivity coefficient, which describes the heat conduction capacity. is the heat conduction term, which indicates the diffusion of heat. H is the heat source term per unit volume, is the viscous dissipation term, which represents the heat generated by fluid deformation, α is the thermal expansion coefficient, which describes the volume expansion caused by temperature change, and ΔS is the entropy change rate; Equation (1.4) represents the conservation equation of chemical composition, which describes the change of chemical composition over time, including convection and diffusion. is the concentration of the i-th chemical component, is the convection term, which represents the transport effect of the velocity field on the chemical components, and qi is the chemical reaction or external source term; Ω represents the computational domain, i.e., the modeling area, and in Ω means that this set of equations is solved inside the computational domain Ω; When the model is established, the model is divided into appropriate grids to ensure a balance between the accuracy of the calculation results and the calculation efficiency. The next step is the stress field calculation stage. After the model is initialized, boundary conditions and loading conditions are applied to the model, and the calculation is run through the numerical method. The main principles are as follows: Under the action of gravity, inertial centrifugal force, Coriolis force, and pressure gradient force, the equation of motion is: （1.5） The equations of motion are: （1.6） The continuity equation is: （1.7） in, V is the velocity vector, G is the effective gravity, including the combined force of gravity and inertial centrifugal force, Ω is the angular velocity vector of the earth's rotation, r is the position vector, ρ is the fluid density, is the pressure gradient, For other external forces, ν is the kinematic viscosity coefficient, ΔV is the Laplace operator of velocity, which represents the viscous diffusion term, w , u , v For speed z , x , y The direction of the component, x , y , z is the spatial coordinate, is the component of the Coriolis parameter, g is the acceleration due to gravity, is the horizontal turbulent diffusion coefficient; Finally, the three-dimensional stress tensor distribution of each point in the model is solved.

Images