Geothermal resource comprehensive geophysical exploration method and system

Abstract

The invention provides a geothermal resource comprehensive geophysical exploration method, which comprises the following steps of: acquiring surface thermal anomaly data, and generating a thermal anomaly distribution diagram; collecting regional construction data, and obtaining a fault zone spatial distribution result; micro-crack dynamic signals are collected, and a real-time positioning result is obtained; high-temperature logging data in a well are collected, and a thermal reservoir porosity distribution three-dimensional image is generated. And constructing a thermal reservoir temperature field, porosity and permeability three-dimensional joint distribution model, and outputting an exploration result report. According to the method, multiple types of data are integrated, a multi-parameter continuum data field under a unified coordinate system is constructed, cross gradient constraint inversion and generative adversarial network simulation are combined, an underground medium physical parameter space coupling relation is described, and the problem that a traditional method is insufficient in resolution is solved; geological prior knowledge is embedded into a Bayesian inversion framework, inversion multiplicity is reduced, the physical rationality of an imaging result is improved, and finally, full-chain high-precision exploration capability is formed through visualization, parameter statistics and reserve calculation.

Description

Technical Field

[0001] This invention relates to the field of geothermal resource exploration technology, and in particular to a comprehensive geophysical exploration method for geothermal resources. Background Technology

[0002] Currently, geothermal resource exploration mainly relies on the combined application of single or a few geophysical methods, such as gravity exploration, magnetic exploration, seismic wave detection, and surface temperature monitoring. Traditional methods typically involve independently acquiring and processing data separately, for example, using satellite thermal infrared data to identify surface thermal anomaly areas, or using UAV magnetic / gravity measurements to infer subsurface structural features. These methods can preliminarily delineate the extent of geothermal reservoirs in specific scenarios, but they are limited by the single data dimension, simplified inversion algorithms, and lack of geological constraints, making it difficult to achieve joint inversion of multiple physical parameters and high-precision imaging. For example, conventional gravity-magnetic joint inversion often uses linear stacking models, failing to fully consider the spatial coupling relationship between fracture zones and porosity; while well logging data can provide porosity information, it lacks cross-scale correlation analysis with surface thermal anomalies and microfracture dynamics, resulting in low resolution and multiple solutions in three-dimensional imaging of geothermal reservoirs. Summary of the Invention

[0003] This invention aims to at least solve the technical problems of low resolution and multiple solutions in the existing three-dimensional imaging of geothermal reservoirs, and innovatively proposes a comprehensive geophysical exploration method for geothermal resources.

[0004] To achieve the above-mentioned objectives of this invention, this invention provides a comprehensive geophysical exploration method for geothermal resources, the method comprising: S1. Collect surface thermal anomaly data of the target area, and use satellite and calibration algorithms to invert the surface temperature based on the surface thermal anomaly data to generate a thermal anomaly distribution map; S2. The UAV magnetic / gravity exploration system collects regional tectonic data of the target area. Based on the regional tectonic data, diurnal variation correction is performed and then combined with gravity data to obtain the spatial distribution results of the fault zone. S3. Collect dynamic signals of microcracks in the target area to obtain real-time microcrack location results; S4. Collect high-temperature logging data from wells in the target area and generate a three-dimensional image of the porosity distribution of the thermal reservoir. S5. Based on the joint inversion algorithm, the spatial distribution results of the fault zone, the real-time location results of micro-fractures and the three-dimensional image of the porosity distribution of the thermal reservoir are fused. The joint distribution of the physical parameters of the underground medium is simulated by the generative adversarial network. Furthermore, geological prior knowledge is introduced to constrain the multiple solutions of the inversion and a three-dimensional joint distribution model of the temperature field, porosity and permeability of the thermal reservoir is constructed. S6. Based on the three-dimensional joint distribution model, output the exploration results report.

[0005] In another aspect, the present invention also provides a comprehensive geothermal resource exploration system, the system comprising: processor; Memory used to store processor-executable instructions; The processor is configured to implement a comprehensive geophysical exploration method for geothermal resources when executing executable instructions.

[0006] The beneficial effects of this invention are as follows: By integrating surface thermal anomalies, spatial distribution of fault zones, microfracture dynamics, and well logging data, this invention constructs a multi-parameter continuum data field under a unified coordinate system. Combined with cross-gradient constraint inversion and generative adversarial network simulation, it effectively characterizes the spatial coupling relationship of physical parameters (temperature, porosity, permeability) of the subsurface medium, solving the problem of insufficient resolution caused by isolated data dimensions in traditional methods. By embedding geological prior knowledge into a Bayesian inversion framework, it dynamically balances data fitting and prior constraints, significantly reducing the ambiguity of inversion and improving the physical rationality of imaging results. Finally, through three-dimensional visualization, parameter statistical histograms, and resource reserve calculation, it forms a high-precision exploration capability across the entire chain from data fusion to output results.

[0007] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0008] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which: Figure 1 This is a flowchart of a comprehensive geophysical exploration method for geothermal resources according to the present invention. Detailed Implementation

[0009] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0010] Example 1 like Figure 1 As shown, A comprehensive geophysical exploration method for geothermal resources, the method comprising: S1. Collect surface thermal anomaly data of the target area, and use satellite and calibration algorithms to invert the surface temperature based on the surface thermal anomaly data to generate a thermal anomaly distribution map; S2. The UAV magnetic / gravity exploration system collects regional tectonic data of the target area. Based on the regional tectonic data, diurnal variation correction is performed and then combined with gravity data to obtain the spatial distribution results of the fault zone. S3. Collect dynamic signals of microcracks in the target area to obtain real-time microcrack location results; In step S3, it is necessary to explain in detail that when collecting dynamic signals of microcracks in the target area, high-precision distributed optical fiber sensing technology is used. By deploying an optical fiber sensing network in the target area, the strain changes along the optical fiber are monitored in real time, and the dynamic propagation process of microcracks is inverted using the quantitative relationship between strain and crack aperture. In specific implementation, the optical fiber deployment scheme is first determined according to the geological structure characteristics of the target area. A spiral-encircling deployment is used in key locations such as fault zone intersections and thermal reservoir boundaries to enhance signal sensitivity, while a linear deployment is used in homogeneous rock areas to reduce costs. The raw strain signals collected by the optical fiber sensing network need to be temperature compensated and corrected to eliminate measurement errors caused by environmental temperature fluctuations. The corrected strain data is then processed by wavelet transform to denoise and extract high-frequency components related to microcrack activity. A strain-crack aperture conversion model is established based on elasticity theory, and a correction coefficient is introduced to consider the non-homogeneity of the rock. This coefficient is obtained through indoor triaxial core experiments. Finally, the existing particle filtering algorithm is used to achieve real-time positioning of microcracks with a positioning accuracy of up to the centimeter level. The positioning results are output in three-dimensional coordinate form, and key parameters such as crack direction, dip angle, and aperture change rate are marked.

[0011] In this embodiment, the method for establishing a strain-crack aperture conversion model based on elasticity theory is as follows: First, the fundamental relationship between strain and displacement in elasticity is clarified, that is, there is a specific differential relationship between strain components and displacement components. For rock media in three-dimensional space, considering the stress generated by microcrack activity, a stress-strain relationship is established according to Hooke's law. Combining the elastic modulus, Poisson's ratio, and other elastic parameters of the rock, an initial theoretical model between strain and crack aperture is derived. Then, the non-homogeneity of the rock is considered. Due to the presence of structural defects such as pores and joints inside the actual rock, its mechanical properties vary. Since spatial differences exist, a correction coefficient related to rock heterogeneity is introduced. This correction coefficient is calibrated through indoor triaxial core experiments. In the experiments, different stresses are applied to cores with different degrees of heterogeneity, and the resulting strain and simulated crack aperture changes are measured. The quantitative relationship between the correction coefficient and rock heterogeneity indices (such as porosity and joint density) is determined by data fitting. Finally, the correction coefficient is substituted into the initial theoretical model to obtain a strain-crack aperture conversion model that considers rock heterogeneity. This model can more accurately reflect the conversion relationship between strain and crack aperture in actual rocks.

[0012] S4. Collect high-temperature logging data from wells in the target area and generate a three-dimensional image of the porosity distribution of the thermal reservoir. S5. Based on the joint inversion algorithm, the spatial distribution results of the fault zone, the real-time location results of micro-fractures and the three-dimensional image of the porosity distribution of the thermal reservoir are fused. The joint distribution of the physical parameters of the underground medium is simulated by the generative adversarial network. Furthermore, geological prior knowledge is introduced to constrain the multiple solutions of the inversion and a three-dimensional joint distribution model of the temperature field, porosity and permeability of the thermal reservoir is constructed. S6. Based on the three-dimensional joint distribution model, output the exploration results report.

