Geothermal energy prediction method based on multi-source data

By integrating multi-source data fusion and feature gradient analysis, combined with spatial connectivity and regional growth algorithms, the problems of accuracy and continuity in geothermal energy distribution prediction were solved, enabling precise geothermal energy exploration and utilization in complex geological contexts.

CN121480783BActive Publication Date: 2026-06-09CHINA UNIV OF MINING & TECH (BEIJING)

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
CHINA UNIV OF MINING & TECH (BEIJING)
Filing Date
2025-09-15
Publication Date
2026-06-09

Smart Images

  • Figure CN121480783B_ABST
    Figure CN121480783B_ABST
Patent Text Reader

Abstract

The present application provides a kind of geothermal energy prediction method based on multi-source data, it is related to geothermal energy exploration technical field, the present application is by in target area according to preset grid scale carries out regular grid division, and with grid center point synchronous acquisition ground temperature observation value and ground heat flow observation value, utilize minimum-maximum normalization method to standardize multi-source observation data, construct two-dimensional feature vector reflecting geothermal characteristic state, realized the fusion and space standardization processing of multi-source data, based on the feature vector difference of non-boundary grid in grid array, adopt the automatic calculation of spatial feature gradient using Euclidean distance criterion to determine edge recognition threshold, filter spatial mutation area as edge candidate point, and further combined with spatial eight connectivity analysis, remove isolated distribution edge candidate point, obtain spatial continuous and physically reasonable boundary point set, using region growing algorithm to expand geothermal distribution area realizes the recursive generation of geothermal distribution area.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of geothermal energy exploration technology, specifically a geothermal energy prediction method based on multi-source data. Background Technology

[0002] Geothermal energy, as an important clean energy source, has broad application prospects in new energy development and regional energy planning. Especially in areas with frequent geological activity and abundant geothermal resources, scientifically and accurately identifying the distribution areas and boundaries of geothermal resources is of great significance for resource evaluation, development and utilization, and environmental protection. For a large area with great geothermal resource potential, such as basins, rift zones, the periphery of geothermal fields, or new urban areas, due to the vast area, complex geological structure, and uneven distribution of geothermal anomalies, it is necessary to scientifically and efficiently predict the specific distribution range and boundaries of geothermal anomalies throughout the entire region. Only in this way can we achieve accurate exploration, rational development and efficient utilization of resources, avoid resource waste and environmental risks, and provide basic spatial data support for subsequent drilling site layout and engineering design.

[0003] Most existing technologies rely on single physical quantities such as surface temperature or heat flow and expert experience to manually delineate geothermal areas, or use simple threshold methods for spatial division. These methods are difficult to balance multi-source data fusion, spatial feature changes, and automated processing. They are also susceptible to measurement errors, data noise, and human subjectivity, resulting in inaccurate identification of anomaly boundaries, discontinuous spatial distribution, and difficulty in adapting to the prediction needs of complex geological backgrounds and large-scale areas. To address the above practical needs and technical problems, a technical method is needed that can fully integrate multi-source data such as surface temperature and surface heat flow to identify geothermal area boundaries, realize spatial connectivity analysis, and perform regional recursive expansion, so as to improve the accuracy of geothermal energy distribution prediction.

[0004] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] (a) Technical problems to be solved

[0006] To address the shortcomings of existing technologies, this invention provides a geothermal energy prediction method based on multi-source data. Through multi-source data normalization processing and two-dimensional feature vector construction, it fully utilizes the spatial distribution information of surface temperature and surface heat flow to identify geothermal energy anomaly zones. Employing spatial abrupt change detection based on feature gradients, it can objectively identify geothermal distribution boundaries, effectively avoiding the subjective delineation and misjudgment problems of single thresholds in existing technologies. Simultaneously, through spatial octet connectivity analysis and outlier removal, it improves the spatial continuity of the boundary point set. The recursive expansion mechanism of the region growing algorithm, combined with the Euclidean distance criterion of feature vectors, ensures that geothermal distribution areas are physically consistent and spatially coherent, significantly improving the automation and accuracy of geothermal energy distribution prediction and solving the technical problems described in the background art.

[0007] (II) Technical Solution

[0008] To achieve the above objectives, the present invention provides the following technical solution:

[0009] A geothermal energy prediction method based on multi-source data, comprising the following steps:

[0010] The target area to be predicted is divided according to the preset grid scale. The target area is divided into regular grids and the grids are numbered to generate a two-dimensional grid array. Surface temperature and surface heat flow observations are obtained at the center point of each grid.

[0011] The surface temperature and surface heat flow observations of all grids are normalized respectively. The normalized surface temperature observations and normalized surface heat flow observations of each grid are combined to form a two-dimensional feature vector that reflects the geothermal characteristics of the grid.

[0012] For non-boundary grids in a two-dimensional grid array, based on two-dimensional feature vectors and using Euclidean distance, feature vector differences are generated to reflect the degree of feature change between the grid and each adjacent grid. The feature gradient value with the largest feature vector difference is selected as the feature gradient value of the grid. The feature gradient value is used to measure the degree of spatial abrupt change of the grid in geothermal features.

[0013] The average level and fluctuation range of all feature gradient values ​​are statistically obtained. The linear weighted sum of the average level and fluctuation range is used as the edge identification threshold. Each feature gradient value is compared with the edge identification threshold. If the feature gradient value corresponding to the grid is greater than or equal to the threshold, the grid is determined to be a candidate edge point of the geothermal area.

[0014] For all grids identified as edge candidates, connectivity analysis is performed based on spatial octet connectivity. Isolated edge candidates are removed to obtain a set of boundary points. Based on the boundary points, the geothermal distribution area is expanded using a region growing algorithm. The two-dimensional feature vector of the grid is used as the criterion to determine whether adjacent grids meet the region growing conditions, thus generating the predicted geothermal distribution area.

[0015] Furthermore, the target area is a ground-view projection of the surface area to be predicted, and the grid is a square grid. The method of dividing the target area into regular grids is based on a preset grid scale, and the space is evenly divided according to the direction of rows and columns to generate each grid in sequence, and all grids completely cover the target area to generate a two-dimensional grid array.

[0016] At the boundary of the target area, if there are grids that exceed the boundary of the target area, the grids that exceed the target area are removed, and only the grids that completely cover the target area are retained. The direction of the rows and columns is based on the geographical direction, with the northernmost grid as the first row and the row number increasing from north to south, and the westernmost grid as the first column and the column number increasing from west to east.

[0017] When numbering the grids, the northernmost and westernmost grids in the target area are determined as the starting points of the grid array.

[0018] Furthermore, when acquiring surface temperature and surface heat flow observations, surface temperature and surface heat flow observations at all grid center points are collected simultaneously. The collection method is as follows:

