Geological engineering area three-dimensional crustal stress numerical analysis method

Abstract

The invention provides a geological engineering area three-dimensional ground stress numerical analysis method, and belongs to the technical field of underground stress analysis, and the method comprises the steps: carrying out the geological analysis of an engineering area and a surrounding area, and collecting geological data; a hydraulic fracturing stress measurement method and a hollow inclusion stress relief method are adopted for measurement at a plurality of measurement points of the engineering area; the crustal stress state of the engineering area is deeply analyzed by analyzing the Chinese crustal stress environment basic database and the geologic structure and stress field characteristics in the engineering area; obtaining a three-dimensional stress calibration point of spatial interpolation by using the regional boundary crustal stress data and the three-dimensional crustal stress measured data point and adopting a machine learning method; a three-dimensional geological finite element model is constructed by using ASPECT geological modeling software, boundary conditions and loading conditions are applied to the model, and three-dimensional stress tensor distribution of each point in the model is calculated through operation of a numerical method; and analyzing the three-dimensional stress field state of the engineering area on the basis of the crustal stress data and characteristics obtained by the steps.

Description

Technical Field

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

[0002] Hydraulic fracturing is one of the more advanced direct stress measurement technologies that has been recognized internationally in recent years. Its core advantage is that it can directly measure the current geostress in underground strata without the need for detailed knowledge of rock mechanics parameters, and can obtain multiple stress parameters. It is easy to operate, can be tested continuously or repeatedly at any depth, has a fast measurement speed, and reliable measurements. Therefore, it has been widely used in recent years and has achieved a lot of results.

[0003] However, traditional hydraulic fracturing relies on some assumptions that do not fully conform to the actual geological environment, such as the rock being linear elastic, isotropic and non-permeable, which limits the applicability of this method and can only provide limited stress parameters, and cannot fully describe the stress field characteristics in complex geological environments. These limitations make it impossible for traditional methods to effectively solve stress measurement problems under complex geological conditions. Summary of the invention

[0004] In view of the problems existing in the prior art, the present invention provides a three-dimensional geostress numerical analysis method for a geological engineering area. By using the measured point data and taking into account the influence of terrain, a three-dimensional finite element model is established. The three-dimensional stress state of the engineering area is further calculated through data assimilation, and the measured data is fused with the model data to improve the accuracy and reliability of the calculation results.

[0005] To achieve the above object, the present invention provides a three-dimensional geostress numerical analysis method for a geological engineering region, characterized in that it comprises the following steps: 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.

[0006] Optionally, 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 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.

[0007] Optionally, 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.

[0008] Optionally, 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.

[0009] Optionally, 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.

[0010] After adopting the above technical solution, the present invention has at least the following beneficial effects: The present invention significantly optimizes the application of traditional hydraulic fracturing under complex geological conditions by introducing three-dimensional finite element models and data assimilation technology, can provide more accurate and reliable underground stress measurement results, can combine multi-dimensional data, and overcome the limitations of traditional methods through detailed stress field simulation. It shows unique advantages, especially in complex environments with multi-level and multi-directional stress fields, and overcomes the disadvantage that traditional methods can only provide limited stress parameters. It not only improves 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

[0011] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0012] Figure 1 A schematic diagram of a flow chart of a method for numerically analyzing three-dimensional geostress in a geological engineering region provided by an embodiment of the present disclosure; Figure 2 A pressure-time record curve for hydraulic fracturing stress measurement; Figure 3 It is a three-dimensional visualization of the longitudinal section location map of the project area; Figure 4 This is the three-dimensional geostress distribution map of the project area. DETAILED DESCRIPTION

[0013] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0014] The present invention is described in detail with specific embodiments below: refer to Figure 1 The embodiment of the present disclosure provides a three-dimensional geostress numerical analysis method for a geological engineering region, characterized by comprising the following steps: 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; In this step, hydraulic fracturing in-situ stress measurements were carried out in two boreholes with a designed depth of 90-110 m, and stress relief measurements were carried out 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 horizontal maximum principal stress were arranged in 6 sections. The number and position of the measurement points were increased to obtain more measurement data. The pressure-time record curve was drawn based on the measurement data, such as Figure 2 As shown in the figure, the fracture pressure of rock is extracted through the pressure-time recording curve , instantaneous closing pressure and re-tensioning pressure , burst pressure It is obtained by determining the peak pressure of the first cycle. In the hydraulic fracturing test, as the fluid pressure gradually increases, when the rock reaches its tensile strength limit, cracking will occur. The peak pressure recorded at this time is the fracture pressure, and the re-tension pressure Take the point where the slope of the pressure-time curve changes significantly, usually the pressure value of the third cycle or the average of the second to fourth cycles, and the instantaneous closing pressure It is equal to the horizontal minimum principal stress. Its determination methods include the inflection point method and the single tangent method. 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 line. The single tangent method is often used to ensure accuracy and obtain the rupture pressure. , instantaneous closing pressure and re-tensioning pressure To calculate the maximum and minimum horizontal principal stresses and the in-situ tensile strength of the rock.

