Terrestrial heat and earthquake symbiosis mechanism analysis method and system, electronic equipment and storage medium

Abstract

The invention belongs to the technical field of geophysical exploration of geothermal resources and seismic activities, and provides a geothermal and seismic symbiosis mechanism analysis method and system, electronic equipment and a storage medium, and the method comprises the steps of three-dimensional joint inversion, computational domain grid division, attribute addition, boundary condition addition, reference temperature addition, simulation calculation and correlation analysis. According to the method, three-dimensional joint inversion is performed on the electromagnetic data and the gravity data, so that the model better conforms to actual underground geologic structure characteristics, idealized constraints of a theoretical model are separated, and the accuracy of a simulation result is improved; through attribute addition, boundary condition addition, reference temperature addition, simulation calculation and correlation analysis, a complete earthquake activity and geothermal resource deep control and symbiosis theory system is constructed, and extraction of internal correlation of geothermal activity and earthquake occurrence is realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the geophysical exploration technical field of geothermal resources and seismic activity, and particularly relates to a geothermal and seismic symbiotic mechanism analysis method and system, an electronic device and a storage medium. BACKGROUND

[0002] Global tectonic activity is the result of relative motion between plates. The compression and collision between plates make the interface area of the modern lithospheric plate the most active area, which is characterized by abnormal display of high heat flow area, volcanic activity, magmatic intrusion, orogenic movement, seismic activity and metamorphism, etc., which provides sufficient heat source and upward channel for the formation and distribution of high-temperature geothermal fluid. The explosion of deep fluid upward migration and escape produces earthquakes of different magnitudes and focal depths in the earth's interior. Fluid in the crust plays a crucial role in the release of seismic and thermal energy. Faults, as the migration channel of fluid and groundwater or the important storage space of thermal material, are the main control elements of deep geothermal resources and seismic activity. The open faults under active seismic activity allow groundwater to circulate to the deep high-temperature area and flow back to the surface in the form of geothermal fluid to form geothermal resources. The volatile substances produced in the fluid migration process enter the thermal reservoir or the surrounding rock layer, causing the increase of pore pressure and the rupture of rock layer, thereby triggering earthquakes. Therefore, geothermal resources and seismic disasters are closely related in terms of formation mechanism. The characterization method of geothermal and seismic symbiotic mechanism is mainly geochemical analysis and geophysical exploration. Geochemical methods can discuss the relationship between geothermal and seismic by analyzing the source and genesis of geothermal fluid, thermal reservoir temperature, gas genesis, and mantle fluid release intensity. Geophysical exploration obtains the underground geological structure by inversion, providing spatial information for the discussion of the relationship between geothermal resources and seismic activity. However, the external cause of seismic activity and geothermal resources is the result of tectonic regime, and the internal cause depends on the material and structural properties. The current research method only stays in the observation characteristics and conversion relationship of the two, and lacks a complete deep control and symbiotic system of seismic activity and geothermal resources and related theories. Establishing a unified framework is a shortcut to solve the problems of resources and disasters, and is also a key scientific problem for the strategic goals of a livable earth and safe supply of new energy.

[0003] Traditional heat-flow-solid coupling models have double limitations in application scenarios and research scales. On the one hand, their applications are mostly focused on the field of enhanced geothermal production, such as exploring the influence of different fracture types on heat transfer by simulating the actual production injection process, or focusing on the microscopic simulation of rock skeleton and pore at the laboratory scale to explain the coupling mechanism of multiple physical fields, which has not broken through the scale bottleneck and expanded to the macroscopic simulation at the crust level, so it is difficult to analyze the interaction between actual geothermal resources and earthquake disasters. On the other hand, due to the level of technology development, the current research on the relationship between geothermal and earthquake still has significant defects: most researches do not fully consider the fluid thermal convection effect and excessively rely on idealized theoretical models, resulting in a large deviation between research results and real geological processes, and the internal relationship between geothermal and earthquake phenomena cannot be accurately described. SUMMARY

[0004] In order to overcome the shortcomings of the prior art, the purpose of the present application is to provide a geothermal and earthquake symbiotic mechanism analysis method, system, electronic equipment and storage medium, which solves the problems that the traditional model is difficult to analyze the interaction between actual geothermal resources and earthquake disasters and excessively relies on idealized theoretical models.

[0005] To achieve the above purpose, the present application provides the following scheme:

[0006] A geothermal and earthquake symbiotic mechanism analysis method, comprising:

[0007] Three-dimensional joint inversion is performed on the pre-acquired electromagnetic data and gravity data to obtain a three-dimensional resistivity and density model;

[0008] The three-dimensional resistivity and density model is divided into a plurality of calculation domain grids;

[0009] Fluid field, heat transfer field, deformation field, heat source attribute, thermal physical parameter and mechanical parameter are added to each calculation domain grid after the division to obtain an attribute grid;

[0010] Boundary conditions are added to the fluid field, the heat transfer field and the deformation field of each attribute grid to obtain a simulation grid;

[0011] Reference temperature is added to the simulation grid according to the distribution position of the simulation grid to obtain a heat-flow-solid coupling model;

[0012] The heat-flow-solid coupling model is simulated under different heat source conditions to obtain a temperature field in the crust and a strain field in the crust;