[0013] The principle of a comprehensive geothermal resource exploration method in this embodiment is as follows: First, multi-source heterogeneous data of the target area are collected, including surface thermal anomalies, regional structures, microfracture dynamics, and well logging data, to construct a continuum data field under a unified coordinate system. Surface thermal anomaly data is used to generate a temperature distribution map via satellite inversion; regional structural data is used to obtain fault zone distribution through magnetic / gravity joint inversion; microfracture signals are located at the centimeter level using distributed fiber optic sensing technology; and well logging data is used to generate a three-dimensional porosity image. Subsequently, a cross-gradient constraint inversion algorithm is used, with fault zone distribution as structural constraints, microfracture location as dynamic boundaries, and porosity distribution as the initial model. Multi-physics coupling inversion is achieved by minimizing the objective function. To address the problem of multiple solutions in the inversion, prior geological knowledge is introduced to construct a Bayesian inversion framework, dynamically balancing data fitting and prior constraints, ultimately generating a three-dimensional joint distribution model including temperature field, porosity, and permeability. This model integrates the distribution of thermal anomalies, the spread of fault zones, and the location of micro-fractures through three-dimensional visualization technology. It combines parameter statistical histograms with resource reserve calculations to form a complete exploration chain from data fusion to output results, which significantly improves the accuracy and reliability of geothermal resource exploration.

[0014] As an optional embodiment of the present invention, generating the thermal anomaly distribution map in step S1 may include: S101. Collect raw data of surface thermal radiation in the target area; In step S101, it should be noted that when collecting raw surface thermal radiation data of the target area, a high-precision infrared thermal imager must be used for full-coverage scanning. During the scanning process, the flight altitude must be kept constant (usually set to 50-100 meters) and the flight path spacing must not exceed 30% of the instrument's field of view to ensure that the spatial resolution of the data is better than 0.5 meters. The raw surface thermal radiation data includes multiple parameters such as surface temperature field, emissivity, and atmospheric transmittance. Atmospheric correction must be performed using the radiative transfer equation to eliminate the attenuation effect of absorption bands such as water vapor and carbon dioxide on the thermal radiation signal.

[0015] S102. Perform radiometric calibration and atmospheric correction on the raw surface thermal radiation data to eliminate response differences and the influence of atmospheric transmittance, and obtain the true surface thermal radiation brightness value. In step S102, it should be noted that when performing radiometric calibration on the raw surface thermal radiation data, a radiometric calibration table needs to be established based on the instrument parameters of the infrared thermal imager. This table describes the quantitative relationship between the digital signal output by the instrument and the actual thermal radiation intensity of the surface. Calibration experiments are conducted using a standard blackbody radiation source with known emissivity to obtain the instrument's response curves at different temperatures, thereby determining the calibration coefficients. Atmospheric correction processing requires combining real-time measured atmospheric parameters (such as water vapor content, carbon dioxide concentration, aerosol optical thickness, etc.) and using an atmospheric radiative transfer model (such as the MODTRAN or 6S model) to simulate the absorption and scattering effects of the atmosphere on thermal radiation. The atmospheric influence component is subtracted from the raw data to finally obtain the true surface thermal radiation brightness value.

[0016] S103. Based on the actual surface thermal radiation brightness value and combined with the dynamic estimation results of surface emissivity, the split window algorithm is used to retrieve surface temperature and generate continuous surface temperature spatial distribution data. In step S103, it should be noted that when retrieving surface temperature based on the true surface thermal radiation brightness value, the existing split-window algorithm constructs a temperature inversion equation by utilizing the difference in atmospheric water vapor absorption characteristics between two adjacent thermal infrared bands. Specifically, the first step is to obtain a dynamic estimation result of the surface emissivity of the target area. This result can be obtained through multispectral data inversion or ground-based measurements, and its accuracy directly affects the accuracy of the temperature inversion. Then, the true surface thermal radiation brightness value is substituted into the split-window algorithm equation, which includes parameters such as surface emissivity, atmospheric transmittance, and band-average emissivity. The equation is solved using the least squares method to obtain the initial inverted surface temperature value. Finally, the initial inversion value is verified and corrected for accuracy. This can be done by comparing and analyzing measured temperature data from ground meteorological stations, or by using cross-validation methods to assess the reliability of the inversion results. Ultimately, continuous spatial distribution data of surface temperature is generated, with a spatial resolution of up to 30 meters and a temperature accuracy better than ±0.5℃.

[0017] S104. Spatial overlay analysis is performed on continuous surface temperature spatial distribution data, elevation data, and land cover type data. Continuous surface temperature raster images are generated through interpolation processing. Surface thermal anomaly data are identified based on temperature gradient thresholds to obtain thermal anomaly areas. In step S104, it should be noted that when performing spatial overlay analysis on continuous surface temperature spatial distribution data with elevation data and land cover type data, different data sources must first be unified to the same coordinate system and resolution. Elevation data can be obtained using a Digital Elevation Model (DEM), and land cover type data can be acquired through remote sensing image classification. Using spatial analysis tools on a Geographic Information System (GIS) platform, temperature data is overlaid with elevation and land cover type data to generate a composite data layer containing multi-dimensional information. Subsequently, the existing Kriging interpolation method is used to perform spatial interpolation processing on discrete temperature measurement points. This method, by considering spatial autocorrelation, can generate continuous surface temperature raster images, with a resolution that can be set according to actual needs (typically 30-100 meters). To identify areas of surface thermal anomalies, screening based on a temperature gradient threshold is required. This threshold can be determined by statistically analyzing historical temperature data of the target area or dynamically adjusted in conjunction with geological structural background (such as fault zones and thermal reservoir distribution). In practice, the temperature gradient value of each grid cell is first calculated, which is the rate of temperature change per unit distance. Then, it is compared with a preset threshold. Areas exceeding the threshold are marked as potential thermal anomaly areas. Finally, the elevation and land cover type data are combined for further screening to exclude non-geological thermal anomalies caused by topographic relief or human activities (such as the urban heat island effect). Finally, the distribution range of surface thermal anomalies associated with underground thermal reservoirs is obtained.

[0018] S105. Based on the geological structural background, classify and visualize the thermal anomaly areas to generate a thermal anomaly distribution map that conforms to geological exploration specifications.

[0019] In step S105, it should be noted that when classifying thermal anomaly regions based on geological structural background, a comprehensive analysis of geological factors such as fault zone distribution, reservoir boundaries, and lithological variations is necessary. First, based on the spatial relationship between the thermal anomaly region and the fault zone, thermal anomalies located at fault zone intersections or along fault zone extension directions are designated as Level 1 regions. These regions are typically closely related to the upward flow channels of deep heat sources. Second, thermal anomalies distributed on both sides of a single fault zone or near reservoir boundaries are designated as Level 2 regions, whose heat conduction processes are significantly controlled by geological structures. Finally, thermal anomalies far from fault zones and distributed in homogeneous rock mass areas are designated as Level 3 regions, primarily reflecting changes in the shallow geothermal field. The classification criteria need to be dynamically adjusted based on historical exploration data of the target area. For example, in a known geothermal field, the temperature gradient threshold for Level 1 regions can be set to ≥3℃ / 100m, for Level 2 regions to 1.5-3℃ / 100m, and for Level 3 regions to 0.5-1.5℃ / 100m. The visualization mapping process employs a layered color scheme. First-level areas use a red gradient to represent high-temperature anomalies, second-level areas use an orange gradient to represent medium-temperature anomalies, and third-level areas use a yellow gradient to represent low-temperature anomalies. Auxiliary information such as fault zone strikes, thermal reservoir boundaries, and temperature gradient contour lines are also labeled on the map. The final generated thermal anomaly distribution map includes map title, scale, directional indicators, legend, and responsibility label. The spacing between temperature contour lines is set to 0.5-2℃ based on data accuracy, and the contour line values ​​are labeled. This distribution map organically integrates geological structural constraints with thermal anomaly characteristics.

[0020] As an optional embodiment of the present invention, optionally, obtaining the spatial distribution result of the fault zone in step S2 includes: S201. Use an unmanned aerial vehicle (UAV) magnetic exploration system to collect magnetic gradient data of the target area and obtain raw magnetic anomaly data; In step S201, it should be noted that when using an unmanned aerial vehicle (UAV) magnetic exploration system to collect magnetic gradient data of the target area, a high-precision magnetometer is required. This magnetometer must have low noise and high sensitivity characteristics, capable of capturing weak magnetic anomaly signals. The UAV must maintain a stable attitude during flight, with the flight altitude set between 50-150 meters. The flight path spacing is determined based on the magnetometer's field of view, typically not exceeding 70% of its effective detection width, to ensure complete data coverage and a spatial resolution better than 10 meters. The collected raw magnetic anomaly data includes multiple parameters such as the total geomagnetic field strength and gradient changes; auxiliary information such as flight trajectory and timestamps must be recorded.