[0019] The center point of each grid is determined as the spatial reference for data collection. Surface temperature is observed and recorded at the center point of each grid, and the surface temperature value at that point is taken as the surface temperature observation value of the grid center point. At the same grid center point, the surface heat flow data at that point is continuously recorded, and the average value is taken as the surface heat flow observation value of the grid. The surface temperature observation value and surface heat flow observation value collected at each grid center point are mapped to the grid number.

[0020] Furthermore, for all grids, the minimum and maximum values ​​of surface temperature observations and surface heat flow observations are statistically analyzed. The minimum-maximum normalization method is used to standardize the surface temperature observations and surface heat flow observations of each grid to between 0 and 1.

[0021] The normalized surface temperature observations and surface heat flow observations are combined to form a two-dimensional vector, denoted as the first vector. Line number The two-dimensional eigenvectors of the column grid are The combination method is as follows: ,in, and They represent Line number Surface temperature and surface heat flow observations after grid normalization.

[0022] Furthermore, the non-boundary mesh is a mesh located inside a two-dimensional mesh array, in which there are directly adjacent meshes in all four directions: up, down, left, and right; that is, a non-boundary mesh refers to a mesh that is not the outermost mesh, but is surrounded by other meshes and has a complete four-neighborhood relationship.

[0023] The logic for generating the feature vector difference between the grid and its neighboring grids is as follows:

[0024] For each non-boundary grid, first determine the row and column number of the grid in the array, and then extract the two-dimensional feature vectors of the grid and its four adjacent grids above, below, to the left and to the right in sequence according to the grid number;

[0025] The two-dimensional feature vector of the target grid is compared with the two-dimensional feature vector of its neighboring grids, and the difference in geothermal feature state between the target grid and its neighboring grids is calculated. The difference value is calculated by calculating the numerical difference between the target grid and its neighboring grids in the two dimensions of normalized surface temperature and normalized surface heat flux, and then taking the square root of the sum of the squares of the differences in the two dimensions to obtain the degree of feature change between the target grid and each neighboring grid, which is labeled as the feature vector difference value.

[0026] After obtaining the feature vector difference values ​​in four directions for each non-boundary grid, the largest value among the feature vector difference values ​​in the four directions is selected as the feature gradient value of that grid.

[0027] Furthermore, the logic for obtaining the edge recognition threshold by statistically analyzing all feature gradient values ​​is as follows:

[0028] The characteristic gradient values ​​of all non-boundary grids in the two-dimensional grid array are statistically analyzed, and the average level and standard deviation of all characteristic gradient values ​​are calculated. The average level and fluctuation amplitude of the characteristic gradient values ​​are combined in a linear weighted manner to obtain the identification threshold for judging the edge of the geothermal area.

[0029] The feature gradient value of each grid is compared with the edge recognition threshold. For grids whose feature gradient value is greater than or equal to the edge recognition threshold, the grid is determined to be a candidate edge point of the geothermal area.

[0030] Furthermore, when removing spatially isolated edge candidate points to obtain boundary points, all grids that are determined to be edge candidate points are traversed. For each edge candidate point, its spatial position in the two-dimensional grid array is checked, and its adjacency relationship with other edge candidate points is identified.

[0031] For each edge candidate point, analyze its neighborhood in eight directions, namely the positions of the eight adjacent grids: above, below, left, right, and upper left, upper right, lower left, and lower right, and determine whether there are other grids in these neighborhood positions that are identified as edge candidate points.

[0032] If no other edge candidate points are found in the neighborhood of a certain edge candidate point in any of its eight directions, then the point is considered to be spatially isolated. Isolated edge candidate points are eliminated, and the remaining edge candidate points are marked as boundary points.

[0033] Furthermore, for each boundary point, check whether there are other boundary points in its neighborhood in eight directions, and group all boundary points that can be connected to each other through their neighborhood in eight directions into the same connected component. After traversing all boundary points, obtain all connected components.

[0034] For each connected component, if each boundary point in the connected component has at least two other boundary points connected to it in its neighborhood in eight directions, then the connected component forms a closed boundary. If there is a boundary point in the connected component that is connected to only one other boundary point in its neighborhood in eight directions, then the connected component is a non-closed component.

[0035] For all non-closed branches, identify endpoint boundary points, where each endpoint boundary point is a boundary point connected to only one other boundary point. In each non-closed branch, extract all endpoint boundary points.

[0036] For all endpoint boundary points, select the endpoint boundary point with the smallest spatial distance as a pair of endpoint boundary points. For each pair of endpoint boundary points, select the grid on the shortest spatial path between the endpoint boundary points and mark it as a boundary point. Count all boundary points to form a boundary point set.

[0037] Furthermore, the method for expanding the geothermal distribution area based on boundary points using a region growing algorithm to generate the predicted geothermal distribution area is as follows:

[0038] All boundaries in the set of boundary points are used as starting seed points. For each starting seed point, all neighboring grids that have not yet been marked as boundary points in its eight directions are identified as current expansion candidate grids to be determined. For each expansion candidate grid, its two-dimensional feature vector is extracted, and the Euclidean distance between it and the two-dimensional feature vector of the current starting seed point is calculated.

[0039] If the Euclidean distance is less than the preset growth threshold, the candidate grid is determined to meet the regional growth conditions and can be included in the geothermal distribution area. If the Euclidean distance is greater than the growth threshold, the candidate grid is determined not to be included in the geothermal distribution area for the time being.

[0040] For all candidate grids that meet the growth conditions, they are assigned to the geothermal distribution area and used as new growth seed points. For the newly assigned growth seed points, the above steps of expanding candidate grid determination, calculating the Euclidean distance between two-dimensional feature vectors, and determining the regional growth conditions are repeated recursively until no new grids meet the regional growth conditions during the regional growth process. At this point, the regional growth algorithm terminates. All assigned grids and boundary points in the boundary point set together constitute the predicted geothermal distribution area.

[0041] (III) Beneficial Effects

[0042] This invention divides the target area into regular grids according to a preset grid scale, and simultaneously acquires surface temperature and surface heat flow observations at the grid center point. It standardizes the multi-source observation data using the minimum-maximum normalization method, constructs a two-dimensional feature vector reflecting the geothermal characteristic state, and realizes the fusion and spatial standardization processing of multi-source data. Based on the feature vector differences of non-boundary grids within the grid array, it automatically calculates the spatial feature gradient using the Euclidean distance criterion, and determines the edge identification threshold through statistical methods, accurately screening spatial abrupt change areas as edge candidate points. This solves the problems of strong subjectivity and inaccurate boundary identification in the prior art.

[0043] This invention further combines spatial octane connectivity analysis to remove isolated edge candidate points and obtain a spatially continuous and physically reasonable set of boundary points. On this basis, a region growing algorithm is used to expand the geothermal distribution area, and the region affiliation is dynamically determined by the Euclidean distance of the feature vector, thus realizing the recursive generation of the geothermal distribution area. Attached Figure Description