[0013] The coupling relationship between the geothermal and the earthquake is obtained by using the GIS spatial analysis technology to superimpose and compare the spatial distribution data of the temperature field in the crust, the strain field in the crust, and the spatial distribution data of the surface geothermal resources and the underground earthquake disasters, and to analyze the spatial coupling relationship.

[0014] Preferably, the pre-acquired electromagnetic data and gravity data are subjected to three-dimensional joint inversion to obtain a three-dimensional resistivity and density model, including:

[0015] The electromagnetic data are subjected to three-dimensional MT inversion to obtain an MT inversion result; the objective function of the three-dimensional MT inversion is: ; wherein, ; is the objective function value of the three-dimensional MT inversion; is a data fitting functional; represents a data type, is 1 for MT, is 2 for gravity; is a model parameter vector; is a regularization factor; is a model stability functional; is a constraint term;

[0016] The MT inversion result is converted into a structure reference model according to a physical property prior information conversion method; an expression of the structure reference model includes: , ; wherein, is a fitting functional for the electromagnetic data; is a regularization factor for the electromagnetic data; is a model parameter vector for the electromagnetic data; is an objective function value of the structure reference model conversion; is a parameter vector for converting the electromagnetic model into a gravity model; is a prior rock physical property correlation function of density and resistivity;

[0017] The structure reference model is added to a preset wide-range constraint, and three-dimensional gravity joint inversion is performed on the gravity data under the wide-range constraint to obtain a gravity inversion result; a calculation formula of the three-dimensional gravity inversion is: ; wherein, is the gravity inversion result; is a fitting functional for the gravity data; is a model parameter vector for the gravity data; is a regularization factor for the gravity data; is a physical property constraint term of Gramian constraint;

[0018] iterating the MT inversion result and the gravity inversion result to obtain the three-dimensional resistivity and density model.

[0019] Preferably, the three-dimensional resistivity and density model is divided into a plurality of calculation domain grids, including:

[0020] According to the collected geological information, the three-dimensional resistivity and density model is divided into different calculation domain blocks by using an interpolation function;

[0021] Defining the geological structure of each calculation domain block; the geological structure includes: thickness, geometric shape, lithology;

[0022] Using the finite element method to perform tetrahedral meshing on each calculation domain block to obtain a plurality of calculation domain grids.

[0023] Preferably, the thermal physical parameters include: density, thermal conductivity, heat generation rate, constant pressure heat capacity, thermal expansion coefficient, specific heat rate; the mechanical parameters include: Young's modulus, Poisson's ratio, porosity, permeability.

[0024] Preferably, the boundary condition addition is performed on the fluid field, the heat transfer field, and the deformation field of each attribute grid to obtain a simulation grid, including:

[0025] The inlet of the fluid field is set to a water head of 15km, the outlet of the fluid field is set to a water head of 6km, the upper boundary temperature of the aquifer of the fluid field is set to 623.15K, and the lower boundary temperature of the aquifer of the fluid field is set to 787.15K;

[0026] The upper boundary annual average ground surface temperature of the heat transfer field is set to 293.15K, and the lower boundary condition of the heat transfer field is set to a mantle heat flow of 40mW / m 2 , the side boundary condition of the heat transfer field is set to adiabatic;

[0027] The lower boundary of the deformation field is set to a fixed boundary, the upper boundary of the deformation field is set to a free surface, and the side boundary of the deformation field is set to a symmetric boundary to obtain the simulation grid.

[0028] Preferably, the reference temperature addition is performed on the simulation grid according to the distribution position of the simulation grid to obtain a thermal-fluid-solid coupling model, including:

[0029] The upper crust ground surface temperature of the simulation grid is set to 293.15K;

[0030] The lower crust rock mass part temperature of the simulation grid is set to 773.15K;

[0031] The temperature of the magma part of the simulation grid is set to 873.15K, and the thermo-hydro-mechanical coupling model is obtained.

[0032] Preferably, the thermo-hydro-mechanical coupling model is simulated under different heat source conditions to obtain the temperature field in the crust and the strain field in the crust, including:

[0033] A heat conduction equation is constructed; the heat conduction equation is: ; wherein, is the effective thermal conductivity; is the temperature field variable; is the fluid density; is the fluid constant-pressure specific heat capacity; is the seepage velocity; is the volumetric heat source; is the viscous dissipation heat;

[0034] A fluid seepage equation control equation is constructed; the fluid seepage equation control equation includes: , ; wherein, is the mass source term; is the permeability tensor; is the fluid viscosity; is the pore fluid pressure; is the gravitational acceleration;

[0035] A thermal stress control equation is constructed; the thermal stress control equation includes: , , , ; wherein, is the total stress tensor; is the solid position vector; is the thermal strain; is the linear thermal expansion coefficient; is the reference temperature; is the unit tensor; is the stress tensor; is the stiffness matrix; is the body force;

[0036] The thermo-hydro-mechanical coupling model is simulated by using the heat conduction equation, the fluid seepage equation control equation, and the thermal stress control equation, and the temperature field in the crust and the strain field in the crust are obtained.

[0037] Preferably, a geothermal and earthquake coexistence mechanism analysis system includes:

[0038] A joint inversion model is used to perform three-dimensional joint inversion on pre-acquired electromagnetic data and gravity data to obtain a three-dimensional resistivity and density model;

[0039] a mesh division module, configured to divide the three-dimensional resistivity and density model into a plurality of calculation domain meshes;

[0040] an attribute adding module, configured to add fluid field, heat transfer field, deformation field, heat source attribute, thermal physical parameter and mechanical parameter to each of the calculation domain meshes after the division, to obtain attribute meshes;

[0041] a boundary condition adding module, configured to add boundary conditions to the fluid field, the heat transfer field and the deformation field of each of the attribute meshes, to obtain simulation meshes;

[0042] a reference temperature adding module, configured to add reference temperature to the simulation meshes according to the distribution positions of the simulation meshes, to obtain a thermo-hydro-mechanical coupling model;

[0043] a simulation calculation module, configured to simulate the thermo-hydro-mechanical coupling model under different heat source conditions, to obtain an internal temperature field and an internal strain field of the crust;

[0044] a correlation analysis module, configured to superimpose and compare the internal temperature field and the internal strain field of the crust with the spatial distribution data of surface geothermal resources and underground seismic disasters by using GIS spatial analysis technology, and analyze the spatial coupling relationship, to obtain a geothermal and seismic coupling relationship analysis result.

[0045] Preferably, an electronic device comprises at least one processor and a memory connected with the processor in communication; the memory stores instructions executable by the processor, and the instructions are executed by the processor to enable the processor to execute the foregoing geothermal and seismic symbiotic mechanism analysis method.

[0046] Preferably, a non-transitory computer-readable storage medium stores computer instructions for enabling a computer to execute the foregoing geothermal and seismic symbiotic mechanism analysis method.

[0047] The present application discloses the following technical effects:

[0048] The present application provides a geothermal and seismic symbiotic mechanism analysis method, system, electronic device and storage medium, which solves the problem of excessive dependence on idealized theoretical models of traditional models by performing three-dimensional joint inversion on electromagnetic data and gravity data, and improves the accuracy of simulation results; through attribute addition, boundary condition addition, reference temperature addition, simulation calculation and correlation analysis, the problem that traditional models are difficult to analyze the interaction relationship between actual geothermal resources and seismic disasters is solved, and the extraction of the internal correlation between geothermal activity and earthquake occurrence is realized. BRIEF DESCRIPTION OF DRAWINGS

[0049] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the accompanying drawings needed in the embodiments will be briefly introduced as follows. Obviously, the accompanying drawings in the following description only represent some embodiments of the present application, and all other drawings obtained by those of ordinary skill in the art without creative effort based on these drawings belong to the protection scope of the present application.

[0050] Figure 1 A geothermal and earthquake symbiotic mechanism analysis process schematic diagram is provided for the embodiments of the present application.

[0051] Figure 2 A geological model diagram is provided for the embodiments of the present application, (a) is a three-dimensional geological model schematic diagram, and (b) is an N-S direction profile diagram selected from the three-dimensional model.

[0052] Figure 3 Numerical simulation calculation result schematic diagrams are provided for the embodiments of the present application, (a) is a temperature simulation result diagram, (b) is an equivalent stress simulation result diagram, (c) is a shear stress simulation result diagram, and (d) is a distance distribution schematic diagram.

[0053] Figure 4 Temperature field and stress field result profile diagrams are provided for the embodiments of the present application, (a) is a temperature field N-S direction result profile diagram, (b) is an equivalent stress N-S direction result profile diagram, and (c) is a shear stress N-S direction result profile diagram.

[0054] Figure 5 Temperature curve and well measured temperature curve comparison diagrams are provided for the embodiments of the present application, (a) is a C-2 well measured and simulated temperature curve comparison diagram, and (b) is a C-3 well measured and simulated temperature curve comparison diagram. DETAILED DESCRIPTION

[0055] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments only represent some embodiments of the present application, and all other embodiments obtained by those of ordinary skill in the art without creative effort based on these drawings belong to the protection scope of the present application.

[0056] The present application aims to provide a geothermal and earthquake symbiotic mechanism analysis method, system, electronic device and storage medium, to solve the problems that the traditional model is difficult to analyze the interaction relationship between the actual geothermal resources and the earthquake disasters and excessively relies on idealized theoretical models.

[0057] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application will be further described in detail below with reference to the accompanying drawings and specific embodiments. In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0058] Figure 1 A geothermal and earthquake symbiotic mechanism analysis process schematic diagram provided for an embodiment of the present application is shown in FIG. 1, and the present application provides a geothermal and earthquake symbiotic mechanism analysis method, including: Figure 1

[0059] Step 100: performing three-dimensional joint inversion on pre-acquired electromagnetic data and gravity data to obtain a three-dimensional resistivity and density model;

[0060] Step 200: dividing the three-dimensional resistivity and density model into a plurality of calculation domain grids;

[0061] Step 300: adding fluid field, heat transfer field, deformation field, heat source attribute, thermal physical property parameter, and mechanical parameter to each of the calculation domain grids after the division to obtain an attribute grid;

