Underground medium heat transfer model-based bottom heat flow and deep layer temperature inversion method

Abstract

The invention discloses a bottom heat flow and deep layer temperature inversion method based on an underground medium heat transfer model, and relates to the field of deep stratum temperature field prediction. The underground medium heat transfer model-based bottom heat flow and deep temperature inversion method mainly comprises the following steps of: constructing an electromagnetic data three-dimensional axial anisotropy inversion parameter target functional to carry out fine inversion on an electromagnetic data electrical structure of a target area to obtain an electrical structure model and an underground heat storage stratum structure model of the target area; a deep stratum heat transfer model is obtained in combination with regional geothermal geological information; and utilizing a Newton-Raphson optimization algorithm to perform inversion on the bottom prior heat flow value distribution boundary condition to obtain temperature field distribution of the deep stratum heat transfer model. According to the bottom heat flow and deep layer temperature inversion method based on the underground medium heat transfer model, the electrical characteristics of the underground medium can be converted into temperature field distribution based on a small number of borehole logging constraints, the prediction precision is high, and the prediction range is wide.

Description

Technical Field

[0001] The present invention relates to the field of deep stratum temperature field prediction, and more specifically, to a bottom heat flow and deep temperature inversion method based on an underground medium heat transfer model. Background Art

[0002] The accuracy of deep temperature prediction directly affects the accuracy of geothermal resource exploration and evaluation and the benefits of development and utilization. Therefore, it is extremely important to estimate the distribution characteristics of the underground temperature field as accurately as possible.

[0003] At present, there are two main ways to obtain the temperature inside the earth: direct measurement and indirect detection. The first type of direct measurement method is mainly based on the existing borehole temperature interpolation to obtain the regional temperature field, but if the number of borehole temperature measurements is too small, it often leads to considerable errors, especially in areas with complex geological structures; and if there are too many boreholes, although the temperature field prediction accuracy can be improved to a certain extent, the cost is too high. The second type of indirect detection method is mainly based on geochemical or geophysical means to indirectly predict the geothermal field. Among them, the geochemical method mainly uses geochemical geothermometers to predict the temperature field, that is, the temperature characteristics are inferred from various empirical formulas such as geochemical temperature scales. This type of indirect geothermometer can predict the distribution range of heat reservoir temperature, but it cannot provide regional temperature distribution and it is difficult to match the predicted temperature with depth. The geophysical method is mainly based on the relationship between geophysical parameters and temperature to infer the temperature field, that is, the empirical relationship formula between physical properties and temperature is used to estimate the formation temperature, but the parameters in this type of empirical formula are assumed to be constant in any geological environment and depth, which often causes large prediction errors.

[0004] Therefore, in the existing temperature field prediction methods, it is difficult to balance cost and accuracy for borehole interpolation; it is difficult for geochemical temperature scales to accurately match temperature with depth; and it is difficult for geophysical indirect geothermometers with purely empirical formulas to accurately estimate deep stratum temperatures. Therefore, the existing temperature estimation methods are difficult to accurately predict the temperature distribution at depths not reached by the borehole and in the space between wells, and it is even more difficult to accurately provide regional deep geothermal field distribution characteristics. Summary of the invention

[0005] The purpose of the present invention is to provide a bottom heat flow and deep temperature inversion method and computer equipment based on an underground medium heat transfer model, which can convert the electrical characteristics of the underground medium into temperature field distribution based on a small amount of borehole logging constraints, with high prediction accuracy and a wide prediction range.

[0006] The present invention provides a bottom heat flow and deep temperature inversion method based on an underground medium heat transfer model, which is characterized by comprising the following steps: Step 1: constructing a three-dimensional axial anisotropy inversion parameter target functional of electromagnetic data according to electromagnetic detection data and prior information of a target area; Step 2: performing a fine inversion of the electrical structure of electromagnetic data of the target area according to the three-dimensional axial anisotropy inversion parameter target functional of the electromagnetic data to obtain an electrical structure model; Step 3: obtaining an underground heat reservoir formation structural model of the target area according to the electrical structure model, regional geothermal geology and formation information and prior information on thermophysical parameters of different blocks in the deep formation of the target area; Step 4: The underground heat reservoir stratum structural model of the target area is combined with regional geothermal geological information to obtain a deep stratum heat transfer model; the deep stratum heat transfer model includes boundary conditions, and the boundary conditions include upper surface average temperature boundary conditions, bottom prior heat flow value distribution boundary conditions and adiabatic boundary conditions on both sides of the deep stratum heat transfer model section; Step 5: According to the prior information of thermophysical property parameters of different blocks in the deep stratum of the target area, the upper surface average temperature boundary conditions and the adiabatic boundary conditions on both sides of the deep stratum heat transfer model section, the Newton-Raphson optimization algorithm is used to invert the bottom prior heat flow value distribution boundary conditions to obtain the temperature field distribution of the deep stratum heat transfer model.

[0007] Furthermore, step 1 specifically includes: according to the electromagnetic detection data and prior information of the target area, a three-dimensional axial anisotropy model functional is used to constrain the data target functional at the same time, and a three-dimensional axial anisotropy inversion parameter target functional of the electromagnetic data is constructed, such as the formula: , in, Represents the target functional of the three-dimensional axial anisotropy inversion parameters of electromagnetic data; represents the model parameters; It is a data mismatch; and are the weights of the smoothness constraint term and the anisotropy penalty term respectively; and They are the smoothness constraint term and the anisotropy penalty term respectively; is the data covariance matrix; Represents a given model m forecast data; represents observation data; is the model weight matrix; is the reference model; is the number of electromagnetic data observation points, , and The inversion units are , and The principal resistivity in the direction.