[0021] S202. Perform diurnal variation correction on the original magnetic anomaly data to eliminate the influence of diurnal variation of the geomagnetic field and obtain the corrected magnetic anomaly data. In step S202, it should be noted that when performing diurnal variation correction on the original magnetic anomaly data, a geomagnetic diurnal variation observation station needs to be established near the target area, and a high-precision geomagnetic instrument should be used to continuously record the variation curve of the geomagnetic field over time. This curve reflects the diurnal variation law of the geomagnetic field, including short-term fluctuations such as solar diurnal variation and geomagnetic storms. During data processing, the original magnetic anomaly data is first segmented into time series according to the UAV flight trajectory. Then, the magnetic field variation values ​​for the corresponding time period are obtained from the diurnal variation observation station. The diurnal variation data is matched to each measurement point using a linear interpolation method. Finally, the diurnal variation component is subtracted from the original magnetic anomaly to obtain corrected magnetic anomaly data that only reflects the magnetic differences of underground geological bodies. This processing process needs to consider the spatial distance between the UAV and the diurnal variation station. When the distance exceeds 50 kilometers, a geomagnetic field gradient correction model needs to be introduced to eliminate the influence of spatial differences and ensure that the diurnal variation correction accuracy is better than ±1 nT. The corrected magnetic anomaly data should retain the spatial distribution characteristics of the original data while eliminating false anomalies caused by solar radiation, ionospheric disturbances, etc.

[0022] S203. Simultaneously use the UAV gravity exploration system to collect gravity anomaly data of the target area, and perform terrain correction and Bouguer correction to obtain gravity gradient data; In step S203, it should be noted that when simultaneously using the UAV gravity survey system to collect gravity anomaly data of the target area, a high-precision gravimeter must be installed. This gravimeter must possess high stability and low drift characteristics, and be able to accurately measure minute changes in the gravity field. The UAV must also maintain a stable attitude during flight, with the flight altitude set between 100-300 meters. The flight path spacing is determined based on the gravimeter's detection range, typically not exceeding 60% of its effective detection width, to ensure data continuity and a spatial resolution better than 50 meters. The collected gravity anomaly data includes multiple parameters such as gravity field intensity and gradient changes, and auxiliary information such as flight trajectory and altitude must be recorded simultaneously. When performing terrain correction, digital elevation model (DEM) data must be used to eliminate the influence of terrain undulations on gravity measurements. This process involves calculating the difference between the terrain elevation at each measurement point and the gravimeter's observation height, and then using a terrain correction formula to subtract the gravity anomaly component caused by the terrain. Bouguer correction further considers factors such as Earth's curvature and altitude correction, using the Bouguer correction formula to normalize gravity data to the same reference surface, eliminating gravity variations caused by altitude differences. After topographic correction and Bouguer correction, the obtained gravity gradient data can more accurately reflect the density differences of underground geological bodies.

[0023] S204. The corrected magnetic anomaly data and gravity gradient data are jointly inverted. The least squares inversion algorithm is used to fuse the magnetic and gravity data to generate a model of underground density and magnetic susceptibility distribution. The expression for the least squares inversion algorithm is: ; , ; in, This represents the least squares objective function, which comprehensively measures the goodness of fit between magnetic and gravity data and the smoothness of the model. This indicates the total number of observation points collected by the UAV magnetic exploration system. Indicates the first Corrected magnetic anomaly data for each observation point Indicates the first Forward modeling values ​​of magnetic anomalies at each observation point This indicates the total number of observation points collected by the UAV gravity exploration system. Indicates the first Gravity gradient data at each observation point Indicates the first The forward modeling values ​​of the gravity gradient at each observation point This represents the regularization coefficient, used to balance the data fitting term and the model smoothing term. Represents the spatial gradient of the density model. This represents the spatial gradient of the magnetic susceptibility model. This represents a model of underground magnetic susceptibility distribution. Indicates model parameters (density) in the underground medium ,magnetic susceptibility The spatial coordinates of the location, Represented as the magnetic dipole field kernel function, Represents the spatial coordinates of a surface observation point. Represents the underground density distribution model. This represents the kernel function of the Newtonian gravitational gradient.

[0024] In step S204, it should be noted that when jointly inverting the corrected magnetic anomaly data and gravity gradient data, the least squares inversion algorithm achieves coordinated constraints between magnetic and gravity data by constructing an objective function that includes a data fitting term and a model smoothing term. The data fitting term requires minimizing the sum of squared residuals between the observed data (corrected magnetic anomaly and gravity gradient) and the forward-modeled values ​​(magnetic anomaly and gravity gradient simulated based on the current model parameters), ensuring that the model can interpret the actual observation results. The model smoothing term, by constraining the spatial gradients of the density model and the magnetic susceptibility model, suppresses numerical oscillations during the inversion process and generates a continuous distribution model that conforms to geological laws. The regularization coefficient λ is a key parameter, and its value needs to be determined through cross-validation or the L-curve method. When λ is too small, the model overfits the observed noise, leading to the amplification of local anomalies; when λ is too large, the model is over-smoothed, losing true geological information. The magnetic dipole field kernel function and the Newtonian gravitational gradient kernel function describe the field effects generated by the subsurface magnetic body and density body on the surface, respectively. Their integral forms ensure comprehensive coverage of the subsurface three-dimensional space during the inversion process. By iteratively optimizing the objective function, the resulting subsurface density and magnetic susceptibility distribution models can simultaneously reflect the density differences and magnetic characteristics of geological bodies, providing multi-physics constraints for the accurate identification of geological structures such as fault zone distribution and thermal reservoir boundaries. For example, in a known fault zone region, the density model may exhibit a linear gradient change, while the magnetic susceptibility model shows high-value anomalies due to the enrichment of magnetic minerals within the fault zone. The degree of agreement between the two spatial locations can effectively verify the reliability of the inversion results.

[0025] S205. Based on the underground density and magnetic susceptibility distribution model, linear anomaly zones are identified through edge detection algorithms. Combined with geological structural feature analysis, a three-dimensional image of the spatial distribution of fault zones is generated, and the spatial distribution results of fault zones are output.

[0026] The expression for the edge detection algorithm is: in, This represents the subsurface density / magnetic susceptibility model at a spatial point. gradient strength at that point The density / magnetic susceptibility model is represented in The value at that location, and This represents the kernel function of the Sobel operator. Represents a three-dimensional convolution operation; In step S205, it should be noted that when performing edge detection based on the underground density and magnetic susceptibility distribution model, the 3D model data needs to be divided into regular grid cells, with the density and magnetic susceptibility values ​​of each cell serving as input parameters. The gradient components in the x, y, and z directions are calculated using the Sobel operator kernel function. The horizontal direction (xy plane) uses a 2D Sobel operator extended to 3D space, while a dedicated kernel function is designed for the vertical direction (z-axis) to capture vertical changes. The 3D convolution operation traverses the entire model space using a sliding window, calculating the gradient intensity value at each grid point. This value comprehensively reflects the rate of change of density / magnetic susceptibility in 3D space. To highlight linear anomaly features, non-maximum suppression processing is applied to the gradient intensity image, retaining only local maxima along the gradient direction, while a gradient threshold is set to filter out minor fluctuations caused by noise. In the gradient image after thresholding, continuous high-value regions represent potential linear anomaly zones, whose orientation closely matches the extension direction of the fault zone. When analyzing geological structural features, the identified linear anomaly zones need to be compared and verified with regional geological maps, seismic profiles, and other data. The spatial relationship between the anomaly zones and known fault zones and thermal reservoir boundaries should be analyzed. For example, when a linear anomaly zone crosses contact zones of different lithologies or overlaps with areas of low gravity or high magnetic strength, it can be further confirmed as a manifestation of a fault zone. If the anomaly zone intersects regional tectonic lines at a small angle or extends parallel to them, it is necessary to consider whether it is a secondary fault or associated structure. The final generated 3D image of the fault zone's spatial distribution uses a transparency overlay display technique to simultaneously display the density model, magnetic susceptibility model, and edge detection results. Density anomaly areas are represented by blue tones, magnetic susceptibility anomaly areas by red tones, and linear anomaly zones are outlined with yellow highlighted lines. The image output should include auxiliary information such as 3D coordinate axes, scale bar, and legend, and support viewing at any angle of rotation to intuitively present the 3D morphology of the fault zone in underground space. The result can be exported to a common 3D model format (such as .obj or .stl).