[0062] Step 400: adding boundary condition to the fluid field, the heat transfer field, and the deformation field of each of the attribute grids to obtain a simulation grid;

[0063] Step 500: adding reference temperature to the simulation grid according to the distribution position of the simulation grid to obtain a thermal-fluid-solid coupling model;

[0064] Step 600: simulating and calculating the thermal-fluid-solid coupling model under different heat source conditions to obtain an intracrustal temperature field and an intracrustal strain field;

[0065] Step 700: using GIS spatial analysis technology to superimpose and compare the intracrustal temperature field, the intracrustal strain field, and the spatial distribution data of surface geothermal resources and underground earthquake disasters, and to analyze the spatial coupling relationship to obtain a geothermal and earthquake coupling relationship analysis result.

[0066] Specifically, the three-dimensional joint inversion on the pre-acquired electromagnetic data and gravity data to obtain the three-dimensional resistivity and density model includes:

[0067] performing three-dimensional MT inversion on the electromagnetic data to obtain an MT inversion result; the objective function of the three-dimensional MT inversion is: ; wherein, ; is the objective function value of the three-dimensional MT inversion; is a data fitting functional; represents the data type, is 1 for MT, is 2 for gravity; is a model parameter vector; is a regularization factor; is a model stability functional; is a constraint term;​

[0068] transforming the MT inversion result into a structure reference model according to a physical property prior information transformation method; an expression of the structure reference model comprises: 、 ; wherein, is a fitting functional for the electromagnetic data; is a regularization factor for the electromagnetic data; is a model parameter vector for the electromagnetic data; is a target function value of the structure reference model conversion; is a parameter vector of electromagnetic model conversion into a gravity model; is a prior rock physical property correlation function of density and resistivity;

[0069] adding the structure reference model into a preset wide range constraint, and performing three-dimensional gravity joint inversion on the gravity data under the wide range constraint to obtain a gravity inversion result; a calculation formula of the three-dimensional gravity inversion is: ; wherein, is the gravity inversion result; is a fitting functional for the gravity data; is a model parameter vector for the gravity data; is a regularization factor for the gravity data; is a physical property constraint term of Gramian constraint;

[0070] iterating the MT inversion result and the gravity inversion result to obtain the three-dimensional resistivity and density model.

[0071] Further, the three-dimensional resistivity and density model is divided into a plurality of calculation domain grids, comprising:

[0072] according to collected geological information, the three-dimensional resistivity and density model is divided into different calculation domain blocks by using an interpolation function;

[0073] defining a geological structure of each of the calculation domain blocks; the geological structure comprises: thickness, geometric shape, lithology;

[0074] using a finite element method to perform tetrahedral mesh partitioning on each of the calculation domain blocks to obtain a plurality of the calculation domain grids.

[0075] Specifically, the thermal physical property parameters comprise: density, thermal conductivity, heat generation rate, constant pressure heat capacity, thermal expansion coefficient, specific heat rate; the mechanical parameters comprise: Young's modulus, Poisson's ratio, porosity, permeability.

[0076] Further, boundary condition adding is performed on the fluid field, the heat transfer field and the deformation field of each of the attribute grids to obtain a simulation grid, including:

[0077] The inlet of the fluid field is set to a water head of 15 km, the outlet of the fluid field is set to a water head of 6 km, the upper boundary temperature of the aquifer of the fluid field is set to 623.15 K, and the lower boundary temperature of the aquifer of the fluid field is set to 787.15 K;

[0078] The upper boundary annual average ground surface temperature of the heat transfer field is set to 293.15 K, and the lower boundary condition of the heat transfer field is set to a mantle heat flow of 40 mW / m 2 The side boundary condition of the heat transfer field is set to adiabatic;

[0079] The lower boundary of the deformation field is set to a fixed boundary, the upper boundary of the deformation field is set to a free surface, and the side boundary of the deformation field is set to a symmetric boundary to obtain the simulation grid.

[0080] Specifically, reference temperature adding is performed on the simulation grid according to the distribution position of the simulation grid to obtain a heat-flow-solid coupling model, including:

[0081] The upper crust ground surface temperature of the simulation grid is set to 293.15 K;

[0082] The lower crust rock mass part temperature of the simulation grid is set to 773.15 K;

[0083] The magma part temperature of the simulation grid is set to 873.15 K to obtain the heat-flow-solid coupling model.

[0084] Further, simulation calculation is performed on the heat-flow-solid coupling model under different heat source conditions to obtain a temperature field in the crust and a strain field in the crust, including:

[0085] A heat conduction equation is constructed; the heat conduction equation is: wherein, is an effective thermal conductivity; is a temperature field variable; is a fluid density; is a fluid constant-pressure specific heat capacity; is a seepage velocity; is a volumetric heat source; is viscous dissipation heat;

[0086] A fluid seepage equation control equation is constructed; the fluid seepage equation control equation includes: , wherein, is a mass source term; is a permeability tensor; is a fluid viscosity; is a pore fluid pressure; is a gravitational acceleration;

[0087] constructing a thermal stress governing equation; the thermal stress governing equation comprises: , , , ; wherein, is a total stress tensor; is a solid displacement vector; is a thermal strain; is a linear thermal expansion coefficient; is a reference temperature; is a unit tensor; is a stress tensor; is a stiffness matrix; is a body force;

