A parameter optimization method for deep-shallow coupled three-dimensional geostress model based on data assimilation

Abstract

The present invention belongs to the field of earth science and engineering technology, and provides a method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation. The method is implemented by the following steps: first, collecting topographic data of the surface engineering area, analyzing seismic data to obtain underground geological structure information, and measuring three-dimensional geostress data of some observation points in the engineering area; secondly, constructing a geological model, a mechanical model and a mathematical model; then, setting loading conditions and performing finite element calculations to solve stress, displacement and deformation; then, introducing measured data for data assimilation and calibration of calculation results, by updating the initial stress parameter set of the model and performing iterative adjustments; the present invention realizes the deep fusion of observation data and numerical models through data assimilation technology, significantly improving the prediction precision and accuracy of the three-dimensional geostress model, and providing more accurate data support for engineering design and geological disaster prediction.

Description

Technical Field

[0001] The invention belongs to the field of earth science and engineering technology, and specifically relates to a method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation. Background Art

[0002] Geostress refers to the various forces that the rocks inside the earth are subjected to. Its sources are diverse, including compression at the boundaries of continental plates, changes in geothermal gradients, gravitational effects, mantle thermal convection, and magma intrusion. From the perspective of engineering, a relatively mature and accurate geostress measurement method is the hollow inclusion core stress relief method. This method is an advanced form of the casing stress relief method, which performs precise measurements by selecting appropriate original rock stress measurement points. A major advantage of this method is that it can obtain three-dimensional geostress data, accurately determine the magnitude and direction of the three principal stresses, and provide key information for engineering design. However, its application is also subject to certain limitations. Due to factors such as high economic costs and complex geological structures, the depth and number of measurement points in actual operations are often limited, which may affect the later acquisition of three-dimensional geostress data for the entire geological engineering field, and thus affect the stability of geological engineering and the further calculation of earthwork excavation.

[0003] The derivation process of the finite element method is rigorous, and it can obtain higher accuracy and effectively reduce numerical errors in the discretization process. In addition, this method has good applicability and mesh adaptability. It effectively handles problems with complex geometric shapes and free surfaces by selecting appropriate basis functions, overcomes the limitations of other methods in dealing with complex geometries and variable physical fields, and can be well used in engineering and scientific fields; and allows the use of grid units of different sizes and shapes. When dealing with problems with complex boundaries and irregular geometric shapes, the accuracy of the solution can be improved by automatically adjusting the grid and increasing the order of the elements. When conducting surface fluid-solid coupling analysis, with the help of numerical analysis theory and computer technology, the measured geostress data, engineering geological data and mathematical statistics theory are combined, and all aspects of data are fully utilized. Comprehensively consider various factors affecting geostress, so as to ensure that the stress field reproduced by the inversion is as close to the actual situation as possible, and solve complex surface simulation problems well.

[0004] However, despite the finite element method's excellent performance in dealing with fluid mechanics problems, it faces challenges in integrating complex geological data, especially the lack of effective constraints on observational data. In the traditional geological modeling process, analytical solutions to the equations of motion are difficult to obtain, and geophysical exploration methods can only cover partial areas and are costly. Therefore, data assimilation methods can effectively improve the accuracy of numerical simulation methods in numerically estimating the internal state of geological systems, allowing us to obtain detailed three-dimensional geological information.

[0005] As a key tool for understanding engineering geological conditions and making predictions, numerical models based on data assimilation play a vital role in integrating and assimilating observational data into models. This model allows us to consider dynamic processes both on the surface and deep underground, providing a comprehensive perspective for studying material movement. However, due to the uncertainty of model parameters and the fact that traditional data assimilation methods mainly focus on model state estimation rather than parameter optimization, this may lead to deviations in prediction results, thus affecting the accuracy and reliability of the model. Summary of the invention

[0006] The object of the present invention is to provide a method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation, comprising the following steps:

[0007] S1. Geological data acquisition: collect surface data, analyze seismic data to obtain underground geological structure information, and measure three-dimensional geostress data;

[0008] S2, construction of geological model, mechanical model and mathematical model;

[0009] S21. Using topographic data, underground spatial morphology and fault system characteristics, a multi-level modeling approach is used to accurately describe the regional stratigraphic and fault structures and convert them into a computable geological model;

[0010] S22, introduce rock rheological behavior, set model parameters based on rock properties and geological model data, build a temperature field model in combination with geothermal gradient, and apply external forces and boundary conditions to complete the mechanical model construction;

[0011] S23, dividing the continuous geological body into finite element units, and realizing the numerical simulation of geological phenomena by discretizing the equilibrium equations and continuous field quantities and constructing the stress-strain relationship matrix;

[0012] S3, loading and solving;

[0013] S31, Loading condition setting: Set the thermal boundary conditions according to the thermal structure of the lithosphere in the project area, set the initial temperature, set the top of the model as a free surface and the bottom as a vertical normal constraint, and apply horizontal tectonic forces in the negative directions of the X and Y axes;

[0014] S32, Finite element calculation: Perform finite element calculation on the loaded model to solve stress, displacement and deformation; use a set of equations to describe the properties of rock materials and the movement of high-density fluids;

[0015] S33, 3D geostress calibration: introduce measured data to calibrate the calculation results, update the model parameter set and make iterative adjustments, apply the adaptive parameter correction mechanism to ensure parameter adaptability, perform convergence tests and repeat iterative optimization until the accuracy requirements are met;

[0016] S4. Result analysis and geological verification: Analyze the instantaneous vertical displacement field of the surface to evaluate the terrain changes, calculate the stress intensity distribution to determine the stress concentration area, compare the results predicted by the model, and check whether the displacement field and stress field conform to the geological laws.

[0017] In a preferred solution, the specific process of step S1 includes:

[0018] S11. Collect surface data and obtain regional three-dimensional terrain information to provide high-precision initial conditions for the establishment of subsequent geostress models. The specific steps include using high-precision measurement equipment to collect data on regional terrain and spatial morphology to ensure that the data covers the entire project area; combining the collected data to generate a regional three-dimensional terrain map that meets the project accuracy requirements as the basic input for subsequent models;