[0044] Figure 1 This is a schematic diagram of the overall method flow of the present invention;

[0045] Figure 2 This is a schematic diagram illustrating the division of the target region to be predicted in this invention;

[0046] Figure 3 This is a schematic diagram of the grid extending beyond the target area boundary in this invention;

[0047] Figure 4 This is a schematic diagram of the four-neighbor and eight-neighbor neighbor grids of the grid in this invention;

[0048] Figure 5 This is a schematic diagram of the division of connected branches in this invention. Detailed Implementation

[0049] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0050] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0051] Example:

[0052] Please see Figures 1 to 5 The present invention provides a technical solution:

[0053] A geothermal energy prediction method based on multi-source data, comprising the following steps:

[0054] S1: Divide the target area to be predicted according to the preset grid scale, divide the target area into regular grids, number the grids and generate a two-dimensional grid array, and obtain the surface temperature observation value and surface heat flow observation value at the center point of each grid.

[0055] In this embodiment, the target area is a ground-view projection of the surface area to be predicted. The surface area to be predicted refers to an area that is preliminarily determined to have a high probability of having geothermal energy resources within the surface space based on preliminary exploration, existing geological data, geophysical anomalies, remote sensing images, historical development records, or other scientific analysis methods. This spatial range is the object used by the method of this invention for subsequent geothermal energy distribution prediction and boundary identification. The grid is a square grid. When dividing the target area into regular grids, a preset grid scale is used as a reference. The space is evenly divided according to the direction of rows and columns, and each grid is generated in sequence. All grids completely cover the target area to generate a two-dimensional grid array.

[0056] At the boundary of the target area, if there are grids that exceed the boundary of the target area, the grids that exceed the target area are removed, and only the grids that completely cover the target area are retained. The direction of the rows and columns is based on the geographical direction, with the northernmost grid as the first row and the row number increasing from north to south, and the westernmost grid as the first column and the column number increasing from west to east.

[0057] When numbering the grid, the northernmost and westernmost grids of the target area are designated as the starting points of the grid array. Ground-view projection refers to projecting the actual spatial extent of a surface area from a direction perpendicular to the ground, resulting in a two-dimensional spatial extent. The actual surface area may have hills, valleys, and other undulating terrain. Ground-view projection only focuses on the projected extent of this area on the map, such as using polygons or rectangles to delineate it, regardless of elevation; only the distribution on the two-dimensional plane is considered. Using ground-view projection as the target area is beneficial for spatial data grid division, numbering, data acquisition, and subsequent analysis. It simplifies the complexity of spatial modeling and facilitates the use of standard spatial analysis methods such as two-dimensional grid arrays for geothermal energy prediction.

[0058] The target area refers to the top-down projection of the surface area to be predicted. When dividing the target area into grids, the spatial boundary of the target area is first determined according to the actual engineering requirements and the scope of the survey. Then, according to the pre-set grid scale, the entire target area is uniformly divided into several square grids within the ground projection range. The division process is based on spatial coordinates, generating rows and columns according to the north-south and east-west directions respectively to ensure that the grid arrangement is regular and the spatial coverage is complete. When dividing the grid, the scheme with the largest sum of the areas of all retained grids is selected as the final grid division result. During the grid division process, there may be multiple different grid arrangement methods. Each arrangement method will result in a different number of grids completely located inside the target area. In order to make the target area as completely covered by the grid as possible, this embodiment tries multiple grid division methods and calculates the total area of ​​all retained grids under each method. Finally, the grid division method that maximizes the sum of the areas of these retained grids is selected as the final grid division result.

[0059] In other words, by optimizing the grid arrangement, more grids can completely cover the target area, reducing waste or omissions at the boundaries, thereby improving the spatial integrity and data utilization of prediction and analysis. This approach helps ensure more comprehensive and efficient data collection and spatial analysis within the target area.

[0060] Specifically, the northernmost and westernmost boundaries of the target area are used as the starting points of the two-dimensional grid array. The northernmost grid is set as the first row, and the row numbers increase sequentially from north to south. The westernmost grid is set as the first column, and the column numbers increase sequentially from west to east. Following this spatial arrangement logic, all grids are generated sequentially, and each grid is assigned a row and column number to achieve spatially ordered management and data indexing. For grids at the boundary of the target area, if some grids are found to exceed the boundary of the target area, after generating the grid array, these grids that exceed the boundary are removed by a spatial discrimination algorithm, and only grids that completely cover the interior of the target area are retained to ensure the integrity and standardization of the spatial division of the target area.

[0061] When a grid extends beyond the target area's boundary, it means that when the target area is regularly divided according to a preset grid scale, the actual spatial boundary of the target area usually does not strictly coincide with the edge of the grid array. Therefore, at the boundary, some grids may partially or completely extend beyond the target area. The spatial position of these grids is not entirely within the defined range of the target area, and they exceed or cross the target area's boundary. For example, if the target area is an irregular polygon, and the grid is divided into square units, then at the edge, some grids will often have parts of their area inside the target area, while others will be outside, or even some grids may fall entirely outside the area.

[0062] A grid that completely covers the target area refers to a grid whose spatial extent or area boundaries are entirely within the spatial boundaries set for the target area. Geometrically, all four sides and all points within this type of grid fall within the target area, with no part exceeding the target area's boundaries. Only this type of grid can ensure that its spatial data belongs to the target area, guaranteeing the accuracy of data analysis and the scientific validity of its spatial distribution.

[0063] The spatial extent of each grid is compared with the boundary of the target area. If the entire spatial extent of the grid is inside the target area, it is determined to be a grid that completely covers the target area; if the spatial extent of the grid is partially or entirely outside the target area, it is determined to be a grid that exceeds the boundary of the target area and is discarded.

[0064] Furthermore, when acquiring surface temperature and surface heat flow observations, surface temperature and surface heat flow observations at all grid center points are collected simultaneously. The collection method is as follows:

[0065] The center point of each grid is determined as the spatial reference for data collection. Surface temperature is observed and recorded at the center point of each grid, and the surface temperature value at that point is taken as the surface temperature observation value of the grid center point. At the same grid center point, the surface heat flow data at that point is continuously recorded, and the average value is taken as the surface heat flow observation value of the grid. The surface temperature observation value and surface heat flow observation value collected at each grid center point are mapped to the grid number.

[0066] Based on the generated regular two-dimensional grid array, the spatial coordinates of the center point of each grid cell are determined. This center point serves as the reference point for the spatial distribution of the grid and is the foundation for subsequent data acquisition and spatial analysis. All temperature and heat flow data are archived using the center point as the spatial index. For each grid center point, when collecting surface temperature observations, surface temperature information must be obtained from the spatial location corresponding to that center point. Existing geological survey data, remote sensing temperature distribution maps, geophysical measurement results, historical meteorological data, or other multi-source databases can be used to find the surface temperature observation record closest to the center point. If direct observation data is available for the center point, it is directly selected; otherwise, spatial interpolation methods (such as inverse distance weighting, Kriging interpolation, etc.) can be used to estimate and complete the data using surrounding known temperature points to ensure that representative temperature values ​​are obtained for each grid center point.