[0088] simulating the thermal-fluid-solid coupling model by using the thermal conduction equation, the fluid seepage equation governing equation and the thermal stress governing equation, to obtain a temperature field and a strain field in the crust.

[0089] Specifically, a geothermal and earthquake coexistence mechanism analysis method based on thermal-fluid-solid multi-physical field coupling numerical simulation, a fine crust resistivity and density structure is obtained by three-dimensional joint inversion of electromagnetic (MT) and gravity data, and a numerical simulation geometric model is constructed. According to the geological information, the calculation domain is divided and meshed, the material properties, thermal properties and mechanical parameters are added, the heat transfer, fluid and deformation field boundary conditions are set, and a thermal-fluid-solid coupling model is formed. After simulating the temperature field and stress field under different heat sources, the geothermal and earthquake activity rules and causes are revealed, and the coexistence deep driving mechanism of geothermal resources and earthquake activity is explored. The specific implementation method comprises:

[0090] S1: According to the resistivity model, the actual underground formation structure (undulation change, etc.) is obtained, and then combined with the density model and surface geological information, etc. Different lithology is supplemented to the above structure, and finally the underground geological model that fits the actual situation is obtained. Three-dimensional joint inversion of electromagnetic data and gravity data is used to obtain a three-dimensional resistivity and density model of crustal scale, which is highly consistent with the actual geological conditions. The model is used as the geometric model of numerical simulation;

[0091] S2: According to the collected geological information, the geometric model is divided into different calculation domains by using an interpolation function, and the tetrahedral mesh is divided in the calculation domain by using a finite element method;

[0092] S3: Add the corresponding material properties to each calculation domain, i.e., divide into porous medium heat transfer field, fluid field and solid deformation field, and different heat source properties (such as magma heat source and radioactive heat source), and then add the thermal physical parameters and mechanical parameters of the corresponding rock layer;

[0093] S4: Add appropriate boundary conditions to the fluid field, heat transfer field and deformation field according to the actual situation;

[0094] S5: Establish a coupling relationship to form a heat-flow-solid coupling model;

[0095] S6: Simulate to obtain the temperature field and strain field in the crust generated under different heat source conditions;

[0096] S7: Combine the spatial distribution of surface geothermal resources and underground seismic disasters to obtain the law and cause of geothermal and seismic activity, and finally obtain the correlation between the two in the deep control mechanism.

[0097] Further, in S1, three-dimensional wide-range property constraint joint inversion is performed using electromagnetic data and gravity data to obtain a crust-scale resistivity and density model as an initial geometric model for numerical simulation. The core of the inversion method is to use the correlation of rock physical parameters such as resistivity and density to constrain joint inversion, realize the complementarity of different geophysical methods, and the specific inversion strategy is as follows: performing three-dimensional MT inversion on the electromagnetic data to obtain an MT inversion result; the objective function of the three-dimensional MT inversion is: ; wherein, is the MT inversion objective function value; is a data fitting functional; i represents the data type (MT or gravity); is a model parameter vector; is a regularization factor; is a model stability functional; is a constraint term, which is represented by Gramian constraint, and is represented as: ; wherein, is a Gramian constraint coupling term; is a penalty function constraint term; is a weight coefficient of the Gramian constraint term, which is adaptively adjusted using a decay strategy; is a weight coefficient of the penalty function term;

[0098] The MT inversion result is converted into a structure reference model through physical property prior information to guide the subsequent joint inversion of gravity data, and the calculation formula is: ; wherein, is a model parameter vector of electromagnetic data; is a prior rock physical property correlation function of density and resistivity;

[0099] The structure reference model is added to the wide range constraint, and then three-dimensional gravity joint inversion is performed to obtain a gravity inversion result; the calculation formula of the three-dimensional gravity inversion is: ; wherein, is the gravity inversion result; is a fitting functional of gravity data; is a model parameter vector of gravity data; is a regularization factor of gravity data; is a physical constraint term of Gramian constraint;

[0100] The MT inversion result is converted and the gravity joint inversion is iterated to obtain the three-dimensional resistivity and density model.

[0101] Specifically, in S2, a subsequent numerical simulation model is established using a comsol platform, geological information of a research area is collected, a rough structure in a geometric model is subdivided into undulating layers or irregular blocks by using an interpolation function, an actual underground geological structure is simulated, thickness, geometric shape and lithology are defined, and a three-dimensional geological model as shown in Figure 2 (a) and an N-S profile geological structure as shown in (b) are obtained. Figure 2 (b) are obtained.

[0102] Further, in S3, thermal physical parameters of the rock include: density, thermal conductivity, heat generation rate, constant pressure heat capacity, thermal expansion coefficient, specific heat rate; mechanical parameters include: Young's modulus, Poisson's ratio, porosity, permeability.

[0103] Specifically, in S4, in order to simulate an actual underground geological condition, appropriate boundary conditions are set for different calculation domains, an inlet of a fluid field is set to a water head of 15 km, an outlet is set to a water head of 6 km, an upper boundary temperature of an aquifer is 623.15 K, and a lower boundary temperature is 787.15 K; an upper boundary of a porous medium heat transfer field is set to an annual average ground surface temperature of the research area, 293.15 K, a lower boundary condition is a mantle heat flow of 40 mW / m 2 (a fixed base heat flow), and a side boundary condition of the model is adiabatic; a lower boundary of a deformation field is fixed, an upper boundary is a free surface, and a side boundary is a symmetric boundary.