[0019] S12. Analyze seismic data, extract the spatial morphology and fracture system characteristics of the stratum, and provide accurate information on fracture surface distribution for the establishment of a three-dimensional geostress model; specifically, analyze seismic wave propagation data, extract the geometric characteristics of the stratum; identify the spatial distribution, geometric morphology and properties of the fracture system in detail, especially the spatial position and extension direction of the main fracture surface;

[0020] S13. Measure three-dimensional geostress data to obtain regional geostress distribution data; use the hollow inclusion core stress relief method to measure the in-situ rock stress in the work area; ensure that the measured data covers the entire area, including the stress distribution at key fracture locations; and select measurement points in typical strata in the in-situ rock stress zone where joints and fissures are not developed or where the cementation is good.

[0021] In a preferred solution, the specific process of step S2 includes:

[0022] S21. Integrate topographic data, underground spatial morphology and fault system characteristics interpreted from seismic data; adopt a multi-level modeling approach to construct an accurate geological model by analyzing data at different levels, including surface topography, shallow strata and deep strata, and the spatial relationship between them; convert the above geological model into a computable numerical model; convert geological data into a gridded data structure, where each grid cell represents a part of the geological body, and define the relationship between nodes and grid cells so that geological phenomena can be described and calculated through mathematical equations;

[0023] S22. Introduce the rheological behavior of rocks, including viscoelasticity and plasticity, into the mechanical model, and set the parameters of the mechanical model based on the rock's density, viscosity, elastic modulus and Poisson's ratio properties and geological model data; construct a temperature field model in combination with geothermal gradient data; apply external forces and boundary conditions in the mechanical model to simulate the mechanical response in the actual geological environment, including gravity, crustal pressure, tectonic force, external forces, and top free boundary and bottom fixed boundary conditions;

[0024] S23. Divide the continuous geological body into finite element units; each unit represents a part of the geological body, and the units are interconnected through nodes; discretize the equilibrium equations describing the geological phenomena; discretize the continuous field quantities of stress and strain into the node values ​​of each unit, and perform approximate calculations of the field quantities within the unit through interpolation; construct a stress-strain relationship matrix in each unit to describe the relationship between stress and strain within the unit; and achieve numerical simulation of geological phenomena by solving the discretized set of equations.

[0025] In a preferred solution, the specific process of step S32 includes:

[0026] Discretize the regional geological model into multiple small units, namely finite elements, each of which represents a part of the geological body; define unknowns at the nodes of each unit, including velocity, pressure, and temperature, and describe the variation of these unknowns within the unit through basis functions or interpolation functions; transform continuous differential equations into discrete algebraic equations for easy numerical solution;

[0027] Solve equations, including compressible Stokes equations, temperature field equations and material composition field equations, which describe the mechanical behavior of rock materials, fluid motion and heat conduction processes. Among them, the compressible Stokes equations are:

[0028]

[0029] This equation describes the conservation of momentum in a fluid, where:

[0030] is the velocity field, which indicates the position of the fluid in space and time The speed of the next

[0031] is the pressure field, which indicates the position of the fluid in space and time Pressure under

[0032] is the viscosity of the fluid, is the density, is the acceleration due to gravity;

[0033] is the strain rate tensor, which represents the deformation rate of the fluid;

[0034] Represents the divergence operation, which is used to describe the velocity field The divergence of

[0035] Pressure field The gradient of , describing the change in pressure;

[0036] Indicates in the area The inner equation holds true;

[0037] This equation states that the movement of the fluid is driven by gravity and is proportional to the density and pressure of the fluid;

[0038]

[0039] Represents the divergence operation, used to describe mass flow The divergence of

[0040] This equation describes the conservation of mass of the fluid and ensures the continuity of the fluid;

[0041] The temperature field equation is:

[0042]

[0043] in: is the specific heat capacity, is the thermal conductivity; represents the generation of internal heat, either from radioactive decay, frictional heating, adiabatic compression, or phase change heating; is the coefficient of thermal expansion, It is entropy change;

[0044] Material composition field equations:

[0045]

[0046] This equation describes the material composition field changes, including: Indicates The concentration of a substance; Represents the source or sink term of the material component;

[0047] Pressure breakdown:

[0048] The total pressure Decomposed into static pressure and dynamic pressure ,Right now ; Static pressure is the pressure of the medium when it is stationary, i.e., hydrostatic pressure, which is calculated by the following formula:

[0049]

[0050] Assuming static pressure It is known that the momentum equation simplifies to:

[0051]

[0052] This simplified equation is used to calculate the dynamic pressure , and thus further solve the velocity field ;

[0053] By solving the above equations, we can get the velocity field , pressure field , Temperature field and material composition field The distribution of the three-dimensional ground stress field is further calculated to estimate the stress distribution state in the region.

[0054] In a preferred solution, the specific process of step S33 includes:

[0055] S331. Collect three-dimensional geostress measured data and compare the measured data with the preliminary stress distribution results obtained by finite element calculation to quantify the error between the two; calculate the covariance matrix between the initial parameter set, including density and elastic modulus material parameters, and the model estimated observation set, and analyze the uncertainty between the model and the actual observation;

[0056] S332. According to the increment of the three-dimensional geostress observation value, the covariance matrix is ​​corrected; the parameter set of the model is updated, and the updated parameter set is used as the prior parameter set of the next correction cycle; an adaptive parameter correction cycle mechanism is introduced to perform a distribution check on the updated elastic modulus parameter set; if it is found that the internal variability of the parameter set is insufficient, the expansion coefficient α is applied to enlarge and adjust the parameter prior set to generate the expanded parameters to ensure that the parameter set can effectively adapt to the new observation value; the parameter update and model solution process is repeated, and the model is optimized in combination with the new observation value in each iteration, and the parameter set is corrected again;