[0008] Furthermore, step 2 specifically includes: Use the Gauss-Newton method to minimize the objective function, as shown in the formula: , in, is the data covariance matrix; denotes the Jacobian matrix; superscript Represents the transpose of a matrix; is the weight of the smooth constraint; is the model weight matrix; The weight of the anisotropy penalty term; represents the Hessian matrix in the anisotropic penalty term; represents the model update amount, Represents a given model forecast data; represents observation data; represents the weight of the smooth constraint; represents the model parameters; represents the reference model; Use Armijo linear search to get the step size and update the electrical structure model, such as the formula: , in, Represents the next iteration model; Represents the current iteration model; Indicates the search step length; Indicates the model update amount; The electrical structure model is corrected by iterative inversion so that the overall root mean square error between the final inversion result and the measured data is less than or equal to the preset fitting error value, and the electrical structure model is obtained, as shown in the formula: , in, is the overall root mean square error, represents the maximum number of iterations, and Respectively represent Observation data and Forward modeling data.

[0009] Furthermore, step Step3 specifically includes: Step31: obtaining the underground stratum structure distribution characteristics according to the electrical structure model and the regional geothermal geology and stratum information; Step32: obtaining the thermal physical parameters of different strata in the target area according to the underground stratum structure distribution characteristics and the prior information of the thermal physical parameters of different blocks in the deep strata of the target area; Step33: obtaining the underground thermal reservoir stratum structure model of the target area according to the thermal physical parameters of different strata in the target area.

[0010] Furthermore, the thermophysical parameters of different formations in the above-mentioned target area include formation density, thermal conductivity and specific heat capacity.

[0011] Furthermore, step 5 includes: according to the prior information of thermophysical parameters of different blocks in the deep strata of the target area, the upper surface average temperature boundary condition, the adiabatic boundary conditions on both sides of the deep stratum heat transfer model section and the bottom prior heat flow value distribution boundary condition, taking the borehole temperature measurement as the response constraint condition, constructing the bottom boundary heat flow distribution optimization inversion objective function, such as the formula: , in, is the fitting degree, i.e., the inverse optimization objective function; is the regression sum of squares, is the total sum of squares of deviations, is the number of measured temperature data, is the inversion temperature, that is Temperature field of the heat transfer model during the inversion process Temperature at the verification borehole; is the average measured temperature, is the residual sum of squares; To predict temperature.

[0012] Furthermore, step 5 also includes: using the Newton-Raphson optimization algorithm to generate an initialized random population, such as the formula: , , in, Indicates The individual Dimensional location; express A random number between Indicates the lower boundary; Indicates the upper boundary; is the population matrix of all dimensions; represents the population size, Denotes the dimension that represents the problem.

[0013] Furthermore, step 5 also includes: using the Newton-Raphson optimization algorithm to update the iterative population, such as formula: , , , , in Newton-Raphson optimization algorithm The main components of Represents a normally distributed random number with a mean of 0 and a variance of 1; Represents a A random number for each decision variable; and is used and The positions of the two vectors generated; is the current variable value, and are distinct integers randomly selected from the population, Represents the average value of a set of data. is the process quantity, is the population matrix, The worst position, For the best position; To update The new vector position obtained is is the current vector position, and is a random number between 0 and 1. and are all random position vectors, For consideration , and The new vector position obtained is For consideration , and The new vector position obtained is is the adaptive coefficient, Indicates the current iteration, Indicates the maximum number of iterations.

[0014] Furthermore, step 5 also includes: using the trap avoidance operator algorithm to further optimize the population, such as formula: , , , , in, is the position vector of the nth individual in the population matrix of the next iteration, To obtain the position vector for the next iteration using the trap avoidance operator algorithm, and They are and A uniform random number between and is a random number, For the best position, is the population matrix of the current iteration The position vector of each individual, is the adaptive coefficient, is the average value of a set of data, is the population matrix of the current iteration, For control The determining factor of performance is Represents a binary number 1 or 0. for A uniform random number between .

[0015] The present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, the steps of the bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model are implemented.