[0027] As an optional embodiment of the present invention, optionally, the process of acquiring high-temperature logging data in the target area in step S4 to generate a three-dimensional image of the porosity distribution of the thermal reservoir includes: S401. Synchronously acquire high-temperature logging data in the well. The high-temperature logging data includes wellbore axial temperature profile data and sonic transit time data, and record well depth positioning information. In step S401, it should be noted that when synchronously acquiring high-temperature logging data from the well, a logging instrument combination with high-temperature resistance characteristics must be used. The temperature sensor must be a platinum resistance thermometer or thermocouple type, with a range covering the expected highest temperature in the target area and a 20% safety margin to ensure measurement accuracy of ±0.5℃ even in high-temperature environments. The sonic logging tool must be equipped with a high-frequency transducer (center frequency ≥20kHz) to enhance the response sensitivity to pore fluids. Well depth positioning information is recorded in real time using a high-precision depth encoder with a resolution better than 0.01 meters, and stored synchronously with the logging data. During acquisition, the logging instrument must be lowered / raised at a uniform speed (controlled at 3-5 meters / minute) to avoid data distortion due to speed fluctuations. Simultaneously, the fluid circulation status within the wellbore must be monitored. Acquisition should be paused when mud discharge fluctuations exceed ±10%, and operations should resume after the fluid stabilizes. The wellbore axial temperature profile data should include continuous temperature measurements from the wellhead to the bottom, with the sampling interval dynamically adjusted according to temperature gradient changes (typically densified to 0.1 meters in areas of rapid temperature variation). Acoustic transit time data should record the arrival time of the first wave at each depth point for subsequent porosity calculations. Well depth positioning information must strictly correspond to the logging data, synchronized via timestamps or depth pulse signals, ensuring spatial matching accuracy of temperature, acoustic, and depth parameters better than 0.05 meters in subsequent processing.

[0028] S402. Perform wellbore compensation and formation invasion correction on the acoustic time difference data, and calculate the in-situ acoustic velocity of the thermal reservoir by combining the temperature profile data. In step S402, it should be noted that when performing wellbore compensation correction on acoustic transit time data, the impact of wellbore diameter changes on the acoustic propagation path must be considered. Increasing the wellbore diameter leads to an increase in the acoustic propagation distance, resulting in a larger measured acoustic transit time; conversely, decreasing the wellbore diameter reduces the acoustic propagation distance, resulting in a smaller measured acoustic transit time. By establishing a correction model between wellbore diameter and acoustic transit time, the actual wellbore diameter data is used to correct the acoustic transit time, eliminating errors caused by wellbore diameter changes. Formation invasion correction addresses the phenomenon of drilling fluid intrusion into the formation during drilling. Drilling fluid intrusion alters the properties and distribution of formation pore fluids, thus affecting the propagation speed of acoustic waves within them. By analyzing the depth and extent of drilling fluid intrusion into the formation, and combining parameters such as formation porosity and permeability, a formation invasion correction model is established to correct the acoustic transit time data, obtaining a more accurate in-situ acoustic velocity. When calculating the in-situ acoustic velocity of thermal reservoirs using temperature profile data, the effect of temperature on acoustic propagation speed must be considered. Increased temperature accelerates the propagation speed of acoustic waves in the medium. Based on the experimentally obtained temperature-sound wave velocity relationship curve, the corrected sound wave transit time data and corresponding temperature data are substituted into the formula to calculate the in-situ sound wave velocity of the thermal reservoir (the calculation method is existing technology and will not be elaborated here). This velocity can more realistically reflect the physical properties of the thermal reservoir formation.

[0029] S403. Based on the in-situ acoustic velocity of the thermal reservoir, an acoustic velocity-porosity conversion function is established, and the discrete point acoustic velocity data is extrapolated to three-dimensional space using the Kriging interpolation algorithm. The expression for the acoustic velocity-porosity conversion function is: ; in, Indicates underground point Porosity at that location The velocity of sound waves in the rock matrix is ​​obtained through core experiments and serves as the reference parameter for the transformation function. Indicates the in-situ acoustic velocity of the thermal reservoir. The velocity of sound waves in pore fluid is represented by well temperature gradient data. Dynamic adjustment This represents the temperature sensitivity coefficient, determined by fitting measured temperature-velocity data from the well, and is used to quantify the effect of temperature on sound wave velocity. This represents the axial temperature gradient of the wellbore.

[0030] In step S403, it should be noted that the acoustic velocity of the rock matrix is ​​obtained through core experiment calibration. This process requires selecting representative core samples and measuring the acoustic velocity under strictly controlled experimental conditions to obtain accurate acoustic velocity values ​​of the rock matrix, which serve as the reference parameters for the conversion function. The in-situ acoustic velocity of the thermal reservoir is obtained after correction and calculation in the previous steps, and it truly reflects the physical properties of the thermal reservoir formation. The pore fluid acoustic velocity is dynamically adjusted based on the well temperature gradient data. Since temperature affects the pore fluid acoustic velocity, the pore fluid acoustic velocity is reasonably adjusted according to its changes by monitoring the well temperature gradient data to make it more consistent with the actual situation. The temperature sensitivity coefficient is determined by fitting the measured temperature-velocity data in the well. Temperature and corresponding acoustic velocity data at different locations in the well are collected, and the temperature sensitivity coefficient is obtained using mathematical fitting methods to quantify the influence of temperature on the acoustic velocity. The axial temperature gradient of the wellbore can be measured by a temperature sensor.

[0031] When extrapolating discrete-point acoustic velocity data to three-dimensional space using the Kriging interpolation algorithm, a geostatistical method based on spatial autocorrelation, the algorithm considers the spatial distribution characteristics of known discrete-point acoustic velocity data and the correlation between them, describing this spatial autocorrelation by constructing a semivariance function. Based on the semivariance function model, acoustic velocities in unknown areas are estimated and interpolated, extending the discrete acoustic velocity data to the entire three-dimensional space, thus obtaining a continuous acoustic velocity distribution of the thermal reservoir. This method fully utilizes the information from known data and considers spatial correlation, making the interpolation results more consistent with actual geological conditions.

[0032] S404. Based on the acoustic velocity-porosity conversion function, the temperature gradient abrupt point and the acoustic velocity anomaly region are fused, and the measurement error interference is eliminated by Bayesian joint inversion to generate a porosity probability distribution model. In step S404, it should be noted that when integrating temperature gradient abrupt change points and acoustic velocity anomaly zones, gradient abrupt change points in the temperature profile data must first be identified using numerical analysis methods. These abrupt change points typically correspond to stratigraphic lithological change interfaces or fluid migration channels, and their spatial location can be determined by calculating the first or second derivative of the temperature data. When the derivative value exceeds a preset threshold, it is determined to be an abrupt change point. Simultaneously, the identification of acoustic velocity anomaly zones is based on a three-dimensional acoustic velocity distribution model. By setting a velocity threshold range or velocity gradient change rate, anomalies that differ significantly from the surrounding area are screened out. These anomaly zones may reflect abnormal porosity development or changes in fluid properties. In the Bayesian joint inversion, the prior probability distribution is set based on geological exploration experience (e.g., porosity usually follows a log-normal distribution), while the likelihood function is constructed by describing the error relationship between the observed data and model parameters. During the inversion process, the existing Markov Chain Monte Carlo (MCMC) algorithm is used to sample the posterior probability distribution, and iterative optimization is used to gradually approximate the true porosity distribution. This algorithm effectively suppresses the ambiguity of inversion from a single data source by introducing spatial constraints between abrupt temperature gradient changes and anomalous acoustic velocity regions. Simultaneously, it quantifies porosity uncertainty through a probability distribution. For example, when abrupt temperature gradient changes and low acoustic velocity regions spatially overlap, the inversion result assigns a higher porosity probability value to that region. Conversely, if the spatial deviation between the two is significant, the joint probability weight is reduced to mitigate anomalous interference. The final porosity probability distribution model is presented as three-dimensional volumetric data, with each grid point storing the expected porosity value and its standard deviation. A smaller standard deviation indicates a more reliable porosity estimate at that location.

[0033] S405. Construct a radial basis function grid centered on the borehole, convert the porosity probability distribution model into three-dimensional volume data of porosity isosurface, and output a three-dimensional image of the porosity distribution of the thermal reservoir.

[0034] In step S405, it should be noted that when constructing the radial basis function grid centered on the borehole, the selection of the radial basis function must balance computational efficiency and fitting accuracy. Commonly used are existing Gaussian radial basis functions (with the central decay rate determined through cross-validation optimization) or existing thin-plate spline radial basis functions (suitable for smooth, continuous distributions). The grid nodes are densified exponentially along the borehole radial direction (near-well spacing ≤ 0.5 meters, gradually widening to 2 meters further away), and the vertical layer thickness is dynamically adjusted according to the thermal reservoir thickness (usually each layer thickness ≤ 1 meter). The expected value of the porosity probability distribution model is used as input data, and the porosity value of each grid node is calculated through radial basis function interpolation to generate a continuous three-dimensional porosity field. Boundary constraints must be considered during the conversion process: when grid points exceed the study area, background porosity values ​​obtained from geological exploration (usually 5%-15%) are extrapolated and filled; for areas with missing data, spatial continuity is ensured through joint modeling using Kriging interpolation and radial basis functions. The final generated 3D porosity isosurface data employs a layered color-coding technique, dividing porosity intervals into 5-10 color levels (e.g., 5%-10% blue, 10%-15% green, 15%-20% yellow, 20%-25% orange, and >25% red). The isosurface rendering accuracy is better than 0.1 meters. When outputting the image, borehole trajectory lines (solid red lines), formation boundaries (dashed black lines), and temperature gradient isosurfaces (dashed gray lines) are overlaid. A porosity statistical histogram (showing the frequency of porosity distribution across the entire area) and confidence intervals (porosity fluctuation range at a 95% confidence level) are also labeled.