[0057] S333. Monitor the calculation residuals and check whether the residuals gradually decrease and tend to be stable; compare solutions with different grid accuracy to verify the grid independence of the solution; compare with the analytical solution to evaluate the accuracy of the numerical solution; guide parameter adjustment and solution strategy optimization based on the results of the convergence test; if the convergence test shows that the numerical solution has not yet met the accuracy requirements, repeat the parameter update and iterative adjustment, and convergence test steps; further adjust the model parameters, boundary conditions and initial conditions; through multiple iterations, gradually achieve the convergence and accuracy of the numerical solution until the predetermined accuracy requirements are met.

[0058] In a preferred solution, the specific process of step S4 includes:

[0059] S41. Post-process the model calculation results to extract the instantaneous vertical displacement field data of the surface; analyze the displacement field data to identify areas with significant displacement changes, which represent hot spots of terrain changes; evaluate the trend of terrain changes, including subsidence or uplift, based on the distribution of the displacement field, and predict their impact range;

[0060] S42. Calculate the stress intensity distribution within the region, including principal stress and shear stress; identify areas with higher stress intensity, which are often locations of stress concentration; analyze the location and extent of stress concentration areas and assess their impact on geological stability and engineering safety;

[0061] S43. Compare the hidden fault characteristics predicted by the model with the actual geological survey results to evaluate the prediction accuracy; compare the displacement field and stress field predicted by the model with the actual observation data to analyze the differences and consistency between the two; based on the comparison results, evaluate the prediction ability of the model and identify possible sources of error;

[0062] S44. In combination with the regional geological background, check whether the distribution of displacement field and stress field conforms to the geological structure characteristics and geological evolution laws; analyze whether there are abnormal phenomena in the model results that are inconsistent with geological laws, including unreasonable stress concentration or displacement direction; adjust the model parameters and boundary conditions according to the results of the geological rationality check.

[0063] Compared with the prior art, the present invention has the following beneficial effects:

[0064] First, the present invention designs a complete process from geological data collection, geological and mechanical model construction, loading and solving to result analysis. This process realizes the deep integration of observation data and numerical models through data assimilation technology, significantly improving the prediction accuracy and accuracy of the three-dimensional geostress model. With the continuous iterative adjustment of model parameters, the model can better fit the actual geological conditions, greatly reduce model deviations, enhance the applicability and reliability of the model in complex geological environments, and provide more accurate data support for key areas such as engineering design and geological disaster prediction, which helps to reduce decision-making risks and improve the success rate of project implementation.

[0065] Second, the present invention designs specific steps for geological data collection, model construction, loading solution and three-dimensional geostress calibration. By introducing an adaptive parameter correction mechanism, it ensures that the parameter set maintains sufficient internal variability during the iterative optimization process, effectively avoiding the problem of model rigidity caused by rapid parameter convergence. The enhanced parameter optimization capability enables the model to respond quickly and adjust its own parameters when faced with new observation data, thereby maintaining a high degree of flexibility and adaptability. This not only improves the practicality and accuracy of the model, but also provides a more robust solution for simulation prediction under complex geological conditions.

[0066] Third, this method provides key geological information support for engineering design by accurately simulating the three-dimensional geostress distribution, which helps to optimize the project layout and construction plan, reduce the engineering risks caused by the uncertainty of geological conditions, and improve the safety and stability of the project. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 is a flow chart of the method of the present invention;

[0068] Figure 2 is a schematic diagram of constructing a three-dimensional model mesh in step S2;

[0069] Figure 3 is the initial geothermal cloud map of the calculation area constructed by the mechanical model in step S2;

[0070] Figure 4 It is a longitudinal section cloud map of a water conservancy project (nearly north-south direction) in a certain area in Example 1;

[0071] Figure 5 is the regional in-situ stress map calculated before data assimilation in Example 1;

[0072] Figure 6 It is the regional in-situ stress map calculated after data assimilation in Example 1. DETAILED DESCRIPTION

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

[0074] This paper proposes an innovative data assimilation scheme, which aims to correct the three-dimensional geostress by assimilating geological observation data and thus improve the accuracy of numerical simulation. The implementation of this scheme can be divided into the following four carefully designed steps:

[0075] S1. Geological data collection and analysis;

[0076] Geological data collection and analysis is the basic link of the whole process, and its main purpose is to provide high-precision input data for the establishment of subsequent geostress models. The present invention completes data collection and analysis through the following three key steps:

[0077] Surface data measurement of the work area: In order to obtain the three-dimensional terrain information of the area, high-precision measurement equipment is used to collect the surface data of the work area. The present invention completes the collection and processing of surface data through the following process:

[0078] 1) Use high-precision measurement equipment to collect data on regional terrain and spatial morphology to ensure that the data covers the entire project area;

[0079] 2) Combine the collected data to generate a regional three-dimensional topographic map that meets the engineering accuracy requirements, including longitude and latitude and corresponding elevation information as the basic input for subsequent models.

[0080] This step can provide high-resolution surface topography data, which provides key initial conditions for stress distribution in the area and its simulation.

[0081] Seismic data interpretation: Seismic data provides an important basis for analyzing underground geological structures. The spatial morphology and fracture system characteristics of the strata are accurately extracted through the following processes:

[0082] 1) Analyze seismic wave propagation data and extract the geometric characteristics of the internal strata;

[0083] 2) Detailed identification of the spatial distribution, geometric shape and properties of the fault system, especially the spatial location and extension direction of the main fault surface.

[0084] The accurate interpretation of seismic data provides accurate information on the distribution of fault surfaces, including the dip, inclination, extension direction and length of the fault surface, for the establishment of a three-dimensional geostress model. It also lays the foundation for subsequent calibration and verification of the model.

[0085] Three-dimensional geostress measurement Geostress data collection: When obtaining regional geostress distribution data, the hollow inclusion core stress relief method is used to measure the original rock stress in the work area, and ensure that the measured data covers the entire area, especially the stress distribution at the key fracture position. To ensure the representativeness and accuracy of the data, the selection of measurement points must follow the following principles:

[0086] Representative stratigraphic areas: The measuring points should be selected in typical stratigraphic areas to ensure that the measuring results can reflect the overall distribution characteristics of regional geostress;