[0016] The bottom heat flow and deep temperature inversion method and computer equipment based on the underground medium heat transfer model provided by the present invention have the following beneficial effects: The present invention is based on the joint constrained inversion of borehole temperature measurement-electromagnetic detection-intelligent algorithm, combines the electromagnetic detection data of the target area with the logging data, and carries out the anisotropic electrical structure inversion of electromagnetic data based on the logging constraints; wherein, the three-dimensional axial anisotropy model functional is used to constrain the data target functional at the same time, and the three-dimensional axial anisotropy inversion parameter target functional of the electromagnetic data is established, and the electrical structure of the target area is finely inverted based on the pre-processed nonlinear conjugate gradient algorithm, which can greatly improve the accuracy of the electrical structure model; The present invention is based on converting the inverted underground stratum structure distribution into the distribution of thermophysical parameters of different strata in the target area, integrating the geothermal geological information of the target area to construct a deep medium heat transfer model; integrating the thermophysical parameters of different strata in the target area, the upper surface average temperature boundary conditions, and the adiabatic boundary conditions on both sides of the deep stratum heat transfer model section, and using the Newton-Raphson optimization algorithm to achieve the bottom heat flow value distribution and the inversion of the deep stratum temperature field; The present invention can solve the problem that existing temperature estimation methods are difficult to accurately predict the temperature distribution at depths not reached by the borehole and in the space between wells, and are difficult to accurately provide regional deep geothermal field distribution characteristics. The present invention can convert the electrical characteristics of underground media into temperature field distribution based on a small number of borehole logging constraints, with high prediction accuracy and a wide prediction range. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] The present invention will be further described below with reference to the accompanying drawings and embodiments, in which: Figure 1 It is a flow chart of the bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model provided by the present invention; Figure 2 It is the inversion resistivity profile of the deep formation in the target area provided by the present invention; Figure 3 It is the deep stratigraphic structure and structural interpretation map of the target area provided by the present invention; Figure 4 It is the inversion distribution diagram of the bottom heat flow value of the target area heat transfer model provided by the present invention; Figure 5 It is the inversion profile of the temperature field of deep formation in the target area provided by the present invention; Figure 6 It is a comparison and verification diagram of the target area inversion temperature and the measured temperature provided by the present invention; Figure 7 It is a structural diagram of a deep geothermal field prediction device based on borehole temperature measurement-electromagnetic detection-intelligent algorithm joint constrained inversion provided by the present invention; Figure 8 It is a structural block diagram of the computer device provided by the present invention. DETAILED DESCRIPTION

[0018] In order to have a clearer understanding of the technical features, purposes and effects of the present invention, specific embodiments of the present invention are now described in detail with reference to the accompanying drawings.

[0019] Figure 1 A schematic diagram of the bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model of this embodiment is shown. In this embodiment, the bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model includes the following steps: Step 1: Based on the electromagnetic detection data and prior information of the target area, construct the target functional of the three-dimensional axial anisotropy inversion parameters of the electromagnetic data; In an exemplary embodiment, step Step 1 specifically includes: according to the electromagnetic detection data and prior information of the target area, a three-dimensional axial anisotropy model functional is used to constrain the data target functional at the same time, and a three-dimensional axial anisotropy inversion parameter target functional of the electromagnetic data is constructed, such as the formula: (1) in, Represents the target functional of the three-dimensional axial anisotropy inversion parameters of electromagnetic data; represents the model parameters; It is a data mismatch; and are the weights of the smoothness constraint term and the anisotropy penalty term respectively; and They are the smoothness constraint term and the anisotropy penalty term respectively; is the data covariance matrix; Represents a given model m forecast data; represents observation data; is the model weight matrix; is the reference model; is the number of electromagnetic data observation points, , and The inversion units are , and principal resistivity in direction; It should be noted that the target functional in Step 1 is the target functional for the three-dimensional axial anisotropy inversion of electromagnetic data; Step 2: According to the target functional of the three-dimensional axial anisotropy inversion parameter of the electromagnetic data, a fine inversion of the electrical structure of the electromagnetic data of the target area is performed to obtain an electrical structure model; In an exemplary embodiment, step Step 2 specifically includes: Use the Gauss-Newton method to minimize the objective function, as shown in the formula: (2) in, is the data covariance matrix; denotes the Jacobian matrix; superscript Represents the transpose of a matrix; is the weight of the smooth constraint; is the model weight matrix; The weight of the anisotropy penalty term; represents the Hessian matrix in the anisotropic penalty term; represents the model update amount, Represents a given model forecast data; represents observation data; represents the weight of the smooth constraint; represents the model parameters; represents the reference model; Use Armijo linear search to get the step size and update the electrical structure model, such as the formula: (3) in, Represents the next iteration model; Represents the current iteration model; Indicates the search step length; Indicates the model update amount; The electrical structure model is corrected by iterative inversion so that the overall root mean square error between the final inversion result and the measured data is less than or equal to the preset fitting error value, and the electrical structure model is obtained, as shown in the formula: (4) in, is the overall root mean square error, represents the maximum number of iterations, and Respectively represent Observation data and Forward modeling data; For an exemplary embodiment, in step Step 2, based on the established inversion parameter target functional, a nonlinear conjugate gradient algorithm based on preprocessing is used to solve, and an Armijo line search is used to obtain a step-size update model to achieve fine inversion of the electrical structure of the electromagnetic data in the target area; specifically, based on the established inversion parameter target functional, a nonlinear conjugate gradient algorithm based on preprocessing is used to solve, and an Armijo line search is used to obtain a step-size update model to achieve fine inversion of the electrical structure of the electromagnetic data in the target area, so that the corrected electrical structure model has the best fit with the prior geoelectric model, that is, the overall fitting error between the final inversion result and the measured data is less than or equal to the preset fitting error value ; Realize the fine inversion of anisotropic electrical structure of electromagnetic data in the target area; In an exemplary embodiment, step 2 is specifically as follows: based on the anisotropic inversion target functional of the electromagnetic data constructed in step 1, combined with the existing prior electrical structure model, the model is iteratively updated using formula (2) and formula (3), and the electrical structure model is continuously inverted and corrected, so that the corrected electrical structure model has the best fit with the prior geoelectric model, and the root mean square error RMS (obtained using formula (4)) between the final inversion result and the measured data is 0.0065, which is less than the preset fitting error value , the electrical distribution obtained at this time is the accurate electrical structure model of the target area, such as Figure 2 As shown; Step 3: Obtaining the underground thermal reservoir stratigraphic structure model of the target area based on the electrical structure model, regional geothermal geology and stratigraphic information, and prior information on thermal physical property parameters of different blocks in the deep stratigraphic layers of the target area; In an exemplary embodiment, step Step 3 specifically includes: Step 31: Obtaining underground stratum structural distribution characteristics according to the electrical structure model and regional geothermal geology and stratum information; Step 32: obtaining the thermal physical property parameters of different strata in the target area according to the structural distribution characteristics of the underground strata and the prior information of the thermal physical property parameters of different blocks in the deep strata of the target area; Step 33: obtaining an underground heat reservoir formation structure model of the target area according to the thermal physical property parameters of different formations in the target area; In an exemplary embodiment, the thermophysical property parameters of different formations in the target area include formation density, thermal conductivity and specific heat capacity; For an exemplary embodiment, in step 3, based on the electrical structure model after fine inversion, the regional geological and stratigraphic information is integrated to convert the inverted electrical structure into the stratigraphic structure and structural characteristics of the target area, and define the thermophysical parameters of different stratigraphic structures and structural areas to construct the underground thermal reservoir stratigraphic structural model of the target area; For an exemplary embodiment, step 3 includes: based on the fine inversion of the anisotropic electrical structure of electromagnetic data, the deep electrical structure model of the target area is obtained, the regional geothermal geology and stratigraphic information constraints are integrated, and the structural distribution characteristics of the underground stratigraphic structure (such as Figure 3 As shown); combined with the prior information of thermal physical parameters of different blocks in the deep strata of the target area, the underground strata structural distribution characteristics interpreted by the electrical structure are converted into thermal physical parameters of different strata in the target area, including stratum density ρ , thermal conductivity k and specific heat c The distribution of thermophysical parameters (as shown in Table 1); the above parameters are used to construct the underground thermal reservoir stratigraphic structural model of the target area; Table 1: Table of thermal attribute parameter assignments for each layer of the target area heat transfer model