[0067] The logic for collecting surface heat flow observations is similar to that for temperature collection. For each grid center point, existing surface heat flow data (such as geothermal survey reports, heat flow density distribution maps, and relevant scientific databases) are queried to obtain heat flow data for that point or its neighboring area. If there are direct observations or estimations for the center point, they are used directly. If there is no direct data, spatial interpolation or geological model inference is used, combined with reasonable supplementation from surrounding points with available data, to obtain the surface heat flow observation value for the grid center point. For continuously recorded data, the average value of heat flow data over a certain time interval can be taken to improve the representativeness of the estimation.

[0068] All collected surface temperature and heat flow observations are mapped and stored one-to-one with the spatial coordinates and numbers of the grid center points, forming a standardized data table or data structure. If, due to special terrain, obstacles, or sampling conditions, some grid center points cannot be directly obtained with measured data, the following technical remedial measures can be adopted: First, find valid observation data from neighboring grids around the center point, and estimate the temperature and heat flow values ​​of the center point through spatial interpolation, weighted averaging, or geological model inference. If the surrounding data also cannot meet the requirements, the point can be marked as "missing" or "no data" in the data table, and the missing values ​​can be reasonably processed (such as interpolation to complete) during subsequent analysis to ensure the integrity of the spatial distribution.

[0069] S2: Normalize the surface temperature and surface heat flow observations for all grids. Combine the normalized surface temperature and surface heat flow observations for each grid to form a two-dimensional feature vector that reflects the geothermal characteristics of that grid.

[0070] In this embodiment, for all grids, the minimum and maximum values ​​of surface temperature observations and surface heat flow observations are statistically analyzed. A minimum-maximum normalization method is used to standardize the surface temperature and heat flow observations for each grid to a range of 0 to 1. The formula used is as follows:

[0071]

[0072]

[0073] in, and They represent Line number Surface temperature and surface heat flow observations after grid normalization. and They represent Line number Before grid normalization, i.e., the original measured surface temperature and surface heat flow observations, and These represent the maximum and minimum values ​​of surface temperature observations across all grids, respectively. and These represent the maximum and minimum values ​​of surface heat flow observations across all grids, respectively.

[0074] The normalized surface temperature observations and surface heat flow observations are combined to form a two-dimensional vector, denoted as the first vector. Line number The two-dimensional eigenvectors of the column grid are The combination method is as follows: ,in, and They represent Line number The surface temperature and surface heat flow observations are normalized after grid processing. The surface temperature observation reflects the thermal state of the surface layer at the center point of the grid and is an important basic data for judging the intensity and distribution of geothermal activity. The surface heat flow observation refers to the rate of geothermal energy transfer to the surface per unit area per unit time at the center point of the grid in the target area, obtained through data collection or scientific estimation. It is usually expressed in milliwatts per square meter and represents the heat flow density of the surface at that point. The surface temperature observation reflects the thermal state of the surface, and the surface heat flow observation measures the geothermal energy output capacity. The combination of the two can comprehensively describe the geothermal characteristics of the grid and is the basic data for geothermal energy distribution prediction and resource evaluation.

[0075] Geothermal resource areas are typically characterized by higher surface temperatures and significantly higher surface heat flow densities than surrounding background areas. This is because underground heat energy is released through the surface, causing these physical quantities to exhibit significant spatial anomalies. When surface temperature or heat flow transitions from a high-value area to a background area, there are often abrupt changes or a significant increase in the spatial distribution gradient. Within the resource area, temperature and heat flow changes are relatively gradual, while at the boundary of the anomaly area, the values ​​of these two indicators will change significantly in space.

[0076] In this embodiment, after normalizing the temperature and heat flow observations, the feature vector difference and feature gradient of each grid in the grid array are calculated. If a certain grid has the largest difference in these two indicators with the surrounding grids, then the point is very likely to be located in the transition zone between the geothermal anomaly zone and the background zone, that is, the "boundary" position of the geothermal area.

[0077] S3: For non-boundary grids in a two-dimensional grid array, based on two-dimensional feature vectors and using Euclidean distance, feature vector differences are generated to reflect the degree of feature change between the grid and each adjacent grid. The feature gradient value with the largest feature vector difference is selected as the feature gradient value of the grid. The feature gradient value is used to measure the degree of spatial abrupt change of the grid in geothermal features.

[0078] By calculating the Euclidean distance between each non-boundary grid and its four neighboring grids in the two-dimensional feature space (i.e., normalized surface temperature and normalized surface heat flux), the spatial variation amplitude of local geothermal features can be keenly captured. By selecting the maximum difference as the feature gradient value, areas with abrupt spatial changes can be effectively located, namely the boundaries and transition zones of geothermal anomaly zones. At the same time, considering the two key physical quantities of temperature and heat flux, misjudgments caused by the influence of a single indicator on local anomalies or measurement errors can be avoided.

[0079] In this embodiment, the non-boundary grid is a grid located inside a two-dimensional grid array, which has directly adjacent grids in all four directions: up, down, left, and right; that is, a non-boundary grid refers to a grid that is not the outermost grid, but is surrounded by other grids and has a complete four-neighbor relationship.

[0080] The logic for generating the feature vector difference between the grid and its neighboring grids is as follows:

[0081] For each non-boundary grid, first determine the row and column number of the grid in the array, and then extract the two-dimensional feature vectors of the grid and its four adjacent grids above, below, to the left and to the right in sequence according to the grid number;

[0082] The two-dimensional feature vector of the target grid is compared with the two-dimensional feature vector of its neighboring grids, and the difference in geothermal feature state between the target grid and its neighboring grids is calculated. The difference value is calculated by calculating the numerical difference between the target grid and its neighboring grids in the two dimensions of normalized surface temperature and normalized surface heat flux, and then taking the square root of the sum of the squares of the differences in the two dimensions to obtain the degree of feature change between the target grid and each neighboring grid, which is labeled as the feature vector difference value.

[0083] After obtaining the feature vector difference values ​​in four directions for each non-boundary grid, the largest value among the feature vector difference values ​​in the four directions is selected as the feature gradient value of that grid.

[0084] No. Line number The column grid and the grid above it are the first... Line number The formula used to calculate the difference in eigenvectors between column grids is:

[0085]

[0086] in, Indicates the first Line number The column grid and the grid above it are the first... Line number The difference in eigenvectors between the column grids, and They represent respectively Line number Surface temperature and surface heat flow observations after grid normalization.