[0087] Complete or relatively complete rock mass: The measuring point should be located in the area where joints and fissures are not developed or the cementation is good, and avoid choosing the broken zone or the area with developed fractures;

[0088] Stay away from engineering disturbance areas: The measuring points should avoid areas that are greatly disturbed by engineering, such as tunnels, goafs, and large caves, to ensure that the measuring points are located in the original rock stress area;

[0089] Reasonable arrangement of measuring points: According to geological conditions and engineering requirements, the location and number of measuring points should be reasonably arranged to ensure that the measured stress can accurately reflect the distribution law of ground stress in the entire work area from a spatial perspective.

[0090] The above process can ensure the high accuracy and representativeness of the three-dimensional geostress measured data, and provide strong data support for the calibration and optimization of subsequent models.

[0091] S2. Model establishment;

[0092] Geological model construction:

[0093] The present invention provides a method for constructing a geological model based on multi-level modeling, which is used to accurately describe regional strata and fault structures. A regional geological model is established by comprehensively utilizing terrain data, underground space morphology, and geometric characteristics of the fault system. The geological data is converted into a computable numerical model to further improve the applicability and reliability of the calculation process. In terms of grid division, this method uses two grid forms: Euler grid and Lagrangian grid (such as Figure 2 As shown in Figure 1). The Euler grid is used to simulate the motion characteristics of observation points in a fixed space, while the Lagrangian grid dynamically adjusts the grid nodes as the particles move to achieve the evolution description of the motion variation. The introduction of this dual-grid system not only ensures the flexibility of calculation, but also significantly improves the accuracy and stability of numerical simulation, providing a reliable basic support for the subsequent analysis and optimization of the geostress field.

[0094] Mechanical model construction:

[0095] In order to construct a mechanical model that can accurately reflect the characteristics of geological structures, it is first necessary to introduce the rheological behavior of rocks and materials (such as viscoelasticity and plasticity). This is because these characteristics determine the deformation and failure characteristics of rocks under external conditions such as stress and temperature. The setting of model parameters is mainly based on the main rock properties and mechanical characteristics of different strata. The parameter settings of the mechanical model can be further optimized, including density, viscosity, elastic modulus, internal friction angle, cohesion and Poisson's ratio to ensure the accuracy of the model. Next, based on the geothermal gradient data provided by the geological model, a temperature field distribution model is constructed, and the boundary conditions are adjusted in combination with the temperature change law of deep strata. This is because temperature has a significant influence on the mechanical properties of rocks, and its role runs through the geological model and the mechanical model. It is an important factor in improving the scientific nature of the model (such as Figure 3 On this basis, by enabling the adaptive mesh function of the finite element, the resolution of key areas (such as the basin center and the main fault zone) is further improved, so that the model can accurately reflect local details while maintaining computational efficiency. Finally, combined with the key parameters selected in the previous steps (such as gravity, formation density 2800kg / m³, internal friction angle 30°, cohesion , elastic modulus 3.5 (GPa), Poisson's ratio 0.25), applying horizontal tectonic force and top free boundary conditions to complete the construction of the mechanical model.

[0096] Mathematical model construction:

[0097] The continuous geological body is divided into finite element units, each unit represents a part of the geological body and shares node information, and the continuous field (such as stress, strain) is discretized into the node value of each unit. After discretizing the region, the equilibrium equation is discretized by converting the differential equation into algebraic equations on the nodes, and then the stress-strain relationship matrix is ​​constructed in each unit. The constitutive relationship is discretized.

[0098] S3. Loading solution;

[0099] The loading solution part is the core step to realize the simulation of geostress distribution, and the geological stress distribution is analyzed through finite element calculation.

[0100] (1) Loading condition setting:

[0101] In order to simulate the influence of regional thermal structure and tectonic stress on geological structure, it is necessary to set thermal boundary conditions in the dynamic model in combination with the thermal structure of the lithosphere in the engineering area. Then, initial conditions are applied, among which the key is the setting of initial temperature and initial geostress, and the temperature field distribution is constructed through geothermal gradient data, in order to ensure that the model can accurately reflect the influence of regional geothermal environment on mechanical properties.

[0102] By loading boundary conditions, the model can simulate the mechanical response of the structural unit under the combined action of external force, initial geostress and temperature field, and further analyze the displacement difference, stress distribution and contact behavior between adjacent structural units. In order to set reasonable loading conditions, including fixed thermal boundaries (such as the fixation of the bottom of the rock mass) and stress loading boundaries (such as regional crustal pressure). First of all, it is necessary to combine the geological background of the engineering area. Specifically, the top boundary of the model is set as a free surface boundary to allow it to deform strongly perpendicular to the depth direction; a vertical normal constraint is applied to the bottom of the model to simulate the fixation of deep strata. At the same time, horizontal tectonic forces are applied in the negative direction of the X and Y axes to simulate the influence of actual geological tectonic forces. The setting of this series of loading conditions enables the model to more realistically reflect the coupling between external forces, temperature and stress in the geological environment, and provide a scientific basis for the analysis of the mechanical properties of regional geological structures.

[0103] (2) Finite element calculation:

[0104] The loaded model is calculated using the finite element analysis method: the regional geological model is discretized into small units through meshing, and the stress, displacement and deformation of each unit are solved. A series of equations are solved. The basic equations follow the properties of rock materials, and the equations also describe the movement of high-density fluids. This step is actually to estimate the three-dimensional ground stress state. The specific process is:

[0105] Discretize the regional geological model into multiple small units, namely finite elements, each of which represents a part of the geological body; define unknowns at the nodes of each unit, including velocity, pressure, and temperature, and describe the variation of these unknowns within the unit through basis functions or interpolation functions; transform continuous differential equations into discrete algebraic equations for easy numerical solution;

[0106] Solve equations, including compressible Stokes equations, temperature field equations and material composition field equations, which describe the mechanical behavior of rock materials, fluid motion and heat conduction processes. Among them, the compressible Stokes equations are:

[0107]

[0108] This equation describes the conservation of momentum in a fluid, where:

[0109] is the velocity field, which indicates the position of the fluid in space and time The speed of the next