[0020] It should be noted that in step 3, the underground stratum structural distribution characteristics include deep strata, faults and structures; Step 4: According to the underground heat reservoir formation structure model of the target area, combined with the regional geothermal geological information, a deep formation heat transfer model is obtained; the deep formation heat transfer model includes boundary conditions, which include the upper surface average temperature boundary condition, the bottom prior heat flow value distribution boundary condition and the adiabatic boundary conditions on both sides of the deep formation heat transfer model section; For an exemplary embodiment, in step 4, the heat transfer mode of the underground medium in the target area is analyzed in combination with the profile of the underground heat reservoir formation structure model of the target area and the regional geothermal geological information, the underground heat transfer process of the target area is obtained, and the deep formation heat transfer model is constructed using the heat conduction control equation; For an exemplary embodiment, in step 4, a deep formation heat transfer model is constructed using formula (a4) to formula (a13); For an exemplary embodiment, in step 4, the boundary conditions of the deep formation heat transfer model are defined; The boundary conditions of the deep stratum heat transfer model include: the upper surface average temperature boundary condition , bottom prior heat flux value distribution boundary conditions , and the adiabatic boundary conditions on both sides of the deep formation heat transfer model section ; In this embodiment, the upper surface average temperature boundary condition is the temperature at a depth of 30m in the constant temperature zone of the target area. = 20.0°C, adiabatic boundary conditions on both sides of the deep formation heat transfer model section The bottom prior heat flow value distribution boundary condition is taken as the prior value at each grid node of the deep formation heat transfer model (the target area is meshed during modeling). ; Combined with the thermal reservoir formation density based on the electrical structure model conversion , thermal conductivity and specific heat Parameter distribution (as shown in Table 1); finite element numerical simulation of the temperature field of the heat transfer model Forward response; where is the heat flux value distribution of the bottom boundary condition; is the prior value of the heat flux distribution for the bottom boundary condition; In this embodiment, the forward heat transfer model of the temperature field, the thermal property parameters of the model, and the initial boundary conditions of the model can be used to solve the heat transfer model section in the initial heat transfer model. Temperature field under the condition Forward response; Among them, the differential equation of the underground space temperature field is: (a4) in, is the underground space temperature; The differential equation of the underground temperature field is multiplied by the impulse function And integral: (a5) According to the Hamiltonian operator operation rules and Ogaard formula in field theory, and substituting the boundary conditions, the integral equation of the underground space temperature field is obtained as follows: (a6) in, For the target area, is the upper boundary of the deep formation heat transfer model, is the lower boundary of the deep formation heat transfer model, is the heat flow value, Indicates thermal conductivity, unit , is the temperature, unit , is the divergence, is the gradient, is the impulse function; The integral equation of the underground space temperature field is discretized and linearly interpolated to obtain the grid unit integral of each term in the integral equation: (a7) (a8) (a9) in, is the array of impulse functions in the grid cells, is the array of temperatures within the cell, , is the stiffness matrix within the element, and and the upper boundary condition surface temperature And the bottom heat flow value distribution of the lower boundary condition The source vector in the relevant element is solved to obtain the stiffness matrix and the coefficients in the source vector: (a10) (a11) (a12) (a13) in, , , , The stiffness matrices are , , source vector , Elements within; , , Label the nodes; , , The subunits are , , The length of the side in the direction, in ; Indicates thermal conductivity, unit ; , , For Node Coordinates on the parent unit, indicating the relationship with the child unit , , The direction of the parent unit with the same direction, dimensionless, where the parent unit refers to the local coordinates of the hexahedral unit in the hexahedral mesh of the finite element; the node represents the vertex of the unit; is the upper boundary temperature, in units o C; is the heat flow value, unit ; According to the number of each node of each unit and extended to the element position corresponding to the target area, the linear equations of the synthesized overall stiffness matrix and source vector are obtained: (a14) in, is the target area impulse function array, is the temperature array in the target area, , is the stiffness matrix in the target area, , are all source vectors in the target area, ignoring , after simplification: (a15) in, , ; After solving the linear equations After finite element numerical simulation, the temperature field of underground space is obtained , , Respectively , Axis direction; Step 5: Based on the prior information of thermophysical property parameters of different blocks in the deep formation of the target area, the upper surface average temperature boundary condition and the adiabatic boundary conditions on both sides of the deep formation heat transfer model section, the Newton-Raphson optimization algorithm is used to invert the bottom prior heat flow value distribution boundary condition to obtain the temperature field distribution of the deep formation heat transfer model; It should be noted that in step 5, the bottom boundary heat flux distribution is constructed based on the borehole temperature measurement as the response constraint condition. q down Optimize the inversion objective function (Formula (16)) and continuously modify the iterative model based on the Newton-Raphson optimization algorithm ; For an exemplary embodiment, in step 5, the thermal physical parameters of different formations in the target area obtained in step 3 (including the formation density of the heat transfer model) are integrated. ρ , thermal conductivity k and specific heat c Parameters, as shown in Table 1), the upper surface average temperature boundary condition, the adiabatic boundary conditions on both sides of the deep stratum heat transfer model section, and the temperature field of the heat transfer model in step 4 Forward response algorithm; use formula (5) to construct the bottom boundary heat flux distribution with borehole temperature measurement as the response constraint condition Optimize the inversion objective function; specifically, use the Newton-Raphson-based optimizer (NRBO) algorithm to iteratively invert the bottom optimal heat flux , so that the boundary condition of the optimal bottom heat flux value is The forward response of the temperature at the drill hole on the heat transfer model section is shown below. The measured temperature The highest degree of agreement, and The fitting error is less than or equal to the preset value , the temperature field response obtained at this time That is, the distribution characteristics of the deep geothermal field in the target area; In this embodiment, using Iterative inversion is used to obtain the optimal bottom heat flow value distribution; in the inverted optimal bottom heat flow value distribution (such as Figure 4 The forward modeling response of the temperature at the drill hole on the heat transfer model section under the condition shown in The measured temperature The best fit is obtained. At this time, the optimal inversion objective function (5) is calculated as and The goodness of fit (GOF) value of the two is 0.96, which is greater than the preset value. Finally, the distribution characteristics of the deep temperature field in the target area can be obtained based on the finite element temperature field response forward algorithm in Step 4 (such as Figure 5 As shown); it should be pointed out here that the constraint borehole for the bottom heat flow value and temperature field inversion is D1; ​​the verification borehole for the inversion temperature is D2; In an exemplary embodiment, step 5 includes: according to the prior information of thermophysical parameters of different blocks in the deep formation of the target area, the upper surface average temperature boundary condition and the adiabatic boundary conditions on both sides of the deep formation heat transfer model section and the bottom prior heat flow value distribution boundary condition, with the borehole temperature measurement as the response constraint condition, constructing the bottom boundary heat flow distribution optimization inversion objective function, such as the formula: (5) in, is the fitting degree, i.e., the inverse optimization objective function; is the regression sum of squares, is the total sum of squares of deviations, is the number of measured temperature data, is the inversion temperature, that is Temperature field of the heat transfer model during the inversion process Temperature at the verification borehole; is the average measured temperature, is the residual sum of squares; To predict temperature; It should be noted that the bottom heat flux distribution inversion in Step 5 is based on the Newton-Raphson optimization The algorithm is used to obtain the optimal bottom heat flow value in the inversion. Forward modeling response of temperature at the borehole on the profile of deep formation heat transfer model under distributed boundary conditions The measured temperature The degree of fit is greater than or equal to a preset degree of fit value; It should be noted that the GOF value range is 0-1. The closer the GOF value is to 1, the better the fit between the predicted value and the measured value. In an exemplary embodiment, step Step 5 further includes: using the Newton-Raphson optimization algorithm to generate an initial random population, such as formula: (6) (7) in, Indicates The individual Dimensional location; express A random number between Indicates the lower boundary; Indicates the upper boundary; is the population matrix; represents the population size, Represents the dimension of the problem; In an exemplary embodiment, step Step 5 further includes: using the Newton-Raphson optimization algorithm to update the iterative population, such as formula: , , , , in, It is the main component of the Newton-Raphson optimization algorithm. Represents a normally distributed random number with a mean of 0 and a variance of 1; Represents a A random number for each decision variable; and is used and The positions of the two vectors generated; is the current variable value, and are distinct integers randomly selected from the population, Represents the average value of a set of data. is the process quantity, is the population matrix, The worst position, For the best position; To update The new vector position obtained is is the front vector position, and is a random number between 0 and 1. and are all random position vectors, For consideration , and The new vector position obtained is For consideration , and The new vector position obtained is is the adaptive coefficient, Indicates the current iteration, Indicates the maximum number of iterations; For an exemplary embodiment, in step 5, the inverted bottom heat flow value distribution based on the Newton-Raphson optimization algorithm requires the Newton-Raphson search algorithm to update the iterative population, including: Use Taylor's formula to determine the second derivative: (8) Subtracting or adding the above formula, we can get and The expression is: (9) The updated position is redefined as follows: (10) The Newton-Raphson optimization algorithm converts the position Replace with , identify The position in its neighborhood is better than the position Good; Location By location replace, The position in its neighborhood is better than the position bad; use position ( ) rather than its fitness , recommended It is expressed as: (11) In the formula, Represents a normally distributed random number with a mean of 0 and a variance of 1; The worst position, For the best position; is a random number generated within the dimension of the decision variable; Based on formula (11) , modify formula (10) into the following form: (12) By adding another parameter To improve the use of the Newton-Raphson search algorithm, this parameter guides the population in the right direction. The expression is as follows: (13) in, and is a random number between 0 and 1. and are distinct integers randomly selected from the population; Combining formulas (12) and (13), the vector The current position of is updated by: (14) in By updating The new vector position obtained; in formula (10) The Newton-Raphson method can be further improved to: (15) in and is used and The positions of the two vectors generated; r 1 means and are distinct integers randomly selected from the population; Improved Given in equation (15), using equation (15), equation (14) is updated as follows: (16) The best vector in formula (16) Replace the current vector with The position of to construct a new vector ; (17) The development phase is the main focus of this search direction strategy; in terms of local search, the search method proposed by equation (17) is effective, but has limitations in global search, while the search strategy proposed by equation (16) is effective in global search, but has limitations in local search; however, Equations (16) and (17) are used simultaneously to improve the diversification and reinforcement stages; in the next iteration, the new position vector is expressed in equation (18); (18) in are distinct integers randomly selected from the population; Indicates the current iteration, Indicates the maximum number of iterations; The value of varies between 1 and -1; is the adaptive coefficient of the algorithm, that is, as a rule of thumb, the optimization algorithm must be able to strike a balance between diversification and intensiveness in order to find the optimal solution in the search space and eventually converge to the global solution. Combined with the adaptive parameter , the proposed By introducing random operations during the optimization process, the The diversity of ,avoids falling into the local optimum, and significantly reduces the number of iterations; In an exemplary embodiment, step Step 5 further includes: using a trap avoidance operator algorithm to further optimize the population, such as formula: , , , , in, is the position vector of the nth individual in the population matrix of the next iteration, To obtain the position vector for the next iteration using the trap avoidance operator algorithm, and They are and A uniform random number between and is a random number, For the best position, is the population matrix of the current iteration The position vector of each individual, is the adaptive coefficient, is the average value of a set of data, is the population matrix of the current iteration, For control The determining factors of performance are Represents a binary number 1 or 0. for A uniform random number between ; It should be noted that in Step 5, based on Newton-Raphson optimization The bottom heat flux distribution inversion algorithm needs to use the Trap Avoidance Operator (TAO) algorithm to further optimize the population, including: using the Trap Avoidance Operator to significantly change position, by combining the best position and the current vector position , producing a high-quality solution ; In the above formula, and The population is more diverse, freeing itself from the constraints of the local optimal solution, and improving the diversity of the population. The following equations are used to generate the diversity of the population. If the value of Δ is greater than or equal to 0.5, then The value of is 0; otherwise, the value is 1; For an exemplary embodiment, in Step 5, under the optimal bottom heat flow value distribution boundary condition, the forward response of the temperature at the borehole on the deep formation heat transfer model section is The measured temperature Goodness of fit (GOF) is greater than or equal to the preset fit value ; It should be noted that the inversion temperature uses the goodness of fit (GOF) of formula (5) to measure the deviation between the measured value and the predicted value. The larger the value, the more reliable the prediction result. If the data fit is greater than the preset value ( ), the iteration stops and the current prediction model parameters are output as the optimal prediction result; otherwise, the update is returned. And solve the new prediction model correction; Figure 6 The comparison and verification diagram of the predicted temperature and the measured temperature of the target area provided in this embodiment; In this embodiment, the temperature variation curve with depth at the borehole temperature measurement point in the predicted temperature profile of the target area is extracted and compared with the actual logging temperature; The experimental results show that the D2 borehole (D2 borehole position is as shown in FIG. Figure 5 The goodness of fit between the predicted temperature and the measured temperature is , which shows that the prediction accuracy of the temperature field in this embodiment is 96.57%, and the predicted temperature is highly consistent with the measured temperature.

[0021] like Figure 7 As shown, this embodiment provides a deep geothermal field prediction device, including the following modules: The anisotropic electrical structure inversion module is used to use the three-dimensional axial anisotropy model functional to constrain the data target functional at the same time, establish the three-dimensional axial anisotropy inversion parameter target functional of electromagnetic data, solve it using the preprocessing-based nonlinear conjugate gradient algorithm, use the Armijo line search to obtain the step size update model, and combine the geoelectric section prior information constraints to achieve fine inversion of the electrical structure of the electromagnetic data in the target area; The heat reservoir formation model construction module is used to integrate the geothermal geological information of the target area based on the electrical model after anisotropic structural inversion, convert the underground formation structure and structural characteristics interpreted by the inversion electrical structure model into the distribution of thermophysical parameters of different formations, and construct the underground heat reservoir formation structural model of the target area; The temperature field forward modeling module of the heat transfer model is used to obtain the heat transfer process of the underground medium in the target area and establish a deep formation heat transfer model by combining the underground heat reservoir formation structure model section of the target area and the prior geothermal geological information constraints; the boundary conditions of the deep formation heat transfer model include: the upper surface average temperature boundary condition, the bottom heat flow value distribution boundary condition, and the left and right adiabatic boundary conditions; The bottom heat flow value inversion module is used to integrate the thermophysical parameters of different formations in the target area, the upper surface average temperature boundary conditions, and the left and right adiabatic boundary conditions, to construct the bottom heat flow value distribution inversion objective function with borehole temperature measurement as the response constraint condition, and to use the optimization algorithm to invert the optimal bottom heat flow value distribution boundary conditions; The deep formation temperature field calculation module is used to calculate the temperature field response of the deep formation heat transfer model based on the optimal bottom heat flow value distribution boundary conditions, combined with the thermal physical property distribution characteristics of the underground formation, using the deep formation heat transfer model and forward algorithm.

[0022] This embodiment provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, the steps of the bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model are implemented.

[0023] like Figure 8 As shown, the computer device 120 may include: at least one processor 121, such as a central processing unit (CPU), at least one communication interface 123, a memory 124, and at least one communication bus 122. Among them, the communication bus 122 is used to realize the connection and communication between these components. Among them, the communication interface 123 may include a display screen (Display), a keyboard (Keyboard), and the optional communication interface 123 may also include a standard wired interface and a wireless interface. The memory 124 may be a high-speed random access memory (Random Access Memory, RAM), or a non-volatile memory (non-volatile memory), such as at least one disk storage. The memory 124 may also be at least one storage device located away from the aforementioned processor 121. Among them, the memory 124 stores an application program, and the processor 121 calls the program code stored in the memory 124 to perform any of the above method steps. Among them, the communication bus 122 may be a peripheral component interconnect (PCI) bus or an extended industry standard architecture (EISA) bus, etc. The communication bus 122 can be divided into an address bus, a data bus, a control bus, etc. For ease of representation, Figure 8Only one line is used to represent it, but it does not mean that there is only one bus or one type of bus. Among them, the memory 124 may include a volatile memory (volatile memory), such as a random-access memory (random-access memory, RAM); the memory may also include a non-volatile memory (non-volatile memory), such as a flash memory (flash memory), a hard disk drive (hard disk drive, HDD) or a solid-state drive (solid-state drive, SSD); the memory 124 may also include a combination of the above-mentioned types of memory. Among them, the processor 121 may be a central processing unit (central processing unit, CPU), a network processor (network processor, NP) or a combination of CPU and NP. Among them, the processor 121 may further include a hardware chip. The above-mentioned hardware chip may be an application-specific integrated circuit (application-specific integrated circuit, ASIC), a programmable logic device (programmable logic device, PLD) or a combination thereof. The above-mentioned PLD may be a complex programmable logic device (complex programmable logic device, CPLD), a field-programmable gate array (field-programmable gate array, FPGA), a generic array logic (generic array logic, GAL) or any combination thereof. Optionally, the memory 124 is also used to store program instructions. The processor 121 can call the program instructions to implement the bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model as in this embodiment.

[0024] The embodiments of the present invention are described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific implementation modes, which are merely illustrative rather than restrictive. Under the guidance of the present invention, ordinary technicians in this field can also make many forms without departing from the scope of protection of the present invention and the claims, all of which are within the protection of the present invention.

Claims

1. A bottom heat flow and deep temperature inversion method based on underground medium heat transfer model, characterized in that: The following steps are involved: Step 1: Based on the electromagnetic detection data and prior information of the target area, construct the target functional of the three-dimensional axial anisotropy inversion parameters of the electromagnetic data; Step 2: According to the target functional of the three-dimensional axial anisotropy inversion parameter of the electromagnetic data, a fine inversion of the electrical structure of the electromagnetic data of the target area is performed to obtain an electrical structure model; Step 3: Obtaining the underground thermal reservoir stratigraphic structure model of the target area based on the electrical structure model, regional geothermal geology and stratigraphic information, and prior information on thermal physical property parameters of different blocks in the deep stratigraphic layers of the target area; Step 4: According to the underground heat reservoir formation structure model of the target area, combined with the regional geothermal geological information, a deep formation heat transfer model is obtained; the deep formation heat transfer model includes boundary conditions, which include the upper surface average temperature boundary condition, the bottom prior heat flow value distribution boundary condition and the adiabatic boundary conditions on both sides of the deep formation heat transfer model section; Step 5: Based on the prior information of thermophysical parameters of different blocks in the deep formation of the target area, the upper surface average temperature boundary conditions and the adiabatic boundary conditions on both sides of the deep formation heat transfer model profile, the Newton-Raphson optimization algorithm is used to invert the bottom prior heat flow value distribution boundary conditions to obtain the temperature field distribution of the deep formation heat transfer model.

2. The bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model according to claim 1 is characterized in that: Step 1 specifically includes: according to the electromagnetic detection data and prior information of the target area, a three-dimensional axial anisotropy model functional is used to constrain the data target functional at the same time, and a three-dimensional axial anisotropy inversion parameter target functional of the electromagnetic data is constructed, such as the formula: ， in, Represents the target functional of the three-dimensional axial anisotropy inversion parameters of electromagnetic data; represents the model parameters; It is a data mismatch; and are the weights of the smoothness constraint term and the anisotropy penalty term respectively; and They are the smoothness constraint term and the anisotropy penalty term respectively; is the data covariance matrix; Represents a given model m forecast data; represents observation data; is the model weight matrix; is the reference model; is the number of electromagnetic data observation points, , and The inversion units are , and The principal resistivity in the direction.

3. The bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model according to claim 1 is characterized in that: Step 2 specifically includes: Use the Gauss-Newton method to minimize the objective function, as shown in the formula: ， in, is the data covariance matrix; denotes the Jacobian matrix; superscript Represents the transpose of a matrix; is the weight of the smooth constraint; is the model weight matrix; The weight of the anisotropy penalty term; represents the Hessian matrix in the anisotropic penalty term; represents the model update amount, Represents a given model forecast data; represents observation data; represents the weight of the smooth constraint; represents the model parameters; represents the reference model; Use Armijo linear search to get the step size and update the electrical structure model, such as the formula: ， in, Represents the next iteration model; Represents the current iteration model; Indicates the search step length; Indicates the model update amount; The electrical structure model is corrected by iterative inversion so that the overall root mean square error between the final inversion result and the measured data is less than or equal to the preset fitting error value, and the electrical structure model is obtained, as shown in the formula: ， in, is the overall root mean square error, represents the maximum number of iterations, and Respectively represent Observation data and Forward modeling data.

4. The bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model according to claim 1 is characterized in that: Step 3 specifically includes: Step 31: Obtaining underground stratum structural distribution characteristics according to the electrical structure model and regional geothermal geology and stratum information; Step 32: obtaining the thermal physical property parameters of different strata in the target area according to the structural distribution characteristics of the underground strata and the prior information of the thermal physical property parameters of different blocks in the deep strata of the target area; Step 33: According to the thermophysical parameters of different formations in the target area, the underground heat reservoir formation structure model of the target area is obtained.

5. The bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model according to claim 4 is characterized in that: The thermophysical property parameters of different formations in the target area include formation density, thermal conductivity and specific heat capacity.

6. The bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model according to claim 5 is characterized in that: Step 5 includes: according to the prior information of thermophysical parameters of different blocks in the deep formation of the target area, the upper surface average temperature boundary condition, the adiabatic boundary conditions on both sides of the deep formation heat transfer model section and the bottom prior heat flow value distribution boundary condition, taking the borehole temperature measurement as the response constraint condition, constructing the bottom boundary heat flow distribution optimization inversion objective function, such as the formula: ， in, is the fitting degree, i.e., the inverse optimization objective function; is the regression sum of squares, is the total sum of squares of deviations, is the number of measured temperature data, is the inversion temperature, that is Temperature field of the heat transfer model during the inversion process Temperature at the verification borehole; is the average measured temperature, is the residual sum of squares; To predict temperature.

7. The bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model according to claim 6 is characterized in that: Step 5 also includes: using the Newton-Raphson optimization algorithm to generate an initial random population, such as the formula: ， ， in, Indicates The individual Dimensional location; express A random number between Indicates the lower boundary; Indicates the upper boundary; is the population matrix; represents the population size, Represents the dimension of the problem.

8. The bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model according to claim 7 is characterized in that: Step 5 also includes: using the Newton-Raphson optimization algorithm to update the iterative population, such as formula: ， ， ， ， in, Newton-Raphson optimization algorithm The main components of Represents a normally distributed random number with a mean of 0 and a variance of 1; Represents a A random number for each decision variable; and is used and The positions of the two vectors generated; is the current variable value, and are distinct integers randomly selected from the population, Represents the average value of a set of data. is the process quantity, is the population matrix, The worst position, For the best position; To update The new vector position obtained is is the current vector position, and is a random number between 0 and 1. and are all random position vectors, For consideration , and The new vector position obtained is For consideration , and The new vector position obtained is is the adaptive coefficient, Indicates the current iteration, Indicates the maximum number of iterations.

9. The bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model according to claim 8 is characterized in that: Step 5 also includes: using the trap avoidance operator algorithm to further optimize the population, such as formula: ， ， ， ， in, is the position vector of the nth individual in the population matrix of the next iteration, To obtain the position vector for the next iteration using the trap avoidance operator algorithm, and They are and A uniform random number between and is a random number, For the best position, is the population matrix of the current iteration The position vector of each individual, is the adaptive coefficient, is the average value of a set of data, is the population matrix of the current iteration, For control The determining factor of performance is Represents a binary number 1 or 0. for A uniform random number between .

10. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the steps of the bottom heat flow and deep temperature inversion method based on the underground medium heat transfer model as described in any one of claims 1-9 are implemented.

Images