[0104] Further, in S5, in the process of geothermal generation, unevenly distributed lithospheric temperature and thermal stress generated by thermal structure will accumulate stress in the brittle layer of the middle and upper crust until an earthquake occurs, so a thermal-flow-solid multi-physical field coupling is realized by using this relationship, and in the case of considering thermal expansion, different reference temperatures are set, the upper crust is set to a ground surface temperature of 293.15 K, the lower crust rock part is set to 773.15 K, and the magma part is set to 873.15 K.

[0105] Specifically, in S6, the actual underground situation is simulated to obtain the temperature distribution and thermal strain generated by different heat sources in the steady state. The specific steps are as follows: the temperature distribution is calculated by considering conduction and convection, and the thermal strain is calculated by using the thermal expansion coefficient The thermal strain is introduced into the elastic mechanics equation as the “initial strain”, and the stress tensor and equivalent stress are solved to obtain the stress tensor and equivalent stress. The specific formula includes:

[0106] 1) Heat conduction equation (porous medium, steady state):

[0107]

[0108] This formula describes the heat transport process of saturated fluid in a porous medium, where, is the effective thermal conductivity (considering solid and fluid), is the fluid density, is the specific heat capacity of the fluid at constant pressure, is the seepage velocity, is the temperature field variable, is the volumetric heat source (considering magma heat source and radioactive heat source), is the viscous dissipation heat, which can be ignored in most actual underground water and geothermal modeling, and in the steady state:

[0109]

[0110] 2) Fluid seepage equation control equation:

[0111] Mass conservation equation of porous medium flow:

[0112]

[0113] This formula describes the conservation of local fluid mass, that is, the fluid entering the unit volume is balanced with the fluid source term, where is the fluid density, is the seepage velocity, is the mass source term (such as fluid generation and consumption);

[0114] Darcy's law:

[0115]

[0116] This formula describes how the fluid flows in the porous medium due to the pressure gradient and gravity, where, is the permeability tensor, is the fluid viscosity, is the pore fluid pressure, is the fluid density, is the gravitational acceleration, and the model assumes that the fluid is incompressible, isothermal, and laminar.

[0117] 3) Thermal stress control equation (solid mechanics):

[0118] The control equation is a small deformation linear elastic mechanics equation, considering thermal expansion strain, stress-strain relationship:

[0119]

[0120]

[0121]

[0122] Force balance equation:

[0123]

[0124] where, is the solid position vector, is the body force, mainly gravity.

[0125] Further, Figure 3 (a) to Figure 3 (c) are the temperature field, equivalent stress field, and shear stress field simulated using the above model, respectively. It can be seen that the temperature field has a large change near the brittle-ductile transition zone at about 10 km and the moho surface, generating a large equivalent stress. Earthquakes also occur at the corresponding positions, such as Figure 4 (a) to Figure 4 (b) shown. At the same time, the equivalent stress generated in the northern part of the study area at a depth of 10 km is higher than that in the southern part, and the southern part generates a large equivalent stress at the depth of the moho surface. Reference Figure 3 (c) and Figure 3 (d) The results of shear stress show that at a depth of 10 km, earthquakes are mostly in the compressive zone, and hot springs and geothermal wells (star-shaped and red rectangular) are at the junction of the compressive zone and the tensile zone, and the shape of the tensile zone is consistent with the strike of the Yalu River fault. The profile also shows that earthquakes in the northern part of the study area are all in the compressive zone (reference Figure 4 ). From the temperature curve, it can be clearly seen (reference Figure 5 ), that C-2 (reference Figure 5 (a)) and C-3 wells (reference Figure 5 ) are in the tensile zone.The simulated temperature curve, the measured temperature curve and the calculated temperature curve of the prior art are highly overlapped, proving that the simulation result has strong reliability and credibility. The simulation result proves that the heat source generates heat, which makes the temperature distribution of the lithosphere uneven, affecting the stress state in the crust. Compared with the radioactive heat source of the rock layer, the thermal stress generated by the magma heat source is higher, and the influence on the thermal state is greater. The existence of the aquifer enhances the compressive stress above the brittle-ductile transition zone, which helps to accumulate stress until an earthquake occurs. The enhancement of the tensile stress around the hot spring and the geothermal well helps the heat to rise to the surface to produce thermal anomalies. The above analysis reveals the dynamic coupling mechanism between the heat source, the stress and the earthquake, characterizes the symbiotic relationship between the geothermal and the earthquake, and provides a theoretical basis for the development of geothermal resources and the assessment of seismic risk.

[0126] As an optional implementation, the embodiment further provides a geothermal and earthquake symbiotic mechanism analysis system, comprising:

[0127] A joint inversion model is used to perform three-dimensional joint inversion on the pre-acquired electromagnetic data and gravity data to obtain a three-dimensional resistivity and density model;

[0128] A grid division module is used to divide the three-dimensional resistivity and density model into a plurality of calculation domain grids;