[0110] is the pressure field, which indicates the position of the fluid in space and time Pressure under

[0111] is the viscosity of the fluid, is the density, is the acceleration due to gravity;

[0112] is the strain rate tensor, which represents the deformation rate of the fluid;

[0113] Represents the divergence operation, which is used to describe the velocity field The divergence of

[0114] Pressure field The gradient of , describing the change in pressure;

[0115] Indicates in the area The inner equation holds true;

[0116] This equation states that the movement of the fluid is driven by gravity and is proportional to the density and pressure of the fluid;

[0117]

[0118] Represents the divergence operation, used to describe mass flow The divergence of

[0119] This equation describes the conservation of mass of the fluid and ensures the continuity of the fluid;

[0120] The temperature field equation is:

[0121]

[0122] in: is the specific heat capacity, is the thermal conductivity; represents the generation of internal heat, either from radioactive decay, frictional heating, adiabatic compression, or phase change heating; is the coefficient of thermal expansion, It is entropy change;

[0123] Material composition field equations:

[0124]

[0125] This equation describes the material composition field changes, including: Indicates The concentration of a substance; Represents the source or sink term of the material component;

[0126] Pressure breakdown:

[0127] The total pressure Decomposed into static pressure and dynamic pressure ,Right now ; Static pressure is the pressure of the medium when it is stationary, i.e., hydrostatic pressure, which is calculated by the following formula:

[0128]

[0129] Assuming static pressure It is known that the momentum equation simplifies to:

[0130]

[0131] This simplified equation is used to calculate the dynamic pressure , and thus further solve the velocity field ;

[0132] By solving the above equations, we can get the velocity field , pressure field , Temperature field and material composition field The distribution of the three-dimensional ground stress field is further calculated to estimate the stress distribution state in the region.

[0133] (3) Three-dimensional geostress calibration:

[0134] Introduce measured data and calibrate the calculated results;

[0135] After completing the preliminary estimation of the geostress state, the measured data is introduced into the model calibration process. By comparing the errors of the finite element calculation results with the measured stress distribution, the deviation of the model results is quantified, and the covariance matrix between the initial geostress parameters and the model estimated observation set is calculated. The calculation of covariance is the basis for uncertainty analysis between the model and the actual observations, and can be used to guide parameter updates. The specific process includes:

[0136] Collect three-dimensional geostress measured data and compare them with the preliminary stress distribution results obtained by finite element calculation to quantify the error between the two; calculate the initial geostress parameters and the covariance matrix between the model estimated observation set, and analyze the uncertainty between the model and the actual observations;

[0137] According to the increment of the three-dimensional geostress observation value, the covariance matrix is ​​corrected; the parameter set of the model is updated, and the updated parameter set is used as the prior parameter set for the next correction cycle; an adaptive parameter correction cycle mechanism is introduced to perform a distribution check on the updated geostress parameter set; if the internal variability of the parameter set is found to be insufficient, the expansion coefficient α is applied to amplify and adjust the parameter prior set to generate the expanded parameters to ensure that the parameter set can effectively adapt to the new observation value; the parameter update and model solution process is repeated, and the model is optimized in combination with the new observation value in each iteration, and the parameter set is corrected again;

[0138] Monitor the calculation residuals to check whether the residuals gradually decrease and tend to stabilize; compare solutions with different grid accuracies to verify the grid independence of the solution; compare with the analytical solution to evaluate the accuracy of the numerical solution; guide parameter adjustment and solution strategy optimization based on the results of the convergence test; if the convergence test shows that the numerical solution has not yet met the accuracy requirements, repeat the parameter update and iterative adjustment, and convergence test steps; further adjust the model parameters, boundary conditions, and initial conditions; through multiple iterations, gradually achieve the convergence and accuracy of the numerical solution until the predetermined accuracy requirements are met.

[0139] S4. Parameter update and iterative adjustment;

[0140] First, according to the increment of the three-dimensional geostress observations, the covariance matrix is ​​corrected, and the parameter set of the model is updated accordingly. The updated parameter set is used as the prior parameter set for the next correction cycle to further optimize the model performance. On this basis, an adaptive parameter correction cycle mechanism is introduced to perform a distribution check on the updated geostress parameter set. If it is found that the internal variability of the parameter set is insufficient, it may cause the parameters to lose their diffusion quickly, thereby affecting the flexibility and adaptability of the model. To address this problem, the expansion coefficient α is applied to amplify and adjust the parameter prior set to generate expanded parameters to ensure that the parameter set can effectively adapt to the new observations. In the parameter perturbation stage, we performed parameter perturbations on the measured geostress data. The specific formula is:

[0141] F=F0+α*δF0

[0142] Among them, F0(i,j) is the two-dimensional distribution value in the estimated stress grid characteristic map, and (i,j) corresponds to the longitude and latitude coordinates of the project area. deltaF0=0.1*F0, deltaF0 represents the ground stress disturbance value, which is 10% of the initial value. α is the expansion coefficient, and the value range is -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9. 20 parallel simulation experiments were carried out for each expansion coefficient to obtain the corresponding simulation data.

[0143] Specifically, F represents the change in geostress (i.e., the value in the geostress grid characteristic map, which is the independent variable), and Y corresponds to the geostress value of each point in the forward model calculation (i.e., the dependent variable). Through statistical analysis of the relationship between the independent variable and the dependent variable, we can calculate the Pearson correlation coefficient to measure the linear correlation between the two:

[0144]

[0145] After calculating the slope K, compare the difference ΔY between the in-situ stress value of the simulation result when the expansion coefficient α=0 and the gridded simulation data (the calculation result when the disturbance coefficient α≠0). According to the linear relationship Y=K*X+b, ΔX is solved, that is, the in-situ stress change that needs to be adjusted, and this change is added to the initial in-situ stress parameter.

[0146] Subsequently, the parameter update and model solving process are repeated. To ensure the stability of model convergence, we assimilate 1-2 observations into the model each time, gradually optimize the model with new observations in each iteration, and calibrate the parameter set again. Through multiple iterations, until the parameter correction reaches the convergence standard, all observations are effectively used.