[0035] As an optional embodiment of the present invention, optionally, constructing a three-dimensional joint distribution model of the thermal reservoir temperature field, porosity, and permeability in step S5 includes: S501. Based on the spatial distribution results of fracture zones, real-time location results of micro-fractures, and three-dimensional images of thermal reservoir porosity distribution, a three-dimensional mesh alignment algorithm is used to unify the spatial coordinate system. Discrete point data is converted into a continuous data field through nearest neighbor interpolation to form a multi-source data fusion preprocessing result. The expression for the 3D mesh alignment algorithm is: , ; in, The target mesh point set is generated using a 3D mesh alignment algorithm. This represents a coordinate transformation operator constrained by geological structure. Represents the original discrete point dataset. Represents grid points The k-th attribute value (such as density, porosity) at the location. ∈{ , , The corresponding thermal reservoir temperature field, porosity, and permeability express The number of discrete data points included. This represents the attribute weight coefficient. Indicates the first The attribute values ​​of each data point This represents the Diracdelta function, used to select the nearest neighbor discrete point. express The first in Each specific data point is a nearest neighbor interpolation function. Input variables; In step S501, it is necessary to explain in detail that the 3D mesh alignment algorithm aims to establish a unified spatial coordinate system, mapping discrete point data from different sources and with different resolutions to the same 3D mesh framework. The geological structure-constrained coordinate transformation operator integrates geological structural information such as the fault zone distribution direction and micro-fracture extension trends, performing rotation, scaling, and translation transformations on the original coordinate system to ensure that all data achieves spatial alignment under geological structural constraints. For example, when the fault zone extends in a northeast direction, the coordinate transformation operator will preferentially stretch the data along the fault zone direction to maintain the continuity of the geological structure. The nearest neighbor interpolation function calculates the Euclidean distance between the target mesh point and the original discrete points, selecting the nearest discrete point as the interpolation reference. The setting of attribute weight coefficients needs to consider the reliability of the data source. For example, for micro-fracture location data obtained from high-precision microseismic monitoring, a higher weight coefficient (≥0.8) can be assigned; while for porosity data obtained through well logging data inversion, the weight should be appropriately reduced (∈[0.5,0.7]). The Diracdelta function takes a value of 1 when a discrete point coincides with a target grid point, and 0 otherwise, ensuring that each target grid point's attribute value is contributed only by the nearest discrete point. After processing by this algorithm, the originally scattered spatial distribution data of fracture zones, real-time location data of micro-fractures, and porosity distribution data of thermal reservoirs are unified into a three-dimensional grid of the same resolution (grid size is typically set to 10m × 10m × 5m depending on the study area), forming a continuum data field. This process not only eliminates spatial misalignment between multi-source data but also balances the accuracy differences between different data sources through attribute weighting coefficients.

[0036] S502. Based on the preprocessing results of multi-source data fusion, a joint inversion objective function is constructed using a cross-gradient constraint inversion algorithm. The spatial distribution of the fault zone is used as the construction constraint, the micro-fracture location results are used as the dynamic boundary conditions, and the porosity distribution is used as the initial model. By minimizing the objective function, a multi-physics joint inversion is performed to generate an initial distribution model of the physical parameters of the underground medium. The expression for the cross-gradient constraint inversion algorithm is: ; in, Represents the gradient vector of physical parameters of the underground medium. ={ , , , } Corresponding to density, magnetic susceptibility, porosity, and permeability, Represents the gradient vector constrained by geological structure. The vector cross product operator represents the physical essence of cross gradient constraints. Through the spatial coupling of gradients of different parameters, it forces the spatial correlation between fracture zones and porosity anomaly zones, and between microcracks and permeability mutation zones. The expression for the joint inversion objective function is: in, This represents the joint inversion objective function, which comprehensively measures the goodness of fit of multiphysics data, the strength of cross-gradient constraints, and the degree of geological prior matching. , , The regularization coefficients represent the weights that balance data fitting, cross-gradient constraints, and prior constraints, and are dynamically adjusted using a Bayesian inversion algorithm. Indicates the first Observational data ( =1 represents magnetic data. =2 represents gravity data. =3 represents porosity data). Represents the forward-calculated value, obtained through a kernel function. (Magnetic dipole field) and (Newton's gravitational gradient) calculation, Represents the gradient vector of physical parameters of the underground medium. Represents the gradient vector constrained by geological structure. This represents a priori geological model.

[0037] In step S502, it is necessary to explain in detail that the cross-gradient constraint inversion algorithm forces the establishment of a correlation between geological structures and physical parameters through the spatial coupling between gradients of different physical parameters. Specifically, the geological structure constraint gradient vector is generated by integrating geological information such as the distribution direction of fault zones and the extension trend of microfractures, and its component direction is consistent with the direction of the main structural line; while the gradient vector of physical parameters of the subsurface medium is composed of the spatial variation rate of parameters such as density, magnetic susceptibility, porosity, and permeability. When the vector cross product operator is applied to both, it only produces a non-zero value when there is a significant deviation between the parameter gradient direction and the structural gradient direction, thereby forcing the parameter anomaly zone and the structural zone to spatially coincide through the process of minimizing the objective function. For example, when the angle between the porosity gradient direction and the fault zone structural gradient direction exceeds 30°, the cross gradient term will significantly increase the objective function value, forcing the inversion result to adjust the porosity distribution until the two directions tend to be consistent.

[0038] The joint inversion objective function consists of three parts: a data fitting term, which quantifies the degree of matching between the model and measured data by calculating the L2 norm difference between observed data and forward simulation values; a cross-gradient constraint term, which measures the coupling strength between parameter gradients and tectonic gradients through the aforementioned vector cross product; and a geological prior constraint term, which ensures that the inversion results conform to regional geological patterns by calculating the Mahalanobis distance between the current model and the background model obtained from geological exploration. The dynamic adjustment mechanism of the regularization coefficient is implemented based on a Bayesian framework: in the initial stage, the data fitting term is given a high weight to quickly approximate the true solution; as the number of iterations increases, the weight of the cross-gradient constraint is gradually increased; and in the final stage, the prior constraint is strengthened to eliminate abnormal fluctuations caused by non-geological factors. Forward calculation values ​​are obtained through kernel function convolution: magnetic data uses the magnetic dipole field kernel function, gravity data uses the Newton gravitational gradient kernel function, and porosity data is generated through radial basis function interpolation kernels.

[0039] The algorithm iteratively optimizes the model parameters: each iteration first calculates the forward response of the current model, then adjusts the parameter values ​​based on the gradient descent direction of the objective function, until the convergence condition is met (the rate of change of the objective function value <1% or the number of iterations reaches 200). The final generated initial distribution model of underground medium physical parameters contains four parameter fields: density field (resolution 10m×10m×5m), magnetic susceptibility field (resolution 15m×15m×8m), porosity field (resolution 5m×5m×2m), and permeability field (resolution 8m×8m×4m). The parameter fields are spatially aligned using a 3D mesh alignment algorithm. Cross-validation is used in the model validation phase: the dataset is randomly divided into a training set (70%) and a validation set (30%). After obtaining the model parameters through inversion from the training set, the data fitting residuals (root mean square error <5%) and cross-gradient coupling (average angle <15°) are calculated on the validation set. The final model is output only after both indicators meet the criteria.

[0040] S503. Based on the initial distribution model of underground medium physical parameters, the joint distribution of underground medium physical parameters is simulated using a generative adversarial network. The generator takes the initial model as input to generate a simulated parameter set, and the discriminator evaluates the matching degree between the parameter set and the real data based on actual geological data. The inversion parameter set is optimized through iterative training until the matching degree reaches the standard, and the optimized inversion parameter set is obtained. In step S503, it is necessary to explain in detail that the generative adversarial network consists of two parts: a generator and a discriminator. These two parts continuously optimize parameters through an adversarial training mechanism, ultimately achieving a high degree of matching between the simulated parameter set and real geological data. The generator adopts a three-dimensional convolutional neural network architecture. Its input is an initial distribution model of underground medium physical parameters (including four parameter fields: density, magnetic susceptibility, porosity, and permeability). It extracts spatial features through multiple convolutional kernels (3×3×3 size, stride 1) and uses deconvolutional layers (stride 2) to achieve upsampling, finally outputting a simulated parameter set with the same resolution as the real geological data. To enhance the model's generalization ability, the generator introduces a residual connection mechanism, adding low-level features to high-level features to avoid the gradient vanishing problem. Simultaneously, a batch normalization layer (BatchNorm) is used to accelerate training convergence, and a LeakyReLU activation function (slope 0.2) is used instead of the traditional ReLU to alleviate neuron death.

[0041] The discriminator employs a dual-branch structure, processing simulated parameter sets and actual geological data separately. Each branch contains four convolutional layers (kernel size 5×5×5, stride 2) and two fully connected layers. Max pooling layers (size 2×2) are used to reduce feature dimensionality, ultimately outputting a matching score within the range of 0-1. To improve discrimination accuracy, the discriminator introduces an attention mechanism: a channel attention module is added after the third convolutional layer, generating channel weights through global average pooling to enhance the feature representation of key parameters (such as porosity anomalies); simultaneously, a spatial attention module is used to generate a spatial weight map through convolution operations, focusing on geologically sensitive areas (such as near fault zones).

[0042] The training process employs an alternating optimization strategy: when the generator parameters are fixed, the discriminator learns to distinguish between real and simulated data by minimizing the binary cross-entropy loss function (BCEWithLogitsLoss); when the discriminator parameters are fixed, the generator optimizes the parameter set by maximizing the discriminator output, making the simulated data closer to the real distribution. To stabilize the training process, the Wasserstein distance loss function is used instead of the traditional BCE loss, and the gradient norm of the discriminator (Lipschitz continuity) is constrained by the gradient penalty term to avoid mode collapse. A low learning rate is set in the early stage of training, and it decays exponentially (decay rate 0.95) after every 100 iterations, with a total of 2000 iterations. Training is terminated early when the discriminator score fluctuation range is less than 0.05 and the generator loss value decreases by less than 1% for 50 consecutive iterations.

[0043] The optimized inversion parameter set must meet three quality indicators: For data fit, the root mean square error (RMSE) between the simulated parameter set and the actual geological data must be less than 3%; for spatial correlation, by calculating the cross-covariance matrix of the simulated and real parameter fields, the main diagonal elements (autocorrelation coefficients) must be greater than 0.85, and the off-diagonal elements (cross-correlation coefficients) must be less than 0.3; for geological rationality, 100 parameter sets are generated using the Markov Chain Monte Carlo (MCMC) method, and their means and standard deviations are statistically analyzed. The 95% confidence interval of the optimized parameter set must have an overlap rate of more than 80% with the background model obtained from geological exploration. The final output inversion parameter set includes four optimized parameter fields: density field, magnetic susceptibility field, porosity field, and permeability field. Each parameter field is spatially aligned using a three-dimensional grid alignment algorithm, and the porosity field and permeability field satisfy a power-law relationship between permeability and porosity.

[0044] S504. Based on the optimized inversion parameter set, geological prior knowledge is introduced to constrain the inversion ambiguity. The geological prior probability distribution is determined by combining the thermal anomaly distribution map and the fault zone distribution results. The data fitting and prior constraints are balanced by the existing Bayesian inversion algorithm to eliminate the inversion ambiguity and obtain the eliminated inversion parameter set. In step S504, it should be noted that the geological prior knowledge constraint utilizes known geological information to narrow the inversion solution space and address the problem of multiple solutions in the inversion results of a single geophysical method. Specifically, a geological prior probability distribution model needs to be constructed first: based on the thermal anomaly distribution map, geothermal gradient anomaly areas are extracted; combined with the fault zone distribution results, tectonic activity zones are divided (extending 500 meters on each side of the fault zone); the study area is divided into a high-probability zone (satisfying both thermal anomaly and fault zone conditions), a medium-probability zone (satisfying only one condition), and a low-probability zone (satisfying neither condition), assigned prior probability values ​​of 0.8, 0.5, and 0.2, respectively. This probability distribution is extended to the entire study area using a three-dimensional kriging interpolation algorithm to generate a prior probability field with the same resolution as the inversion parameter set.

[0045] Existing Bayesian inversion algorithms construct a posterior probability density function by fusing the likelihood function of observed data with the geological prior probability distribution: P(m|d)∝P(d|m)·P(m), where m represents the model parameter vector, d represents the observed data vector, P(d|m) is the data likelihood term (measured by the L2 norm as the residual between the observed data and the forward simulation value), and P(m) is the geological prior term.

[0046] The Metropolis-Hastings algorithm is used for Markov chain Monte Carlo sampling during the inversion process: the initial model is randomly selected from the optimized inversion parameter set. In each iteration, when generating candidate models, a Gaussian perturbation (standard deviation of 5% of the parameter mean) is added to the current model parameters. The acceptance probability α = min{1, [P(d|m)·P(m)] / [P(d|m)·P(m)]} is calculated. If the random number < α, the candidate model is accepted; otherwise, the current model is retained. After 5000 iterations (the first 1000 iterations are discarded as a pre-burning period), the model with the highest posterior probability density is selected as the final solution.

[0047] This algorithm eliminates inversion ambiguity through a triple mechanism: data fitting terms force the model to match the observed data, geological prior terms constrain the model to conform to regional geological patterns, and probability weights under a Bayesian framework balance the contributions of both. For example, when the inversion results of physical property parameters in a certain region conflict with the geological prior probability (such as low porosity values ​​in high-probability areas), the posterior probability density will decrease significantly, forcing subsequent iterations to adjust towards the high-probability areas. The final output inversion parameter set after eliminating ambiguity shows that the spatial distribution of the porosity and permeability fields matches the thermal anomaly center and fault zone trend by more than 40%, and the probability density functions of each parameter field all exhibit a unimodal characteristic, indicating that ambiguity is effectively suppressed.

[0048] S505. Based on the eliminated inversion parameter set, a three-dimensional joint distribution model of the temperature field, porosity, and permeability of the thermal reservoir is constructed.

[0049] In step S505, it should be noted that constructing a three-dimensional joint distribution model involves integrating the inversion parameter set after eliminating multiple solutions, and using multi-parameter collaborative constraints to characterize the features of the thermal reservoir. The implementation is divided into three stages: First, based on density and magnetic susceptibility data, a basic geological framework is generated using a three-dimensional kriging interpolation algorithm. This framework determines spatial correlation parameters (nuclear effect value < 0.1, sill value > 0.9) through variogram analysis to ensure the continuous expression of geological structures (such as anticlines and synclines). Second, porosity and permeability fields are used as key physical property parameters, and a power-law relationship between permeability and porosity is introduced for joint constraint. Regional parameters are fitted using nonlinear least squares method to ensure that the spatial distribution of permeability and porosity fields conforms to geological laws (such as high porosity areas corresponding to high permeability areas). Finally, combined with geothermal gradient data in the thermal anomaly distribution map, a radial basis function interpolation method is used to generate a temperature field, which is dynamically corrected using the heat conduction equation (∂T / ∂t = α∇²T, where α is the thermal diffusivity coefficient) to ensure the coupling relationship between the temperature field and the fault zone distribution (obtained through seismic profile interpretation). The construction of the 3D joint distribution model adopts a layered modeling strategy: the bottom layer is the geological structure layer, obtained through inversion of density and magnetic susceptibility data; the middle layer is the physical property parameter layer, integrating porosity and permeability fields; and the top layer is the temperature field layer, generated based on thermal anomaly data. Spatial matching between layers is achieved through a 3D mesh alignment algorithm, and model accuracy is verified through cross-validation: the dataset is divided into a training set (60%), a validation set (20%), and a test set (20%). On the validation set, the root mean square error (RMSE < 4%) and structural similarity index (SSIM > 0.85) are calculated, and on the test set, the model's predictive ability for unknown areas is evaluated (prediction error < 5%). The final generated 3D joint distribution model contains four core parameter fields: temperature field, porosity field, permeability field, and geological structure field (fault zone identification accuracy > 90%). Each parameter field is displayed at multiple scales through 3D visualization techniques (such as volume rendering and isosurface extraction), supporting thermal reservoir resource assessment (e.g., thermal reservoir volume calculation error < 8%) and development scheme design (e.g., well location deployment optimization).

[0050] As an optional embodiment of the present invention, optionally, in step S6, the exploration results report output based on the three-dimensional joint distribution model includes: S601. Based on the three-dimensional joint distribution model, multi-source data fusion visualization is performed, and thermal anomaly distribution map, fracture zone spatial distribution results and microcrack real-time location results are superimposed to generate a three-dimensional visualization result image. In step S601, it is necessary to explain in detail that multi-source data fusion visualization achieves a three-dimensional display of the characteristics of thermal reservoirs by integrating geological information from different sources and at different scales. In specific implementation, firstly, the temperature field, porosity field, permeability field and geological structure field in the three-dimensional joint distribution model are unified in data format, and a three-dimensional grid alignment algorithm is used to ensure that each parameter field is accurately matched in the spatial coordinate system; secondly, the thermal anomaly distribution map (generated through geothermal gradient data) and the spatial distribution results of fault zones (obtained based on seismic profile interpretation) are transformed into coordinates and unified into the same reference frame; finally, the real-time location results of microfractures (obtained through microseismic monitoring or acoustic emission technology) are introduced and superimposed on the model in the form of point clouds. Layered visualization of multi-source data is achieved through color mapping (e.g., red represents high-temperature areas and blue represents low-temperature areas) and transparency adjustment (e.g., fault zones are set to semi-transparent). The final generated three-dimensional visualization image in this embodiment contains four types of information layers: the basic geological layer (showing the geological structure derived from the inversion of the density field and magnetic susceptibility field), the physical property parameter layer (displaying the spatial distribution of the porosity field and permeability field), the thermal field layer (presenting the correspondence between the temperature field and the thermal anomaly center), and the dynamic monitoring layer (marking the real-time location and propagation direction of microcracks). The information fusion between the layers is achieved by overlaying the transparency (e.g., setting the transparency of the physical property parameter layer to 70%).

[0051] S602. Extract key parameters from the three-dimensional joint distribution model, combine them with the porosity probability distribution model for statistical analysis, and generate parameter histograms and spatial distribution rose diagrams. Key parameters include thermal reservoir temperature gradient, porosity distribution range, and permeability spatial variation coefficient. In step S602, it is important to explain in detail that extracting key parameters from the three-dimensional joint distribution model is a crucial step in the quantitative analysis of the characteristics of the thermal reservoir. Specifically, firstly, the temperature gradient of the thermal reservoir is extracted for the temperature field. By calculating the ratio of the temperature difference at different depths to the vertical distance, the rate of temperature change with depth is obtained. Secondly, for the porosity field, its distribution range is determined. The porosity values ​​of all grid nodes are statistically analyzed, and the minimum, maximum, and quartiles are calculated to define the fluctuation range of porosity. Thirdly, for the permeability field, the spatial variation coefficient is calculated. First, the standard deviation of permeability for each grid node is obtained, and then divided by the mean permeability to quantify the spatial dispersion of permeability. When conducting statistical analysis using a porosity probability distribution model, the extracted key parameters are correlated with the model (generated through previous steps, reflecting the probability density of porosity in different geological units). A parameter histogram is plotted with porosity on the horizontal axis and probability density on the vertical axis to visually demonstrate the distribution pattern of porosity (e.g., normal or skewed distribution). Simultaneously, with the study area center as the origin, sector-shaped regions are divided according to azimuth. The mean values ​​of key parameters (e.g., temperature gradient, permeability coefficient of variation) within each sector are statistically analyzed, and a spatial distribution rose diagram is plotted to reveal the spatial anisotropy of the parameters (e.g., temperature gradients near fault zones are significantly higher than in other areas). The final parameter histogram and spatial distribution rose diagram can help determine the heterogeneity of the thermal reservoir (e.g., whether high-porosity areas are concentrated on both sides of the fault zone) and the thermal conductivity of the reservoir (e.g., the correlation between temperature gradient and permeability coefficient of variation).

[0052] S603. Combine the geological structure background of the target area to perform geological interpretation of the three-dimensional joint distribution model, identify the dominant channels and closed boundaries of the thermal reservoir, and calculate the geothermal resource reserves based on the thermal reservoir volume, temperature gradient and fluid properties of the target area. In step S603, it is necessary to explain in detail the geological interpretation of the three-dimensional joint distribution model in conjunction with the geological structural background of the target area. In practice, this first requires a detailed analysis of the geological structural background of the target area, including the distribution of fault zones, stratigraphic lithology, and tectonic activity history. Next, using the geological structural field (such as fault zone identification results) in the three-dimensional joint distribution model, combined with physical property fields (porosity field, permeability field) and temperature field, the dominant channels of the geothermal reservoir are identified. Dominant channels are typically characterized by high porosity, high permeability, and high temperature; these areas are the main pathways for geothermal fluid migration. Simultaneously, by analyzing the sealing characteristics of the geological structural field, such as the closed sections of anticlines or the closed segments of fault zones, the sealing boundaries of the geothermal reservoir are determined. These boundaries limit the migration range of geothermal fluids. After identifying the dominant channels and closed boundaries, the geothermal resource reserves are further calculated using geothermal resource reserve calculation formulas (such as the volumetric method and the thermal reservoir method) by combining the geothermal reservoir volume (obtained through a three-dimensional joint distribution model), temperature gradient (extracted from the temperature field), and fluid properties (such as fluid density and specific heat capacity, which can be obtained through experimental measurements or literature data) of the target area. This calculation process not only considers the geometric characteristics of the geothermal reservoir but also fully considers its thermophysical properties, thus the calculation results are more accurate and reliable.

[0053] S604. Generate exploration results report based on three-dimensional visualization results images, parameter histograms, spatial distribution rose diagrams, and geothermal resource reserves.

[0054] In step S604, it is necessary to explain in detail that, firstly, the three-dimensional visualization of the results is presented as the core part of the report. This image presents the basic geological structure, physical parameter distribution, thermal field characteristics, and dynamic monitoring information (such as the location of microfractures) of the geothermal reservoir in an intuitive three-dimensional way. Secondly, parameter histograms and spatial distribution rose diagrams are used as supplements to the quantitative analysis. The histograms reveal the distribution patterns of key parameters (such as porosity, temperature gradient, and permeability variation coefficient) through statistical methods, while the rose diagrams show the anisotropic characteristics of the parameters through spatial orientation analysis. The combination of the two can effectively assess the heterogeneity and thermal conductivity efficiency of the geothermal reservoir. Finally, the calculated geothermal resource reserves are presented as the key conclusion of the report. This result is based on the geothermal reservoir volume, temperature gradient, and fluid properties obtained from the three-dimensional joint distribution model. During the report writing process, it is necessary to ensure that the data sources of each part are clear, the analysis methods are reasonable, and the conclusions are well-founded. At the same time, a layout combining text and graphics should be used to improve the readability of the report, resulting in the final exploration results report.

[0055] Example 2 A comprehensive geothermal resource exploration system, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to implement a comprehensive geophysical exploration method for geothermal resources when executing executable instructions.

[0056] It should be noted that the computer device includes a processor, a memory, and may also include one or more of a multimedia component, an input / output (I / O) interface, and a communication component.

[0057] The processor controls the overall operation of the computer device to complete all or part of the steps in the above-mentioned comprehensive geophysical exploration method for geothermal resources.

[0058] Memory is used to store various types of data to support the operation of the computer device. This data may include, for example, instructions for any application or method used to operate on the computer device, as well as application-related data. Memory can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read-Only Memory (EPROM), Programmable Read-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk.

[0059] The multimedia component may include a screen and an audio component, wherein the screen may be, for example, a touch screen, and the audio component is used to output and / or input audio signals; for example, the audio component may include a microphone for receiving external audio signals, the received audio signals may be further stored in memory or transmitted via a communication component; the audio component may also include at least one speaker for outputting audio signals.

[0060] I / O interfaces provide interfaces between the processor and other interface modules, such as keyboards, mice, buttons, etc.; these buttons can be virtual buttons or physical buttons.

[0061] The communication component is used for wired or wireless communication between the computer device and other devices; wireless communication, such as Wi-Fi, Bluetooth, Near Field Communication (NFC), 2G, 3G, 4G or 5G, or one or more combinations thereof, and the corresponding communication component may include: Wi-Fi module, Bluetooth module, NFC module, mobile communication module.

[0062] As a preferred embodiment, the computer device may be implemented by one or more application-specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field-programmable gate arrays (FPGAs), controllers, microcontrollers, microprocessors, or other electronic components to execute the above-described geothermal resource integrated geophysical exploration method.

[0063] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A comprehensive geophysical exploration method for geothermal resources, characterized in that, The method includes: S1. Collect surface thermal anomaly data of the target area, and use satellite and calibration algorithms to invert the surface temperature based on the surface thermal anomaly data to generate a thermal anomaly distribution map; S2. The UAV magnetic / gravity exploration system collects regional tectonic data of the target area. Based on the regional tectonic data, after diurnal variation correction, it is jointly inverted with gravity data to obtain the spatial distribution results of the fault zone. S3. Collect dynamic signals of microcracks in the target area to obtain real-time microcrack location results; S4. Collect high-temperature logging data from wells in the target area and generate a three-dimensional image of the porosity distribution of the thermal reservoir. S5. Based on the joint inversion algorithm, the spatial distribution results of the fault zone, the real-time location results of micro-fractures and the three-dimensional image of the porosity distribution of the thermal reservoir are fused. The joint distribution of physical parameters of the underground medium is simulated by a generative adversarial network. Furthermore, geological prior knowledge is introduced to constrain the multiple solutions of the inversion and a three-dimensional joint distribution model of the temperature field, porosity and permeability of the thermal reservoir is constructed. S6. Based on the aforementioned three-dimensional joint distribution model and thermal anomaly distribution map, output the exploration results report.

2. The comprehensive geophysical exploration method for geothermal resources as described in claim 1, characterized in that, Generating the thermal anomaly distribution map in step S1 includes: S101. Collect raw data of surface thermal radiation in the target area; S102. Perform radiometric calibration and atmospheric correction on the raw surface thermal radiation data to eliminate response differences and the influence of atmospheric transmittance, and obtain the true surface thermal radiation brightness value. S103. Based on the actual surface thermal radiation brightness value and combined with the dynamic estimation results of surface emissivity, the split window algorithm is used to retrieve surface temperature and generate continuous surface temperature spatial distribution data. S104. Spatial overlay analysis is performed on continuous surface temperature spatial distribution data, elevation data, and land cover type data. Continuous surface temperature raster images are generated through interpolation processing. Surface thermal anomaly data are identified based on temperature gradient thresholds to obtain thermal anomaly areas. S105. Based on the geological structural background, classify and visualize the thermal anomaly areas to generate a thermal anomaly distribution map that conforms to geological exploration specifications.

3. The comprehensive geophysical exploration method for geothermal resources as described in claim 1, characterized in that, The results of obtaining the spatial distribution of the fault zone in step S2 include: S201. Use an unmanned aerial vehicle (UAV) magnetic exploration system to collect magnetic gradient data of the target area and obtain raw magnetic anomaly data; S202. Perform diurnal variation correction processing on the original magnetic anomaly data to eliminate the influence of diurnal variation of the geomagnetic field and obtain corrected magnetic anomaly data. S203. Simultaneously use the UAV gravity exploration system to collect gravity anomaly data of the target area, and perform terrain correction and Bouguer correction to obtain gravity gradient data; S204. The corrected magnetic anomaly data and gravity gradient data are jointly inverted. The least squares inversion algorithm is used to fuse the magnetic and gravity data to generate a model of underground density and magnetic susceptibility distribution. S205. Based on the aforementioned underground density and magnetic susceptibility distribution model, linear anomaly zones are identified using an edge detection algorithm. Combined with geological structural feature analysis, a three-dimensional image of the spatial distribution of the fault zone is generated, and the spatial distribution result of the fault zone is output.

4. The comprehensive geophysical exploration method for geothermal resources as described in claim 3, characterized in that, The expression for the least squares inversion algorithm is: in, Let represent the least squares objective function. This indicates the total number of observation points collected by the UAV magnetic exploration system. Indicates the first Corrected magnetic anomaly data for each observation point Indicates the first Forward modeling values ​​of magnetic anomalies at each observation point This indicates the total number of observation points collected by the UAV gravity exploration system. Indicates the first Gravity gradient data at each observation point Indicates the first The forward modeling values ​​of the gravity gradient at each observation point Represents the regularization coefficient. Represents the spatial gradient of the density model. This represents the spatial gradient of the magnetic susceptibility model. This represents a model of underground magnetic susceptibility distribution. This represents the spatial coordinates of the model parameters within the subsurface medium. Represented as the magnetic dipole field kernel function, Represents the spatial coordinates of a surface observation point. Represents the underground density distribution model. This represents the kernel function of the Newtonian gravitational gradient.

5. The comprehensive geophysical exploration method for geothermal resources as described in claim 1, characterized in that, The generation of a three-dimensional image of the porosity distribution of the thermal reservoir in step S4 includes: S401. Synchronously acquire high-temperature logging data in the well. The high-temperature logging data includes wellbore axial temperature profile data and sonic transit time data, and record well depth positioning information. S402. Perform wellbore compensation correction and formation invasion correction on the sonic time difference data, and calculate the in-situ sonic velocity of the thermal reservoir in combination with the temperature profile data. S403. Based on the in-situ acoustic velocity of the thermal reservoir, establish an acoustic velocity-porosity conversion function, and use the Kriging interpolation algorithm to extrapolate the discrete point acoustic velocity data to three-dimensional space. S404. Based on the acoustic velocity-porosity conversion function, the temperature gradient abrupt point and the acoustic velocity anomaly region are fused, and the measurement error interference is eliminated by Bayesian joint inversion to generate a porosity probability distribution model. S405. Construct a radial basis function grid centered on the borehole, convert the porosity probability distribution model into three-dimensional volume data of porosity isosurface, and output a three-dimensional image of the porosity distribution of the thermal reservoir.

6. The comprehensive geophysical exploration method for geothermal resources as described in claim 5, characterized in that, The expression for the acoustic velocity-porosity conversion function is: in, Indicates underground point Porosity at that location Indicates the sound wave velocity in the rock matrix. Indicates the in-situ acoustic velocity of the thermal reservoir. Indicates the sound wave velocity of pore fluid. Indicates the temperature sensitivity coefficient. This represents the axial temperature gradient of the wellbore.

7. The comprehensive geophysical exploration method for geothermal resources as described in claim 1, characterized in that, The construction of a three-dimensional joint distribution model of the temperature field, porosity, and permeability of the thermal reservoir in step S5 includes: S501. Based on the spatial distribution results of the fracture zone, the real-time location results of micro-fractures and the three-dimensional image of the porosity distribution of the thermal reservoir, a three-dimensional mesh alignment algorithm is used to unify the spatial coordinate system. Discrete point data is converted into a continuous data field through nearest neighbor interpolation to form a multi-source data fusion preprocessing result. S502. Based on the multi-source data fusion preprocessing results, a joint inversion objective function is constructed using a cross-gradient constraint inversion algorithm. The spatial distribution of the fault zone is used as the construction constraint, the micro-fracture location results are used as the dynamic boundary conditions, and the porosity distribution is used as the initial model. By minimizing the objective function, a multi-physics joint inversion is performed to generate an initial distribution model of the physical parameters of the underground medium. S503. Based on the initial distribution model of the physical parameters of the underground medium, a generative adversarial network is used to simulate the joint distribution of the physical parameters of the underground medium. The generator generates a simulated parameter set with the initial model as input. The discriminator evaluates the matching degree between the parameter set and the real data based on the actual geological data. The inversion parameter set is optimized through iterative training until the matching degree reaches the standard, and the optimized inversion parameter set is obtained. S504. Based on the optimized inversion parameter set, geological prior knowledge is introduced to constrain the inversion ambiguity. The geological prior probability distribution is determined by combining the thermal anomaly distribution map and the fault zone distribution results. The data fitting and prior constraints are balanced by the Bayesian inversion algorithm to eliminate the inversion ambiguity and obtain the eliminated inversion parameter set. S505. Based on the eliminated inversion parameter set, construct a three-dimensional joint distribution model of the thermal reservoir temperature field, porosity, and permeability.

8. The comprehensive geophysical exploration method for geothermal resources as described in claim 7, characterized in that, The expression for the joint inversion objective function is: in, Describes the joint inversion objective function. , , Represents the regularization coefficient. Indicates the first Observational data, This represents the forward-modeled value. Represents the gradient vector of physical parameters of the underground medium. Represents the gradient vector constrained by geological structure. This represents a priori geological model.

9. The comprehensive geophysical exploration method for geothermal resources as described in claim 1, characterized in that, In step S6, based on the three-dimensional joint distribution model, the exploration results report is output, including: S601. Based on the three-dimensional joint distribution model, perform multi-source data fusion visualization, and overlay the thermal anomaly distribution map, the spatial distribution result of the fault zone and the real-time location result of the microcracks to generate a three-dimensional visualization result image. S602. Extract key parameters from the three-dimensional joint distribution model, perform statistical analysis in conjunction with the porosity probability distribution model, and generate parameter histograms and spatial distribution rose diagrams. Key parameters include thermal reservoir temperature gradient, porosity distribution range, and permeability spatial variation coefficient. S603. Combine the geological structure background of the target area to perform geological interpretation of the three-dimensional joint distribution model, identify the dominant channels and closed boundaries of the thermal reservoir, and calculate the geothermal resource reserves based on the thermal reservoir volume, temperature gradient and fluid properties of the target area. S604. Based on the three-dimensional visualization results image, parameter histogram, spatial distribution rose diagram, and geothermal resource reserves, generate an exploration results report.

10. A comprehensive geophysical exploration system for geothermal resources, characterized in that, The system includes: processor; Memory used to store processor-executable instructions; The processor is configured to implement the geothermal resource comprehensive geophysical exploration method according to any one of claims 1 to 9 when executing the executable instructions.

Images