[0087] In this embodiment, the target area is first divided into a regular two-dimensional grid array, and each grid is numbered by row and column. For non-boundary grids within the grid array, since there are directly adjacent grids in the four directions (up, down, left, and right), they have a complete four-neighborhood relationship. For these non-boundary grids, it is necessary to analyze the changes in geothermal characteristics between them and their adjacent grids in order to determine the boundaries of the geothermal anomaly zone.

[0088] After obtaining the differences in four eigenvectors, the largest one is selected as the feature gradient value for that grid. This maximum value represents the most significant spatial abrupt change in geothermal characteristics across all directions, and is used for subsequent boundary determination and anomaly identification. The spatial boundary of a geothermal anomaly zone often manifests as abrupt changes in temperature and heat flow in a certain direction. Calculating only one direction may miss the true boundary; therefore, it is necessary to examine all four directions to capture the most significant spatial changes. Selecting the maximum difference value as the feature gradient ensures that each grid cell reflects its maximum geothermal characteristic abrupt change across all possible directions, thus better capturing the true boundary of the geothermal region.

[0089] The reason for collecting feature vector differences in four directions is that the boundary of a geothermal anomaly zone can appear in any direction. If only one direction is analyzed, the complexity of the spatial distribution may cause the true boundary points to be missed. By analyzing all four directions, it can be ensured that the most prominent feature change points can be effectively detected regardless of how the boundary extends or curves, thus accurately reflecting the essential characteristics of the spatial distribution.

[0090] S4: Statistically obtain the average level and fluctuation range of all feature gradient values, and use the linear weighted sum of the average level and fluctuation range as the edge identification threshold. Compare each feature gradient value with the edge identification threshold. If the feature gradient value corresponding to the grid is greater than or equal to the threshold, the grid is determined to be a candidate edge point of the geothermal area.

[0091] By statistically analyzing all feature gradient values, calculating the average level and fluctuation amplitude, and using the linear weighted sum of these two indicators as the edge identification threshold, adaptive spatial boundary discrimination is achieved. This method automatically adjusts the discrimination criteria based on the actual data distribution, considering both the general level of overall geothermal feature changes and the discreteness and abnormal fluctuations of data within the region. In this way, grids exhibiting abrupt changes in geothermal feature distribution can be effectively distinguished, avoiding the uncertainties and limitations caused by manual or subjective threshold settings.

[0092] For areas with relatively stable data distribution, the threshold will not be too high, thus ensuring that sensitive areas are not missed; while for areas with abnormal and drastic fluctuations, the threshold will not be too low, thus reducing misjudgments and noise interference. This processing method has strong versatility and adaptability, and can be adapted to various actual geothermal scenarios and multi-source data conditions, making boundary extraction more reliable.

[0093] In this embodiment, the logic for obtaining the edge recognition threshold by statistically analyzing all feature gradient values ​​is as follows:

[0094] The characteristic gradient values ​​of all non-boundary grids in the two-dimensional grid array are statistically analyzed, and the average level and standard deviation of all characteristic gradient values ​​are calculated. The average level and fluctuation amplitude of the characteristic gradient values ​​are combined in a linear weighted manner to obtain the identification threshold for judging the edge of the geothermal area.

[0095] The feature gradient value of each grid is compared with the edge recognition threshold. For grids whose feature gradient value is greater than or equal to the edge recognition threshold, the grid is determined to be a candidate edge point of the geothermal area.

[0096] The formula used to calculate the average gradient values ​​of all features is:

[0097]

[0098] in, This represents the average level of all feature gradient values. This represents the number of gradient values ​​for all features. Indicates the first The formula used to calculate the standard deviation of all feature gradient values ​​is:

[0099]

[0100] in, The standard deviation of all feature gradient values ​​is represented by the formula for the identification threshold, which is obtained by combining the average level and fluctuation range of the feature gradient values ​​in a linear weighted manner.

[0101]

[0102] in, Indicates the recognition threshold. This represents the identification weighting coefficient, which is a commonly used value pre-set based on past geothermal resource surveys, spatial anomaly detection, or engineering experience in related fields. For example, k might be 1, 1.5, or 2. Generally, a k value of 1 can better balance false negatives and false positives, making it suitable for scenarios with relatively stable data distribution. If the difference between the edge and background areas is large, the value of the identification weighting coefficient can be appropriately increased to improve the strictness of edge detection.

[0103] The method of combining the average level and fluctuation amplitude of feature gradient values ​​in a linear weighted manner integrates the statistical indices of central tendency (average level) and dispersion (standard deviation). It can adaptively adjust the discrimination threshold according to the actual data distribution. Specifically, the mean reflects the normal level of overall geothermal feature changes, while the standard deviation reflects the volatility and anomaly in the spatial distribution. By combining the two through linear weighting, it can avoid both overly lenient approaches caused by using only the mean and overly strict approaches caused by using only extreme values. It can effectively screen out those grids with the most significant changes in geothermal feature distribution, organically distinguish the edges from general areas, and improve the accuracy and adaptability of boundary identification.

[0104] S5: Perform connectivity analysis on all grids identified as edge candidates based on spatial octet connectivity, remove spatially isolated edge candidates to obtain a set of boundary points, expand the geothermal distribution area based on the boundary points using a region growing algorithm, and determine whether adjacent grids meet the region growing conditions based on the two-dimensional feature vector of the grid to generate the predicted geothermal distribution area.

[0105] For all grid points identified as edge candidates, spatial octet connectivity analysis is performed to remove spatially isolated edge candidates, thereby obtaining a set of boundary points with physical meaning and spatial continuity. This approach significantly improves the accuracy of anomaly zone boundary identification, ensuring that the starting boundary points of subsequent region growing algorithms possess genuine geothermal anomaly zone distribution characteristics. Traditional geothermal boundary identification often relies on manual experience, single threshold methods, or simple spatial screening, which can easily lead to isolated, broken, or discontinuous boundary points, resulting in inaccurate anomaly zone boundary determination and affecting subsequent resource evaluation and development planning. This patented solution systematically identifies and eliminates spatially isolated edge points through octet connectivity analysis, ensuring that the boundary point set forms a connected band or region distribution in space, greatly improving the rationality of spatial distribution and engineering usability. Simultaneously, the region growing algorithm uses boundary points as seed points and recursively expands by fusing the two-dimensional feature vectors of the grid, making dual judgments in both spatial and feature dimensions, thus achieving the identification of geothermal distribution areas.

[0106] In this example, when removing spatially isolated edge candidate points to obtain boundary points, all grids that are determined to be edge candidate points are traversed. For each edge candidate point, its spatial position in the two-dimensional grid array is checked, and its adjacency relationship with other edge candidate points is identified.

[0107] For each edge candidate point, analyze its neighborhood in eight directions, namely the positions of the eight adjacent grids: above, below, left, right, and upper left, upper right, lower left, and lower right, and determine whether there are other grids in these neighborhood positions that are identified as edge candidate points.