[0015] 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; In this step, in order to more deeply explore the characteristics and distribution patterns of the stress field in the project area, it is necessary to carefully consult and analyze the Chinese Crustal Stress Environment Basic Database, draw a rose diagram of the stress direction, and show the horizontal maximum principal stress direction and specific stress state data range of the region.

[0016] By analyzing the basic database of China's crustal stress environment and the geological structure and stress field characteristics in the project area, the aim is to determine the direction of the regional stress field and 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 provide an important basis for engineering design. The main data include the following: 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.

[0017] 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.

[0018] 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.

[0019] 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.

[0020] 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.

[0021] 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.

[0022] 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; 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.

[0023] The specific process is as follows: 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; 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; 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. 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; 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.

[0024] This step combines geographic information systems, geological engineering and machine learning technology to propose an efficient three-dimensional geostress spatial interpolation method. It uses Python's powerful data processing and machine learning library, especially the Gaussian process regression module in the scikit-learn library, to achieve site Kriging interpolation, effectively filling the spatial gap in geostress data. Through data visualization technology, the comparison between the interpolation results and the observed data is intuitively displayed, verifying the accuracy and practicality of the method. This method is not only suitable for the interpolation of geostress data, but also provides new ideas and technical paths for the analysis and processing of other geographic spatial data.

[0025] 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; In this step, the geological data obtained in the early stage, as well as the boundary conditions and material field settings are obtained. The material field component is a mathematical model used to simulate the interaction between different substances, usually used to track the movement and changes of materials. The calculation of the material field increment usually needs to be based on the reaction rate q (T, c), which involves finding the equilibrium state between different domains to achieve the occurrence of chemical reactions; 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 (unit: Pa·s), which represents the internal friction resistance of the fluid. is the strain rate tensor (unit: ), represents the deformation rate of the velocity field, is the divergence of the velocity field (unit: ), describes the rate of change of the volume of the fluid, p is the pressure (unit: Pa), the pressure inside the fluid, ρ is the fluid density (unit: kg / m³), and g is the gravitational acceleration vector (unit: m / s²), which drives the movement of the fluid; Equation (1.2) represents the mass conservation equation, is the divergence of mass flow (unit: kg / (m³•s)), which represents the conservation of mass per unit volume, ρu is the mass flux (unit: kg / (m²•s)), the product of density and velocity; 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-hand side of the equation also includes the generation of internal heat, such as heat generated by radioactive decay, friction heating, adiabatic compression of materials, phase change heat, etc. is the specific heat capacity at constant pressure (unit: J / (kg•K)), T is the temperature (unit: K), is the convection term of temperature (unit: K / s), which indicates the transport effect of velocity field on temperature distribution. k is the heat conductivity coefficient (unit: W / (m•K)), which describes the heat conduction capacity. is the heat conduction term (unit: W / m³), which indicates the diffusion of heat. H is the heat source per unit volume (unit: W / m³), is the viscous dissipation term (unit: W / m³), which represents the heat generated by fluid deformation, and α is the thermal expansion coefficient (unit: ), describes the volume expansion caused by temperature change, ΔS is the entropy change rate (unit: J / (kg•K•s)), which is related to phase change or chemical reaction; 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 ith chemical component (unit: mol / m³ or mass fraction), is the convection term (unit: mol / (m³•s)), which represents the transport effect of the velocity field on the chemical components, qi is the chemical reaction or external source term (unit: mol / (m³•s)), such as phase change, diffusion or external input; Ω 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 (unit: m / s), G is the effective gravity (unit: m / s²), including the combined force of gravity and inertial centrifugal force, Ω is the earth's rotation angular velocity vector (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 Laplace operator of velocity (unit: m / s³), represents the viscous diffusion term, w , u , v For speed z , x , y Directional component (unit: m / s), x , y , z is the spatial coordinate (unit: m), is the component of the Coriolis parameter (unit: ), g is the acceleration due to gravity (unit: m / s²), is the horizontal turbulent diffusion coefficient (unit: m² / s), which characterizes the turbulent mixing effect; Finally, the three-dimensional stress tensor distribution of each point in the model is solved.

[0026] 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.

[0027] In this step, the geological stability is analyzed based on the calculated stress field distribution, and the risks of geological disasters such as landslides and rock bursts are evaluated. The results are used to guide engineering design, such as tunnel excavation and mining.

[0028] The method described in the above embodiment is used to perform three-dimensional geostress analysis on the measured geostress data of a certain project area to obtain the three-dimensional geostress analysis result of the project area. Figure 3 The longitudinal section position of the project area is displayed, which can clearly show the spatial layout and important structural characteristics of the strata and overall ground stress, providing an intuitive reference for subsequent analysis. Figure 4 The three-dimensional geostress distribution of the project area is displayed, and the spatial distribution of geostress data is more effectively displayed, which can improve the accuracy of three-dimensional geostress distribution to the order of 10 meters, and improve the credibility and application value of the calculation results. At the same time, it also effectively eliminates the deviation between the model and reality, enhances the stress prediction ability of local areas, and provides a more scientific basis for subsequent engineering design, structural optimization, etc.

[0029] Although the present invention has been disclosed as above by way of embodiments, it is not intended to limit the present invention. Any person skilled in the art may 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 determined 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 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