[0147] Convergence test:

[0148] To ensure the accuracy and stability of the numerical solution, the numerical solution of the model is tested for convergence. Convergence means that as the iteration proceeds, the solution gradually stabilizes and the residual gradually decreases until the set convergence criteria are met. Usually, convergence is judged based on the change in the residual, especially the norm of the residual. If the norm is less than a predetermined threshold (such as 1e-6), it indicates convergence. At the same time, convergence can also be judged by monitoring the relative change of the solution or the stability of the physical quantity. For example, if the solution change for several consecutive iterations is less than a certain percentage (such as 0.01%), it can be considered that the solution has converged. The results of the convergence test are directly used to guide parameter adjustment and solution strategy optimization to improve the reliability of the model.

[0149] Iterative solution and accuracy optimization:

[0150] If the convergence test shows that the numerical solution does not meet the accuracy requirements, repeat the above steps to further adjust the model parameters, boundary conditions and initial conditions. Through multiple iterative optimizations, the convergence and accuracy of the numerical solution are gradually achieved.

[0151] Embodiment 1:

[0152] This embodiment selects a water conservancy engineering geological structure in a certain area with complex geological structures, including multiple layers of strata and fault systems. The strata in the area include shallow strata and deep strata, and different strata have different rock properties and mechanical properties, such as density, viscosity, elastic modulus, internal friction angle, cohesion and Poisson's ratio. There are multiple fracture surfaces in the area, and the spatial position and extension direction of these fracture surfaces have an important influence on the distribution of ground stress. By analyzing seismic data, the geometric characteristics and spatial distribution of the fracture system are extracted. This embodiment uses measured ground stress data, which are obtained by the hollow inclusion core stress relief method to ensure the accuracy and representativeness of the data. The measured data covers the entire engineering area, including the stress distribution at key fracture locations, providing a reliable basis for the construction and calibration of the model.

[0153] This embodiment accurately simulates and predicts the distribution of geostress in the region, including the assessment of terrain changes, stress concentration areas, and potential geological disaster risks. High-precision measurement equipment is used to collect regional three-dimensional terrain information and generate a regional three-dimensional topographic map that meets the engineering accuracy requirements. The seismic wave propagation data is analyzed to extract the geometric characteristics of the stratum and the spatial distribution of the fracture system. The hollow inclusion core stress relief method is used to measure the original rock stress in the work area to ensure that the measured data covers the entire area, including the stress distribution at key fracture locations.

[0154] This embodiment uses a multi-level modeling method to integrate terrain data, underground space morphology and fracture system characteristics, build an accurate geological model, and convert it into a computable numerical model. The rheological behavior of rocks is introduced, model parameters are set based on rock properties and geological model data, a temperature field model is constructed in combination with geothermal gradients, and external forces and boundary conditions are applied.

[0155] The continuous geological body is divided into finite element units, and the stress-strain relationship matrix is ​​constructed by discretizing the equilibrium equation and continuous field quantity to realize the numerical simulation of geological phenomena. The thermal boundary conditions are set according to the thermal structure of the lithosphere in the project area, and the temperature field distribution is constructed. The initial temperature and initial geostress are set to ensure that the model can accurately reflect the influence of the regional geothermal environment on the mechanical properties.

[0156] The top of the model is set as a free surface boundary, a vertical normal constraint is applied to the bottom, and a horizontal structural force is applied in the negative direction of the X and Y axes. The loaded model is subjected to finite element calculation to solve stress, displacement and deformation. By solving the compressible Stokes equation, temperature field equation and material composition field equation, the distribution of velocity field, pressure field, temperature field and material composition field is obtained, and the three-dimensional geostress field is further calculated. The measured data is introduced to calibrate the calculation results, and the adaptability of the model parameters is ensured by updating the model parameter set and iteratively adjusting it, and convergence test and repeated iterative optimization are performed until the predetermined accuracy requirements are met.

[0157] According to the increment of the three-dimensional ground stress observation value, the covariance matrix is ​​corrected, the parameter set of the model is updated, and an adaptive parameter correction loop mechanism is introduced to ensure that the parameter set can effectively adapt to the new observation value. Through multiple iterations, the convergence and accuracy of the numerical solution are gradually achieved. First, the model calculation results are post-processed to analyze the instantaneous vertical displacement field distribution of the regional surface to evaluate the terrain change trend. Through this analysis, potential subsidence areas or uplift areas can be identified, providing a scientific basis for geological disaster warning and regional planning. Subsequently, the stress intensity distribution in the region is calculated, the location and range of the stress concentration area are analyzed, and the potential damage area or landslide hazard area is located accordingly, providing a key reference for project site selection and disaster prevention and control. Finally, the geological rationality of the model results is tested, including comparing the hidden fault characteristics predicted by the model with the actual fault distribution, evaluating the prediction accuracy, and combining the regional geological background to check whether the distribution of the displacement field and stress field conforms to the geological law.

[0158] According to the above method, combined with the engineering geological measured data and finite element simulation results in this embodiment, the geostress data assimilation was performed. The geostress results before and after data assimilation are as follows:

[0159] Table 1: Comparison of geostress after data assimilation in Example 1

[0160] Note: In the table, compressive stress is positive and tensile stress is negative.

[0161] Table 1 shows the comparison between the measured stress values ​​of stress components (σx and σy) at different measuring points (ZK208, ZK210, ZKL01, ZKL02), the "observed values" obtained by finite element simulation, and the in-situ stress values ​​after data assimilation of the present invention. It can be seen from Table 1 that at all measuring points, there is a certain difference between the "observed values" of finite element simulation in the prior art and the measured stress. For example, at the ZK208 measuring point, the measured stress of σx is 9.06MPa, while the "observed value" of finite element simulation is 10.5MPa; the measured stress of σy is 6.68MPa, while the "observed value" of finite element simulation is 7.2MPa. This difference shows that the traditional finite element simulation method may not fully and accurately reflect the actual in-situ stress distribution in some cases.