[0108] If no other edge candidate points are found in the neighborhood of a certain edge candidate point in any of its eight directions, then the point is considered to be spatially isolated. Isolated edge candidate points are eliminated, and the remaining edge candidate points are marked as boundary points.

[0109] The process iterates through all grid cells identified as edge candidates. For each candidate, its spatial location is determined by its row and column numbers in the 2D grid array. Then, for each of the eight directions of the candidate, the process checks whether adjacent grid cells are also edge candidates. If no other edge candidates are found in any of the adjacent locations in these eight directions, the candidate is considered spatially isolated, meaning it has no spatial connection with other edge points, and is subsequently removed from the edge candidate set. All points that are not removed and are directly connected to other edge candidates in at least one direction are retained and marked as the final boundary points.

[0110] By determining spatial adjacency and eliminating isolated points, false edge points caused by local anomalies, accidental measurement errors, or noisy data can be effectively eliminated, ensuring that the final set of boundary points has spatial continuity and physical meaning. This avoids misinterpretation. In this patented solution and the field of spatial analysis, "adjacency" refers to the direct spatial adjacency between two mesh elements in a two-dimensional mesh array. Specifically, for any candidate edge point, its adjacency is defined by octet connectivity.

[0111] The eight-way connectivity refers to the relationship between a grid cell and its adjacent grid cells in the eight directions above, below, left, right, upper left, upper right, lower left, and lower right (i.e., the four diagonal directions). If another grid cell is located in any one of these eight directions of the target grid and is also identified as an edge candidate point, then the two grid cells are considered adjacent. The criteria for determining the adjacency are: the spatial positions of the two grid cells differ by at most 1 in row and column numbers, and both grid cells belong to the set of edge candidate points. Only when this condition is met are the two edge candidate points considered spatially adjacent.

[0112] Furthermore, for each boundary point, check whether there are other boundary points in its neighborhood in eight directions, and group all boundary points that can be connected to each other through their neighborhood in eight directions into the same connected component. After traversing all boundary points, obtain all connected components.

[0113] For each connected component, if each boundary point in the connected component has at least two other boundary points connected to it in its neighborhood in eight directions, then the connected component forms a closed boundary. If there is a boundary point in the connected component that is connected to only one other boundary point in its neighborhood in eight directions, then the connected component is a non-closed component.

[0114] For all non-closed branches, identify endpoint boundary points, where each endpoint boundary point is a boundary point connected to only one other boundary point. In each non-closed branch, extract all endpoint boundary points.

[0115] For all endpoint boundary points, select the endpoint boundary point with the smallest spatial distance as a pair of endpoint boundary points. For each pair of endpoint boundary points, select the grid on the shortest spatial path between the endpoint boundary points and mark it as a boundary point. Count all boundary points to form a boundary point set.

[0116] For all the selected boundary points, check in turn whether there are other grids marked as boundary points in the neighborhood of each boundary point in its eight directions (i.e., up, down, left, right and four diagonal directions). By checking the spatial adjacency, all boundary points that can be connected to each other through the eight-neighborhood relationship are grouped into the same connected component. A connected component is a group of boundary points that can be continuously connected together in space along the eight-neighborhood relationship. Connected components are not connected to each other. After traversing all boundary points and completing the grouping, all spatial connected components can be obtained.

[0117] Each connected branch is analyzed to determine whether it constitutes a closed boundary. The determination method is as follows: check the number of direct connections between each boundary point in the branch and other boundary points within its eight neighborhood. If all boundary points in the branch are directly connected to at least two other boundary points within their eight neighborhoods, it means that the connected branch forms a closed loop in space, which is a closed boundary. A closed boundary can completely delineate the spatial range of a geothermal anomaly zone.

[0118] If a boundary point in a connected component is directly connected to only one other boundary point within its eight neighborhoods, this indicates that the component has a spatial discontinuity or opening and cannot form a closed region. Such a connected component is considered a non-closed component. For each non-closed component, it is necessary to further identify its endpoint boundary points. The criterion for determining endpoint boundary points is: they are directly connected to only one other boundary point within their eight neighborhoods. By traversing every point in the non-closed component, all boundary points that satisfy the endpoint condition are extracted.

[0119] The extracted endpoint boundary points are paired. For each endpoint boundary point, its spatial distance to other endpoint boundary points within the branch is calculated, and the endpoint boundary point with the smallest distance is selected as its pairing object, forming a set of endpoint boundary point pairs. For each pair of endpoint boundary point pairs, a spatial shortest path search algorithm is used. Along the shortest path in the grid array, from one endpoint to another, grids along the path are selected sequentially, and these grids are marked as boundary points. This completes and connects the original non-closed branch breaks, effectively connecting all disconnected branches that failed to form closed regions, significantly improving the spatial integrity and closure of the geothermal anomaly zone boundary. Finally, all the completed boundary points are statistically analyzed to form the final boundary point set. This set includes not only the original spatially continuous boundary points but also the break points completed by the shortest path algorithm, providing a more realistic and complete reflection of the spatial boundary of the geothermal anomaly zone.

[0120] In non-closed branches, extreme cases may occur, such as an odd number of endpoint boundary points, extremely uneven distance distribution, and multiple endpoints being equidistant from each other. To ensure the rationality and consistency of the completion operation, the coordinates of all endpoint boundary points (i.e., their row and column numbers in the grid array) can be counted. For all endpoint boundary point pairs, the Euclidean or Manhattan distance between each pair is calculated, and the pairing is performed using the "minimum weight matching" principle. That is, the two closest endpoints are selected for pairing first. After pairing, these two endpoints are removed from the set to be paired, and then the search for the closest pair among the remaining endpoints continues, and so on, until all endpoints are paired or only a single endpoint remains. If a single endpoint cannot be paired, it can be manually or supplemented by combining the overall boundary shape. If there are multiple endpoints with equal distances, the endpoint with the smallest row or column number can be selected for pairing first to ensure the uniqueness and repeatability of the pairing.

[0121] In a two-dimensional grid array, the spatial location of each grid can be uniquely identified by its row number and column number. For each pair of endpoint boundary points, the shortest spatial path between the endpoint boundary points refers to taking one endpoint boundary point as the starting point and the other as the ending point. That is, only horizontal, vertical or diagonal movement is allowed. Move from the starting endpoint to the ending endpoint in sequence. At each step, choose to move one step in the direction that is closer to the target endpoint in the row number or column number. Prioritize shortening the distance along the diagonal direction. The shortest spatial path is the path that passes through the fewest grids from the starting point to the ending point. In the shortest spatial path, record the grid number passed at each step. Finally, mark all grids on the path in sequence as the boundary points to be completed.

[0122] In this embodiment, the method for expanding the geothermal distribution area based on boundary points using a region growing algorithm to generate the predicted geothermal distribution area is as follows:

[0123] All boundaries in the set of boundary points are used as starting seed points. For each starting seed point, all neighboring grids that have not yet been marked as boundary points in its eight directions are identified as current expansion candidate grids to be determined. For each expansion candidate grid, its two-dimensional feature vector is extracted, and the Euclidean distance between it and the two-dimensional feature vector of the current starting seed point is calculated.

[0124] If the Euclidean distance is less than the preset growth threshold, the candidate grid is determined to meet the regional growth conditions and can be included in the geothermal distribution area. If the Euclidean distance is greater than the growth threshold, the candidate grid is determined not to be included in the geothermal distribution area for the time being.

[0125] For all candidate grids that meet the growth conditions, they are assigned to the geothermal distribution area and used as new growth seed points. For the newly assigned growth seed points, the above steps of expanding candidate grid determination, calculating the Euclidean distance between two-dimensional feature vectors, and determining the regional growth conditions are repeated recursively until no new grids meet the regional growth conditions during the regional growth process. At this point, the regional growth algorithm terminates. All assigned grids and boundary points in the boundary point set together constitute the predicted geothermal distribution area.

[0126] Specifically, the set of all boundary points identified through the aforementioned process is used as the starting seed points for the region growth algorithm. Each boundary point has a clear position in the grid array and a corresponding two-dimensional feature vector (composed of normalized surface temperature and surface heat flow data). These boundary points serve as the starting basis for growth, ensuring that the expanded region has consistent geothermal characteristics with the boundary of the identified anomalous area. For each starting seed point, its neighboring grids in eight directions are checked, including the top, bottom, left, right, and four diagonal directions. For each neighboring grid, it is determined whether it has not yet been marked as a boundary point and has not been assigned to a geothermal distribution area. Only grids that meet these conditions are considered as candidate grids for expansion to be determined. The purpose of this operation is to ensure that the growth region does not repeatedly cover the already identified boundary points or previously assigned grids, avoiding invalid expansion and infinite loops.

[0127] For each candidate grid to be expanded, its corresponding two-dimensional feature vector is extracted. Then, the Euclidean distance between the feature vector of the candidate grid and the feature vector of the current seed point is calculated. Specifically, the difference between the two in terms of normalized surface temperature and normalized surface heat flux is calculated, squared, added together, and then the square root is taken to obtain the Euclidean distance. This distance value is used to measure the similarity between the candidate grid and the seed point in geothermal physical characteristics. The Euclidean distance calculated in the previous step is compared with a pre-set growth threshold. If the distance is less than the threshold, it means that the candidate grid is highly similar to the seed point in geothermal characteristics and is determined to meet the regional growth conditions and can be included in the geothermal distribution area. At this time, the candidate grid is marked as belonging to the geothermal distribution area, indicating that it has become part of the geothermal distribution area. Conversely, if the distance is greater than the threshold, the candidate grid does not participate in the regional expansion.

[0128] All new grids that meet the growth conditions in this round are used as new growth seed points. The process of neighborhood expansion, feature vector extraction, Euclidean distance calculation, and growth condition determination is repeated for each new seed point. This recursive iterative process continues until no new grids meet the growth conditions, that is, when the neighborhoods of all seed points in the current round cannot be expanded, the algorithm terminates. After the recursive iteration ends, all grids assigned to the geothermal distribution area, together with the initial set of boundary points, constitute the complete predicted geothermal distribution area. This area is not only spatially continuous, but also highly consistent with the known boundary in terms of physical characteristics such as geothermal temperature and heat flow.

[0129] This region growing method ensures that the spatial expansion of geothermal anomaly zones is entirely based on data-driven approaches and physical feature similarity, greatly improving the scientific rigor and accuracy of the predicted regions. The recursive expansion can automatically adapt to the actual spatial distribution of geothermal zones, avoiding errors caused by subjective delineation, and effectively eliminating neighboring grids with inconsistent physical features.

[0130] For each extended candidate grid, the specific method for calculating the Euclidean distance between two-dimensional feature vectors is as follows: For the two-dimensional feature vectors of the extended candidate grid and its corresponding starting seed point, calculate the difference between the two indicators of surface temperature observation and surface heat flow observation, square them respectively, add them together, and then take the square root. This distance reflects the overall similarity between the two grids in the geothermal feature space. The smaller the value, the closer the features are.

[0131] The growth threshold is a key parameter for determining whether an expanded candidate grid can be included in the geothermal distribution area. Its setting method refers to the experience of previous engineering projects or similar geothermal fields, selecting an Euclidean distance value that can effectively distinguish between anomaly areas and background areas. For example, in common normalized feature spaces, the threshold can be set to 0.2, 0.3, etc. As a similarity measure of multidimensional feature space, Euclidean distance can comprehensively reflect the consistency of candidate grids and boundary points in geothermal physical characteristics, avoiding misjudgment due to single-dimensional anomalies. By reasonably setting the growth threshold, the continuity and purity of the geothermal area can be effectively controlled, ensuring the physical consistency of the predicted area and preventing irrelevant grids from being mistakenly included, thereby improving the accuracy of the overall spatial prediction.

[0132] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0133] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0134] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0135] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a division of some logical functions, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or units may be electrical, mechanical, or other forms.

[0136] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0137] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A geothermal energy prediction method based on multi-source data, characterized by the following steps: include: The target area to be predicted is divided according to the preset grid scale. The target area is divided into regular grids and the grids are numbered to generate a two-dimensional grid array. Surface temperature and surface heat flow observations are obtained at the center point of each grid. The surface temperature and surface heat flow observations of all grids are normalized respectively. The normalized surface temperature observations and normalized surface heat flow observations of each grid are combined to form a two-dimensional feature vector that reflects the geothermal characteristics of the grid. For non-boundary grids in a two-dimensional grid array, based on two-dimensional feature vectors and using Euclidean distance, feature vector differences are generated to reflect the degree of feature change between the grid and each adjacent grid. The feature gradient value with the largest feature vector difference is selected as the feature gradient value of the grid. The feature gradient value is used to measure the degree of spatial abrupt change of the grid in geothermal features. The average level and fluctuation range of all feature gradient values ​​are statistically obtained. The linear weighted sum of the average level and fluctuation range is used as the edge identification threshold. Each feature gradient value is compared with the edge identification threshold. If the feature gradient value corresponding to the grid is greater than or equal to the threshold, the grid is determined to be a candidate edge point of the geothermal area. For all grids identified as edge candidates, connectivity analysis is performed based on spatial octet connectivity. Isolated edge candidates are removed to obtain a set of boundary points. Based on the boundary points, the geothermal distribution area is expanded using a region growing algorithm. The two-dimensional feature vector of the grid is used as the criterion to determine whether adjacent grids meet the region growing conditions, thus generating the predicted geothermal distribution area.