[0129] An attribute adding module is used to add fluid field, heat transfer field, deformation field, heat source attribute, thermal physical parameter and mechanical parameter to each calculation domain grid after the division to obtain an attribute grid;

[0130] A boundary condition adding module is used to add boundary conditions to the fluid field, the heat transfer field and the deformation field of each attribute grid to obtain a simulation grid;

[0131] A reference temperature adding module is used to add reference temperature to the simulation grid according to the distribution position of the simulation grid to obtain a thermal-fluid-solid coupling model;

[0132] A simulation calculation module is used to perform simulation calculation on the thermal-fluid-solid coupling model under different heat source conditions to obtain a temperature field in the crust and a strain field in the crust;

[0133] A correlation analysis module is used to perform superposition comparison and spatial coupling relationship analysis on the temperature field in the crust, the strain field in the crust and the spatial distribution data of the surface geothermal resources and the underground earthquake disasters by using GIS spatial analysis technology to obtain a geothermal and earthquake coupling relationship analysis result.

[0134] As an optional implementation, the embodiment further provides an electronic device, comprising at least one processor and a memory connected with the processor in communication; wherein the memory stores instructions executable by the processor, and the instructions are executed by the processor to enable the processor to execute the geothermal and earthquake symbiotic mechanism analysis method.

[0135] As an optional implementation, the embodiment further provides a non-transitory computer readable storage medium storing computer instructions, the computer instructions being used to enable a computer to execute the geothermal and earthquake symbiotic mechanism analysis method.

[0136] The beneficial effects of the present application are as follows:

[0137] The present application makes the model more consistent with the actual underground geological structure characteristics by three-dimensional joint inversion of electromagnetic data and gravity data, breaks away from the idealized constraints of the theoretical model, and improves the accuracy of the simulation results; through attribute addition, boundary condition addition, reference temperature addition, simulation calculation and correlation analysis, a complete theory system of deep control and symbiosis of seismic activity and geothermal resources is constructed, and the internal correlation between geothermal activity and earthquake occurrence is extracted.

[0138] The embodiments in the specification are described in a progressive manner, and each embodiment focuses on the difference from other embodiments, and the same or similar parts between the embodiments can be referred to each other.

[0139] The principles and implementation modes of the present application are described by applying specific examples herein, and the above description of the embodiments is only used to help understand the method of the present application and its core idea; at the same time, for those skilled in the art, according to the idea of the present application, the specific implementation mode and application range will be changed. In conclusion, the content of the specification should not be understood as a limitation of the present application.

Claims

1. A method for analyzing the mechanism of geothermal and earthquake symbiosis, characterized in that: include: Perform 3D joint inversion of pre-collected electromagnetic and gravity data to obtain 3D resistivity and density models; dividing the three-dimensional resistivity and density model into a plurality of computational domain grids; Adding fluid field, heat transfer field, deformation field, heat source properties, thermophysical parameters, and mechanical parameters to each of the divided computational domain grids to obtain an attribute grid; Adding boundary conditions to the fluid field, the heat transfer field, and the deformation field of each of the attribute grids to obtain a simulation grid; Adding a reference temperature to the simulation grid according to the distribution position of the simulation grid to obtain a thermal-fluid-solid coupling model; Simulating the heat-fluid-solid coupling model under different heat source conditions to obtain the temperature field and strain field in the crust; GIS spatial analysis technology is used to perform superposition comparison and spatial coupling relationship analysis on the temperature field in the crust, the strain field in the crust, and the spatial distribution data of surface geothermal resources and underground earthquake disasters to obtain the analysis results of the geothermal and earthquake coupling relationship.

2. The method for analyzing the geothermal and earthquake symbiosis mechanism according to claim 1, characterized in that: Perform a 3D joint inversion of pre-acquired electromagnetic and gravity data to obtain a 3D resistivity and density model, including: Perform three-dimensional MT inversion on the electromagnetic data to obtain an MT inversion result; the objective function of the three-dimensional MT inversion is: ;in, ; is the objective function value of the three-dimensional MT inversion; Fitting functionals to the data; Indicates the data type, 1 means MT, 2 indicates gravity; is the model parameter vector; is the regularization factor; is the model stabilizing functional; is a constraint item; The MT inversion result is converted into a structural reference model according to a physical property prior information conversion method; the expression of the structural reference model includes: 、 ;in, is the fitting functional for the electromagnetic data; is a regularization factor for the electromagnetic data; is a model parameter vector for the electromagnetic data; The objective function value of the structural reference model conversion; is the parameter vector converted from electromagnetic model to gravity model; is the prior rock property correlation function of density and resistivity; The structural reference model is added to a preset wide range constraint, and a three-dimensional gravity joint inversion is performed on the gravity data under the wide range constraint to obtain a gravity inversion result; the calculation formula for the three-dimensional gravity inversion is: ;in, is the gravity inversion result; is the fitting functional for the gravity data; is a model parameter vector for the gravity data; is the regularization factor for the gravity data; is the physical property constraint of Gramian constraint; The MT inversion result and the gravity inversion result are iterated to obtain the three-dimensional resistivity and density model.