[0162] Through the data assimilation method of the present invention, the obtained ground stress value is closer to the measured stress value. For example, at the ZK208 measuring point, the ground stress value of σx after data assimilation is 10.1MPa, which is closer to the measured value of 9.06MPa; the ground stress value of σy after data assimilation is 6.85MPa, which is also closer to the measured value of 6.68MPa. This shows that the data assimilation method of the present invention can significantly improve the accuracy of ground stress simulation and make the simulation results closer to the actual observation data.

[0163] The data in Table 1 show that the data assimilation method of the present invention can effectively reduce the difference between the simulated values ​​and the measured values ​​at multiple measuring points, thereby improving the prediction accuracy of the geostress model.

[0164] For example, at the ZKL01 measuring point, the measured stress of σx is 8.57MPa, the "observed value" of the finite element simulation is 8MPa, and the ground stress value after data assimilation is 8.71MPa, which is closer to the measured value; the measured stress of σy is 5.49MPa, the "observed value" of the finite element simulation is 2.8MPa, and the ground stress value after data assimilation is 4.89MPa, which is also closer to the measured value.

[0165] In summary, the data in Table 1 show that the data assimilation method of the present invention can significantly improve the accuracy of geostress simulation and make the simulation results closer to the actual observation data, thereby providing more accurate data support for engineering design and geological disaster prediction.

[0166] Reference Figure 4-Figure 6 Through the data assimilation method of the present invention, the measured data and the numerical model are deeply integrated, which significantly improves the accuracy of the ground stress distribution. The data assimilation method ensures that the model can better fit the actual geological conditions and reduce the model deviation by updating the model parameter set and making iterative adjustments.

[0167] The application of the data assimilation method of the present invention significantly improves the prediction accuracy and precision of the geostress model. Through multiple iterations of optimization, the model can gradually achieve the convergence and accuracy of the numerical solution until the predetermined accuracy requirement is met. This enhances the applicability and reliability of the model in complex geological environments.

[0168] By introducing measured data for calibration, the model can better reflect the actual situation. Figure 5 (before data assimilation) and Figure 6 (After data assimilation) It can be seen that after data assimilation, the calculated regional geostress in the lower terrain area shows a more realistic geostress distribution. This shows that the data assimilation method can effectively eliminate the deviation between the model and reality and improve the credibility of the model.

[0169] The model after data assimilation shows stronger ability in stress prediction. The model can more accurately predict stress changes and provide a more scientific basis for subsequent engineering design, structural optimization, etc. For example, when evaluating terrain changes, stress concentration areas, and potential geological disaster risks, the model after data assimilation can provide more reliable prediction results.

[0170] This embodiment not only reflects a more accurate distribution of geostress through measured data and data assimilation methods, improves the credibility and application value of the calculation results, but also effectively eliminates the deviation between the model and reality, and enhances the stress prediction ability of the region. This provides a more scientific and reliable basis for subsequent engineering design, structural optimization, and geological disaster prediction. The method of the present invention not only significantly improves the efficiency of parameter correction, but also enhances the adaptability and responsiveness of the model to observed data, thereby providing an efficient and stable solution for three-dimensional geostress simulation under complex geological conditions.

[0171] It is obvious to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the present invention can be implemented in other specific forms without departing from the spirit or essential features of the present invention. Therefore, the embodiments of the present invention are exemplary and non-restrictive.

Claims

1. A method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation, characterized in that: The steps include: S1. Geological data acquisition: collecting surface data, analyzing seismic data to obtain underground structural information, and measuring three-dimensional geostress data; S2. Model construction, including: S21. Use multi-level modeling to describe the regional stratigraphic and fracture structures and convert them into computational geological models; S22, introduce rock rheological behavior, set mechanical model parameters, construct temperature field model and apply boundary conditions; S23, divide the finite element unit, discretize the equilibrium equation and field quantity, construct the stress-strain matrix, and realize numerical simulation; S3, loading solution, S31, setting thermal boundary conditions, initial temperature and boundary constraints, and applying horizontal structural forces; S32, perform finite element calculation to solve stress, displacement and deformation; S33, calibrating the result with the measured data, and iteratively adjusting the parameters until the accuracy requirements are met; S4. Result analysis and verification: analyze terrain changes and stress concentration areas, compare model prediction results, and verify the rationality of displacement and stress fields.

2. The method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation according to claim 1 is characterized by: The specific process of step S1 includes: S11. Collect surface data and obtain regional three-dimensional terrain information to provide high-precision initial conditions for the establishment of subsequent geostress models. The specific steps include using high-precision measurement equipment to collect data on regional terrain and spatial morphology to ensure that the data covers the entire project area; combining the collected data to generate a regional three-dimensional terrain map that meets the project accuracy requirements, including longitude and latitude and corresponding elevation information as the basic input for subsequent models; S12. Analyze seismic data, extract the spatial morphology and fracture system characteristics of the stratum, and provide accurate information on fracture surface distribution for the establishment of a three-dimensional geostress model; specifically, analyze seismic wave propagation data, extract the geometric characteristics of the stratum; identify the spatial distribution, geometric morphology and properties of the fracture system in detail, especially the spatial position and extension direction of the main fracture surface; S13. Measure three-dimensional geostress data to obtain regional geostress distribution data; use the hollow inclusion core stress relief method to measure the in-situ rock stress in the work area; ensure that the measured data covers the entire area, including the stress distribution at key fracture locations; and select measurement points in typical strata in the in-situ rock stress zone where joints and fissures are not developed or where the cementation is good.