2. The geothermal energy prediction method based on multi-source data according to claim 1, characterized in that: The target area is a top-down projection of the surface area to be predicted. The grid is a square grid. The method of dividing the target area into regular grids is based on a preset grid scale, and the space is evenly divided according to the direction of rows and columns. Each grid is generated in sequence, and all grids completely cover the target area to generate a two-dimensional grid array. At the boundary of the target area, if there are grids that exceed the boundary of the target area, the grids that exceed the target area are removed, and only the grids that completely cover the target area are retained. The direction of the rows and columns is based on the geographical direction, with the northernmost grid as the first row and the row number increasing from north to south, and the westernmost grid as the first column and the column number increasing from west to east. When numbering the grids, the northernmost and westernmost grids in the target area are determined as the starting points of the grid array.

3. The geothermal energy prediction method based on multi-source data according to claim 2, characterized in that: When acquiring surface temperature and surface heat flow observations, surface temperature and surface heat flow observations at all grid center points are collected simultaneously. The collection method is as follows: The center point of each grid is determined as the spatial reference for data collection. Surface temperature is observed and recorded at the center point of each grid, and the surface temperature value at that point is taken as the surface temperature observation value of the grid center point. At the same grid center point, the surface heat flow data at that point is continuously recorded, and the average value is taken as the surface heat flow observation value of the grid. The surface temperature observation value and surface heat flow observation value collected at each grid center point are mapped to the grid number.

4. The geothermal energy prediction method based on multi-source data according to claim 1, characterized in that: For all grids, the minimum and maximum values ​​of surface temperature observations and surface heat flow observations are statistically analyzed. The minimum-maximum normalization method is used to standardize the surface temperature and surface heat flow observations of each grid to between 0 and 1. The normalized surface temperature observations and surface heat flow observations are combined to form a two-dimensional vector, denoted as the first vector. Line number The two-dimensional eigenvectors of the column grid are The combination method is as follows: ,in, and They represent Line number Surface temperature and surface heat flow observations after grid normalization.

5. The geothermal energy prediction method based on multi-source data according to claim 4, characterized in that: The non-boundary grid is a grid located inside a two-dimensional grid array, where there are directly adjacent grids in all four directions: up, down, left, and right; that is, a non-boundary grid is not the outermost grid, but is surrounded by other grids and has a complete four-neighbor relationship. The logic for generating the feature vector difference between the grid and its neighboring grids is as follows: For each non-boundary grid, first determine the row and column number of the grid in the array, and then extract the two-dimensional feature vectors of the grid and its four adjacent grids above, below, to the left and to the right in sequence according to the grid number; The two-dimensional feature vector of the target grid is compared with the two-dimensional feature vector of its neighboring grids, and the difference in geothermal feature state between the target grid and its neighboring grids is calculated. The difference value is calculated by calculating the numerical difference between the target grid and its neighboring grids in the two dimensions of normalized surface temperature and normalized surface heat flux, and then taking the square root of the sum of the squares of the differences in the two dimensions to obtain the degree of feature change between the target grid and each neighboring grid, which is labeled as the feature vector difference value. After obtaining the feature vector difference values ​​in four directions for each non-boundary grid, the largest value among the feature vector difference values ​​in the four directions is selected as the feature gradient value of that grid.

6. The geothermal energy prediction method based on multi-source data according to claim 5, characterized in that: The logic for obtaining the edge recognition threshold by statistically analyzing all feature gradient values ​​is as follows: The characteristic gradient values ​​of all non-boundary grids in the two-dimensional grid array are statistically analyzed, and the average level and standard deviation of all characteristic gradient values ​​are calculated. The average level and fluctuation amplitude of the characteristic gradient values ​​are combined in a linear weighted manner to obtain the identification threshold for judging the edge of the geothermal area. The feature gradient value of each grid is compared with the edge recognition threshold. For grids whose feature gradient value is greater than or equal to the edge recognition threshold, the grid is determined to be a candidate edge point of the geothermal area.

7. The geothermal energy prediction method based on multi-source data according to claim 1, characterized in that: When removing spatially isolated edge candidate points to obtain boundary points, all grids that are identified as edge candidate points are traversed. For each edge candidate point, its spatial position in the two-dimensional grid array is checked, and its adjacency relationship with other edge candidate points is identified. For each edge candidate point, analyze its neighborhood in eight directions, namely the positions of the eight adjacent grids: above, below, left, right, and upper left, upper right, lower left, and lower right, and determine whether there are other grids in these neighborhood positions that are identified as edge candidate points. If no other edge candidate points are found in the neighborhood of a certain edge candidate point in any of its eight directions, then the point is considered to be spatially isolated. Isolated edge candidate points are eliminated, and the remaining edge candidate points are marked as boundary points.

8. The geothermal energy prediction method based on multi-source data according to claim 7, characterized in that: For each boundary point, check if there are other boundary points in its neighborhood in eight directions. Group all boundary points that can be connected to each other through their neighborhoods in eight directions into the same connected component. After traversing all boundary points, obtain all connected components. For each connected component, if each boundary point in the connected component has at least two other boundary points connected to it in its neighborhood in eight directions, then the connected component forms a closed boundary. If there is a boundary point in the connected component that is connected to only one other boundary point in its neighborhood in eight directions, then the connected component is a non-closed component. For all non-closed branches, identify endpoint boundary points, where each endpoint boundary point is a boundary point connected to only one other boundary point. In each non-closed branch, extract all endpoint boundary points. For all endpoint boundary points, select the endpoint boundary point with the smallest spatial distance as a pair of endpoint boundary points. For each pair of endpoint boundary points, select the grid on the shortest spatial path between the endpoint boundary points and mark it as a boundary point. Count all boundary points to form a boundary point set.

9. The geothermal energy prediction method based on multi-source data according to claim 8, characterized in that: The method for expanding the geothermal distribution area based on boundary points using a region growing algorithm to generate the predicted geothermal distribution area is as follows: All boundaries in the set of boundary points are used as starting seed points. For each starting seed point, all neighboring grids that have not yet been marked as boundary points in its eight directions are identified as current expansion candidate grids to be determined. For each expansion candidate grid, its two-dimensional feature vector is extracted, and the Euclidean distance between it and the two-dimensional feature vector of the current starting seed point is calculated. If the Euclidean distance is less than the preset growth threshold, the candidate grid is determined to meet the regional growth conditions and can be included in the geothermal distribution area. If the Euclidean distance is greater than the growth threshold, the candidate grid is determined not to be included in the geothermal distribution area for the time being. For all candidate grids that meet the growth conditions, they are assigned to the geothermal distribution area and used as new growth seed points. For the newly assigned growth seed points, the above steps of expanding candidate grid determination, calculating the Euclidean distance between two-dimensional feature vectors, and determining the regional growth conditions are repeated recursively until no new grids meet the regional growth conditions during the regional growth process. At this point, the regional growth algorithm terminates. All assigned grids and boundary points in the boundary point set together constitute the predicted geothermal distribution area.