3. The method for analyzing the geothermal and earthquake symbiosis mechanism according to claim 2, characterized in that: The three-dimensional resistivity and density model is divided into several computational domain grids, including: Based on the collected geological information, the three-dimensional resistivity and density model is divided into different computational domain blocks using an interpolation function; Defining the geological structure of each computational domain block; the geological structure includes: thickness, geometry, and lithology; Each computational domain block is tetrahedron meshed using the finite element method to obtain a plurality of computational domain meshes.

4. The method for analyzing the geothermal and earthquake symbiosis mechanism according to claim 3, characterized in that: The thermophysical parameters include: density, thermal conductivity, heat generation rate, constant pressure heat capacity, thermal expansion coefficient, and specific heat rate; the mechanical parameters include: Young's modulus, Poisson's ratio, porosity, and permeability.

5. The method for analyzing the geothermal and earthquake symbiosis mechanism according to claim 4, characterized in that: Adding boundary conditions to the fluid field, the heat transfer field, and the deformation field of each of the attribute grids to obtain a simulation grid includes: The inlet of the fluid field is set to a water head of 15 km, the outlet of the fluid field is set to a water head of 6 km, the upper boundary temperature of the aquifer of the fluid field is set to 623.15 K, and the lower boundary temperature of the aquifer of the fluid field is set to 787.15 K; The upper boundary of the heat transfer field is set to an annual average surface temperature of 293.15 K, and the lower boundary condition of the heat transfer field is set to a mantle heat flux of 40 mW / m 2 , setting the side boundary condition of the heat transfer field to adiabatic; The lower boundary of the deformation field is set as a fixed boundary, the upper boundary of the deformation field is set as a free surface, and the side boundaries of the deformation field are set as symmetric boundaries to obtain the simulation grid.

6. The method for analyzing the geothermal and earthquake symbiosis mechanism according to claim 5, characterized in that: Adding a reference temperature to the simulation grid according to the distribution position of the simulation grid to obtain a thermal-fluid-solid coupling model, including: The upper crust surface temperature of the simulation grid is set to 293.15K; The temperature of the lower crust rock mass portion of the simulation grid is set to 773.15K; The temperature of the magma portion of the simulation grid is set to 873.15 K to obtain the thermal-fluid-solid coupling model.

7. The method for analyzing the geothermal and earthquake symbiosis mechanism according to claim 6, characterized in that: The heat-fluid-solid coupling model is simulated and calculated under different heat source conditions to obtain the temperature field and strain field in the crust, including: Construct a heat conduction equation; the heat conduction equation is: ;in, is the effective thermal conductivity; is the temperature field variable; is the fluid density; is the specific heat capacity of the fluid at constant pressure; is the seepage velocity; is the volume heat source; Heat dissipated by viscosity; Construct a fluid seepage equation control equation; the fluid seepage equation control equation includes: 、 ;in, is the mass source term; is the permeability tensor; is the fluid viscosity; is the pore fluid pressure; is the acceleration due to gravity; Construct a thermal stress control equation; the thermal stress control equation includes: 、 、 、 ;in, is the total stress tensor; is the solid potential vector; is thermal strain; is the linear thermal expansion coefficient; is the reference temperature; is the unit tensor; is the stress tensor; is the stiffness matrix; is the volume force; The heat conduction equation, the fluid seepage control equation, and the thermal stress control equation are used to simulate and calculate the heat-fluid-solid coupling model to obtain the temperature field and the strain field in the crust.

8. A geothermal and earthquake symbiotic mechanism analysis system, characterized in that: include: Joint inversion model, used to perform 3D joint inversion of pre-acquired electromagnetic and gravity data to obtain 3D resistivity and density models; A grid partitioning module, for partitioning the three-dimensional resistivity and density model into a plurality of computational domain grids; An attribute adding module is used to add fluid field, heat transfer field, deformation field, heat source attribute, thermophysical property parameters, and mechanical parameters to each of the computational domain meshes after segmentation to obtain an attribute mesh; A boundary condition adding module, used for adding boundary conditions to the fluid field, the heat transfer field, and the deformation field of each attribute grid to obtain a simulation grid; A reference temperature adding module, configured to add a reference temperature to the simulation grid according to the distribution position of the simulation grid to obtain a thermal-fluid-solid coupling model; A simulation calculation module is used to simulate the heat-fluid-solid coupling model under different heat source conditions to obtain the temperature field and the strain field in the crust; The correlation analysis module is used to use GIS spatial analysis technology to perform superposition comparison and spatial coupling relationship analysis on the spatial distribution data of the temperature field in the crust, the strain field in the crust, and surface geothermal resources and underground earthquake disasters, so as to obtain the analysis results of the geothermal and earthquake coupling relationship.

9. An electronic device, characterized in that: include: At least one processor and a memory communicatively connected to the processor; wherein the memory stores instructions executable by the processor, and the instructions are executed by the processor so that the processor can execute a geothermal and seismic symbiosis mechanism analysis method according to any one of claims 1 to 7.

10. A non-transitory computer-readable storage medium storing computer instructions, characterized in that: The computer instructions are used to enable a computer to execute a geothermal and seismic symbiotic mechanism analysis method according to any one of claims 1 to 7.

Images