3. The method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation according to claim 1 is characterized in that: The specific process of step S2 includes: S21. Integrate topographic data, underground spatial morphology and fault system characteristics interpreted from seismic data; adopt a multi-level modeling approach to construct an accurate geological model by analyzing data at different levels, including surface topography, shallow strata and deep strata, and the spatial relationship between them; convert the above geological model into a computable numerical model; convert geological data into a gridded data structure, where each grid cell represents a part of the geological body, and define the relationship between nodes and grid cells so that geological phenomena can be described and calculated through mathematical equations; S22. Introduce the rheological behavior of rocks into the mechanical model, including viscoelasticity and plasticity, and set the parameters of the mechanical model based on the rock's density, viscosity, elastic modulus, internal friction angle, cohesion and Poisson's ratio properties and geological model data; construct a temperature field model in combination with geothermal gradient data; apply external forces and boundary conditions in the mechanical model to simulate the mechanical response in the actual geological environment, including gravity, crustal pressure, tectonic force, external forces, and top free boundary and bottom fixed boundary conditions; S23. Divide the continuous geological body into finite element units; each unit represents a part of the geological body, and the units are interconnected through nodes; discretize the equilibrium equations describing the geological phenomena; discretize the continuous field quantities of stress and strain into the node values ​​of each unit, and perform approximate calculations of the field quantities within the unit through interpolation; construct a stress-strain relationship matrix in each unit to describe the relationship between stress and strain within the unit; and achieve numerical simulation of geological phenomena by solving the discretized set of equations.

4. The method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation according to claim 1 is characterized by: The specific process of step S32 includes: Discretize the regional geological model into multiple small units, namely finite elements, each of which represents a part of the geological body; define unknowns at the nodes of each unit, including velocity, pressure, and temperature, and describe the variation of these unknowns within the unit through basis functions or interpolation functions; transform continuous differential equations into discrete algebraic equations for easy numerical solution; Solve equations, including compressible Stokes equations, temperature field equations and material composition field equations, which describe the mechanical behavior of rock materials, fluid motion and heat conduction processes. Among them, the compressible Stokes equations are: This equation describes the conservation of momentum in a fluid, where: is the velocity field, which indicates the position of the fluid in space and time The speed of the next is the pressure field, which indicates the position of the fluid in space and time Pressure under 1 represents the unit tensor, which is used to represent the isotropic part in tensor operations; is the viscosity of the fluid, is the density, is the acceleration due to gravity; is the strain rate tensor, which represents the deformation rate of the fluid; Represents the divergence operation, which is used to describe the velocity field The divergence of Pressure field The gradient of , describing the change in pressure; Indicates in the area The inner equation holds true; This equation states that the movement of the fluid is driven by gravity and is proportional to the density and pressure of the fluid; Represents the divergence operation, used to describe mass flow The divergence of This equation describes the conservation of mass of the fluid and ensures the continuity of the fluid; The temperature field equation is: in: is the specific heat capacity, is the thermal conductivity; represents the generation of internal heat, either from radioactive decay, frictional heating, adiabatic compression, or phase change heating; is the coefficient of thermal expansion, It is entropy change; Material composition field equations: This equation describes the material composition field changes, including: Indicates The concentration of a substance; Represents the source or sink term of the material component; Pressure breakdown: The total pressure Decomposed into static pressure and dynamic pressure ,Right now ; Static pressure is the pressure of the medium when it is stationary, i.e., hydrostatic pressure, which is calculated by the following formula: Assuming static pressure It is known that the momentum equation simplifies to: Calculate the dynamic pressure using formula 1.5 , and thus further solve the velocity field ; By solving the above equations, we can get the velocity field , pressure field , Temperature field and material composition field The distribution of the three-dimensional ground stress field is further calculated to estimate the stress distribution state in the region.

5. The method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation according to claim 1 is characterized by: The specific process of step S33 includes: S331. Collect three-dimensional geostress measured data and compare the measured data with the preliminary stress distribution results obtained by finite element calculation to quantify the error between the two; calculate the covariance matrix between the initial parameter set, including density and elastic modulus material parameters, and the model estimated observation set, and analyze the uncertainty between the model and the actual observation; S332. According to the increment of the three-dimensional geostress observation value, the covariance matrix is ​​corrected; the parameter set of the model is updated, and the updated parameter set is used as the prior parameter set of the next correction cycle; an adaptive parameter correction cycle mechanism is introduced to perform a distribution check on the updated elastic modulus parameter set; if it is found that the internal variability of the parameter set is insufficient, the expansion coefficient α is applied to enlarge and adjust the parameter prior set to generate the expanded parameters to ensure that the parameter set can effectively adapt to the new observation value; the parameter update and model solution process is repeated, and the model is optimized in combination with the new observation value in each iteration, and the parameter set is corrected again; S333. Monitor the calculation residuals and check whether the residuals gradually decrease and tend to be stable; compare solutions with different grid accuracy to verify the grid independence of the solution; compare with the analytical solution to evaluate the accuracy of the numerical solution; guide parameter adjustment and solution strategy optimization based on the results of the convergence test; if the convergence test shows that the numerical solution has not yet met the accuracy requirements, repeat the parameter update and iterative adjustment, and convergence test steps; further adjust the model parameters, boundary conditions and initial conditions; through multiple iterations, gradually achieve the convergence and accuracy of the numerical solution until the predetermined accuracy requirements are met.

6. The method for optimizing parameters of a deep-shallow coupled three-dimensional geostress model based on data assimilation according to claim 1 is characterized by: The specific process of step S4 includes: S41. Post-process the model calculation results to extract the instantaneous vertical displacement field data of the surface; analyze the displacement field data to identify areas with significant displacement changes, which represent hot spots of terrain changes; evaluate the trend of terrain changes, including subsidence or uplift, based on the distribution of the displacement field, and predict their impact range; S42. Calculate the stress intensity distribution within the area, including principal stress and shear stress; analyze the location and range of stress concentration areas and evaluate their impact on geological stability and engineering safety; S43. Compare the hidden fault characteristics predicted by the model with the actual geological survey results to evaluate the prediction accuracy; compare the displacement field and stress field predicted by the model with the actual observation data to analyze the differences and consistency between the two; based on the comparison results, evaluate the prediction ability of the model and identify possible sources of error; S44. In combination with the regional geological background, check whether the distribution of displacement field and stress field conforms to the geological structure characteristics and geological evolution laws; analyze whether there are abnormal phenomena in the model results that are inconsistent with geological laws, including unreasonable stress concentration or displacement direction; adjust the model parameters and boundary conditions according to the results of the geological rationality check.

Images