A method for determining geothermal anomaly area based on whole rock geochemical data

Abstract

The application relates to the technical field of geothermal exploration information, and discloses a method for judging a geothermal abnormal area based on whole-rock geochemical data, which comprises the following steps: obtaining a whole-rock geochemical original data set, performing cleaning and mapping, and forming a discrete space data set; receiving the discrete space data set, calculating a geothermal activity indicating scalar and a characteristic element concentration ratio variable, synthesizing a comprehensive fluid dynamics zoning index to generate a structured data table; receiving the structured data table, calculating a comprehensive geothermal anomaly index through principal component analysis, converting the comprehensive fluid dynamics zoning index into a spatial gradient vector field; receiving the comprehensive geothermal anomaly index and the spatial gradient vector field, constructing a continuous scalar field to define a multi-level abnormal boundary, combining a divergence value to extract a hydrothermal upwelling dynamics center; combining multi-source prior geological constraint data to perform intersection judgment, delineating a target target area boundary set and outputting a result. The application can quickly complete large-area regional scanning and screening, and improves the accuracy of geothermal exploration selection.

Description

Technical Field

[0001] This invention relates to the field of geothermal exploration information technology, specifically a method for identifying geothermal anomaly zones based on whole-rock geochemical data. Background Technology

[0002] Delineating geothermal anomaly zones in the initial stages of geothermal resource exploration is crucial for reducing exploration risks and improving prospecting efficiency. Traditional evaluation methods rely on geophysical exploration and shallow geothermal measurements. Traditional physical exploration methods are costly, yield results with many interpretations, or only reflect instantaneous surface temperature conditions, failing to provide clear indications of deep reservoir distribution and hydrothermal activity intensity. Whole-rock geochemical survey data contains rich information on material composition and evolution; however, current processing methods lack methods to quantify water-rock interactions, preventing existing regional survey data from being transformed into evaluation indicators with clear geological significance and thus failing to replace physical exploration in reducing initial operational costs.

[0003] Current analyses of geochemical data are mostly limited to discussions of absolute values ​​of single-element concentrations or simple concentration quotients, failing to transform discrete material composition data into continuous spatial mathematical models. Interference from local geological environments can lead to false geochemical anomalies for single elements, resulting in multiple interpretations of the analysis results. Currently, there is a lack of a fusion calculation model of multi-dimensional scalar fields and vector fields reflecting the direction of fluid migration, making it impossible to conduct spatial topological comparisons and difficult to objectively and quantitatively locate hydrothermal upwelling dynamic centers. The objectivity of anomaly center identification is relatively low.

[0004] Current geothermal exploration site selection methods have not yet established a multi-source verification process that combines prior geological constraint data. In the absence of deep and large fault structures or magmatic heat sources, relying solely on compositional anomaly data can easily lead to the delineation of incorrect target area boundaries. Due to the lack of standardized processing steps that sequentially and coherently execute data reading, preprocessing, spatial modeling, and boundary output, rapid large-area area scanning and screening cannot be carried out when faced with massive amounts of basic survey data, thus limiting the data processing efficiency and comprehensive prediction accuracy of geothermal exploration site selection. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for identifying geothermal anomaly zones based on whole-rock geochemical data. This method solves the problems of high operating costs in traditional geothermal exploration methods, multiple solutions in single-element geochemical analysis, and low accuracy in geothermal target area delineation and data processing efficiency due to the lack of multi-source spatial data constraints for verification.

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

[0007] This invention provides a method for identifying geothermal anomaly zones based on whole-rock geochemical data, comprising the following steps:

[0008] The raw whole-rock geochemical dataset was obtained, and data cleaning and spatial mapping were performed to form a standardized discrete spatial dataset.

[0009] Receive the discrete spatial dataset, calculate the geothermal activity indicator scalar and characteristic element concentration ratio variables, synthesize the comprehensive fluid dynamics zonation index and generate a structured data table;

[0010] Receive the structured data table, calculate the comprehensive geothermal anomaly index through principal component analysis, and convert the comprehensive hydrodynamic zonation index into a spatial gradient vector field;

[0011] Receive the comprehensive geothermal anomaly index and the spatial gradient vector field, construct a continuous scalar field and define multi-level anomaly boundaries, and extract the hydrothermal upwelling dynamic center by combining the divergence value of the spatial gradient vector field.

[0012] The system receives the hydrothermal upwelling dynamic center, performs spatial intersection judgment based on the imported multi-source prior geological constraint data, delineates the target area boundary set, and outputs comprehensive geological interpretation results.

[0013] Preferably, in the process of forming a standardized discrete spatial dataset, the preset instrument detection limit parameters are called to compare the concentration values ​​of each element in the original whole-rock geochemical dataset and perform the assignment and transformation operation. An outlier screening mechanism is activated to clear invalid values ​​and complete the data cleaning process. For any missing records, spatial interpolation of missing values ​​is performed, and the standard fraction algorithm is used to implement dimensionless transformation. The corresponding latitude and longitude information is collected and converted into planar projection coordinates, and the planar projection coordinates are associated and bound with the standardized concentration vector.

[0014] Preferably, when calculating the ratio variable of geothermal activity indicator scalars to characteristic element concentrations, the original concentration values ​​of major elements are extracted to calculate the hydrothermal alteration intensity index, and the quotient values ​​of the original concentration values ​​of potassium oxide and sodium oxide are used to generate the alkaline migration index; the enrichment coefficient is calculated by comparing the original concentration values ​​of the target indicator element with the background concentration of the region; the original concentration values ​​of deep high-temperature phase indicator elements are extracted and combined with the original concentration values ​​of shallow low-temperature phase indicator elements to perform quotient calculations to generate the characteristic element concentration ratio variable.

[0015] Preferably, in the stage of synthesizing the comprehensive hydrodynamic zoning index, principal component analysis algorithm is introduced to extract the weight coefficients of feature element pairs; weighted average processing is performed on the concentration ratio variables of multiple feature elements with similar physical meaning to obtain the comprehensive hydrodynamic zoning index; and the generated geothermal activity indicator scalar is packaged and converted together with the comprehensive hydrodynamic zoning index.

[0016] Preferably, in the process of obtaining the comprehensive geothermal anomaly index, the hydrothermal alteration intensity index, alkali migration index, and enrichment coefficient are extracted to establish an initial index matrix. Standardization processing is performed on the column vectors to obtain a standardized index matrix. The correlation coefficient matrix covered by the standardized index matrix is ​​calculated and eigenvalue decomposition is performed to obtain the corresponding eigenvalues ​​and eigenvectors. The eigenvectors corresponding to the first principal component are multiplied with the standardized index matrix to extract the score values ​​of each sampling point in the first principal component dimension as the comprehensive geothermal anomaly index.

[0017] Preferably, in the transformation of the spatial gradient vector field, the ordinary Kriging interpolation algorithm is used to perform gridding processing on the discretely distributed comprehensive fluid dynamics zonation index to obtain the interpolation value; the central difference algorithm is used to solve for the partial derivatives of the horizontal axis projection coordinates and the vertical axis projection coordinates to obtain the gradient components; the complete gradient vectors corresponding to each grid node are synthesized, and the negative gradient direction is selected as the fluid transport direction vector and the magnitude is calculated.

[0018] Preferably, when defining the boundaries of multi-level anomalies, the anisotropic kriging interpolation algorithm is invoked, and the azimuth parameter of the main structural line is set as the guiding parameter to establish a continuous scalar field; the comprehensive geothermal anomaly index contained in all effective sampling points is extracted and the arithmetic mean and standard deviation are calculated, and lower limit thresholds for weak anomalies, obvious anomalies, and core anomalies are set; the continuous scalar field is segmented at the pixel level using multi-level anomaly thresholds, and the outer envelopes of adjacent nodes at the same level are extracted to generate multi-level anomaly boundaries.

[0019] Preferably, the process of extracting the hydrothermal upwelling dynamic center includes: using the central difference algorithm to calculate the divergence value of each grid node; selecting divergence values ​​greater than zero and calculating the arithmetic mean and standard deviation of the positive subset to establish a divergence identification threshold; selecting nodes with divergence values ​​not less than the divergence identification threshold and marking them as strongly divergent nodes; performing a scalar-vector joint spatial topology comparison operation; selecting a set of grid nodes that simultaneously possess both core anomalous node attributes and strongly divergent node attributes and combining them into a continuous spatial polygon.

[0020] Preferably, the spatial intersection judgment step is performed, and unified spatial registration and grid alignment are carried out on the distribution map of deep fault surfaces, the map of magmatic rock outcrops, and remote sensing surface temperature anomaly data; the fault central axis is extracted to generate a buffer polygon as a structural constraint surface, the outer envelope of magmatic rocks is extracted to generate a heat source constraint surface, and a set of continuous pixels with temperature values ​​higher than the background temperature threshold is selected to be converted into a surface thermal anomaly surface; the spatial Boolean intersection algorithm is invoked to determine that the target area is located inside the deep geochemical dynamic anomaly zone and overlaps with the structural constraint surface, heat source constraint surface, or surface thermal anomaly surface in space, and an initial intersection polygon is generated.

[0021] Preferably, before outputting the comprehensive geological interpretation results, the area threshold is removed from the initial intersection polygons and smoothed by the same curve to generate the final geothermal target area polygons as the target target area boundary set; the internal attribute information of the final geothermal target area polygons is extracted to implement the logical inference of the geothermal system's genetic type; and a comprehensive geological interpretation report covering the coordinate parameters of the final geothermal target area polygons and the inferred genetic type is output.

[0022] This invention provides a method for identifying geothermal anomaly zones based on whole-rock geochemical data. It has the following beneficial effects:

[0023] 1. This invention extracts whole-rock geochemical raw data and calculates hydrothermal alteration intensity index, alkali migration index and characteristic element concentration ratio variables to quantitatively express the physicochemical process of water-rock interaction in the geothermal system. Since the basic data required for the calculation comes directly from existing regional survey data, no additional investment in physical exploration is required, effectively controlling the early operation costs. At the same time, the extracted indicators directly correspond to the objective evolution process of fluid alteration of the surrounding rock, making the analysis results have clear geothermal anomaly indication significance.

[0024] 2. This invention utilizes principal component analysis to synthesize a comprehensive geothermal anomaly index scalar field, and employs a central difference algorithm to generate a spatial gradient vector field reflecting the direction of fluid migration. Then, it extracts the core anomaly nodes of the scalar field and performs spatial topological comparison with the strongly divergent nodes of the vector field, achieving quantitative localization of the hydrothermal upwelling dynamic center. The introduction of principal component analysis and spatial divergence calculation transforms discrete material composition data into a continuous spatial mathematical model, eliminating the multiple solutions caused by local environmental interference of a single element and improving the objectivity of anomaly center identification.

[0025] 3. After obtaining the hydrothermal upwelling dynamic center, this invention introduces the influence zone of deep faults, the outcropping range of magmatic rocks, and remote sensing surface temperature anomaly data. By performing multi-source constraint verification through spatial Boolean intersection algorithm and superimposing prior geological entity elements, it can accurately eliminate geochemical false anomalies that lack tectonic or heat source support. The sequential and coherent execution of data reading, preprocessing, spatial modeling, and target area boundary output steps can quickly complete the scanning and screening of large areas, simultaneously improving the accuracy of geothermal exploration area selection and data processing efficiency. Attached Figure Description

[0026] Figure 1 This is a structural block diagram of a system for identifying geothermal anomaly zones based on whole-rock geochemical data according to the present invention.

[0027] Figure 2 This is a flowchart illustrating the overall steps of a method for identifying geothermal anomaly zones based on whole-rock geochemical data according to the present invention.

[0028] Figure 3 This is a graph showing the relationship between the hydrothermal alteration intensity index and the concentration ratio of characteristic elements in this invention and the distance to the geothermal anomaly dynamic center. Detailed Implementation

[0029] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0030] See attached document Figure 1 This invention provides a system for identifying geothermal anomaly zones based on whole-rock geochemical data, including a data integration and preprocessing module, an index calculation engine, a comprehensive modeling module, a spatial analysis and visualization module, and a result output and target area delineation module.

[0031] The data integration and preprocessing module receives regional whole-rock geochemical survey data, performs quality checks and missing value imputation on the raw data containing major and trace elements, eliminates dimensional differences in the concentration data of each element through a standardization algorithm, and then binds the attribute data to the geographic coordinate system to generate a discrete spatial dataset and transmits it to the index calculation engine.

[0032] The index calculation engine receives discrete spatial datasets and calculates the hydrothermal alteration intensity index, alkaline migration index, and volatile element enrichment coefficient corresponding to each spatial coordinate point based on the water-rock reaction mechanism. At the same time, it extracts characteristic element pairs that indicate different evolutionary phases to construct element ratio functions, and then sends various scalar index data and element ratio functions to the integrated modeling module.

[0033] The integrated modeling module receives scalar index data and element ratio functions, uses principal component analysis to reduce the dimensionality of various scalar index data to generate a comprehensive geothermal anomaly index scalar field, and solves partial derivatives in two-dimensional space based on the continuous grid generated by interpolation of the element ratio functions to generate a spatial gradient vector field reflecting the direction of fluid migration. The comprehensive geothermal anomaly index scalar field and the spatial gradient vector field are then transmitted to the spatial analysis and visualization module.

[0034] The spatial analysis and visualization module receives the scalar field of the comprehensive geothermal anomaly index and the spatial gradient vector field, executes a geostatistical interpolation algorithm to transform the discrete comprehensive geothermal anomaly index into a continuous grid surface, delineates geothermal anomaly zones at different levels according to the set statistical distribution threshold, and transforms the spatial gradient vector field into an array of directional arrows superimposed on the grid surface to generate a scalar vector joint topological image and outputs it to the results output and target area delineation module.

[0035] The results output and target area delineation module receives scalar vector joint topological images, imports prior deep and large fault surface region data, magmatic rock outcrop range maps, and remote sensing surface temperature anomaly data, performs intersection analysis of multi-source spatial data based on Boolean logic operations, delineates the boundaries of geothermal exploration target areas that meet the constraints, and outputs target area coordinate parameters.

[0036] See attached document Figure 2 This invention provides a method for identifying geothermal anomaly zones based on whole-rock geochemical data, which, combined with the above-mentioned system architecture, includes the following steps:

[0037] S1. The data integration and preprocessing module acquires systematic whole-rock geochemical data of the study area. This data includes a subset of major elements and a subset of trace elements in oxide form. The data integration and preprocessing module performs data cleaning and standardization, uniformly assigns spatial coordinate parameters, and forms a standardized spatial dataset.

[0038] S2. The index calculation engine extracts geothermal activity indicative indicators at each discrete sampling point location, calculates the hydrothermal alteration intensity index, alkali migration index, and enrichment coefficient of target trace elements, establishes characteristic element pairs indicating deep high-temperature phases and shallow low-temperature phases, and generates corresponding concentration ratio variables.

[0039] S3, the integrated modeling module constructs an integrated scalar field and a spatial vector field for the discrete sampling point locations. It extracts the variance contribution weight of the scalar index through principal component analysis and calculates the integrated geothermal anomaly index. It calls the ordinary Kriging interpolation algorithm to fit the continuous surface and establishes the fluid transport vector components at each coordinate position by solving the directional gradient.

[0040] The S4 spatial analysis and visualization module performs anisotropic interpolation calculations to generate a continuous scalar field of the comprehensive geothermal anomaly index. It defines multi-level anomaly boundaries based on statistical thresholds and extracts the hydrothermal upwelling dynamic center that coincides with the core anomaly node and the strong divergence node.

[0041] S5. The output of results and the target area delineation module superimpose data on the structural fault influence zone, magmatic heat source zone and surface thermal anomaly surface, perform spatial Boolean intersection judgment under multi-dimensional constraints, and define the target area boundary set located around the main heat-controlling fault zone or above the heat source and satisfying geochemical dynamics.

[0042] The data integration and preprocessing module is equipped with a data interface for reading raw whole-rock geochemical datasets from external storage media. These datasets record the spatial coordinates and corresponding chemical element test values ​​of multiple discrete geological sampling points. Based on predefined program instructions, the data integration and preprocessing module converts external non-standard format data files into a standardized array that the system can process, and then executes the following steps.

[0043] S101, the data integration and preprocessing module reads the elemental attribute fields from the original dataset and divides the extracted concentration values ​​into major element subsets and trace element subsets according to their geochemical properties. The major element subset includes rock-forming components existing in the form of oxides, specifically covering silicon dioxide, aluminum oxide, potassium oxide, sodium oxide, calcium oxide, magnesium oxide, and ferric oxide.

[0044] The trace element subset contains trace components that record the characteristics of hydrothermal fluid activity, specifically including lithium, rubidium, cesium, strontium, barium, fluorine, chlorine, arsenic, antimony, and boron. The extraction of these specific elements provides a fundamental data source for the subsequent construction of hydrothermal alteration intensity indices and alkaline element migration indices.

[0045] In step S102, the data integration and preprocessing module traverses the major element subset and the trace element subset, performing element concentration detection limit verification. Due to limitations in the physical detection capabilities of laboratory analytical instruments, the actual concentration of trace elements at some sampling points is lower than the effective measurement limit. Therefore, the system calls the preset instrument detection limit parameters for each element in the database and compares the element concentration values ​​extracted from the sampling points one by one.

[0046] When a specific element concentration at a particular sampling point is found to be below the corresponding instrument detection limit, half of the instrument detection limit is extracted and assigned to the concentration variable of that specific element at that sampling point. This assignment operation completes the numerical conversion of data below the detection limit, preventing system errors such as division by zero overflow during subsequent index ratio calculations.

[0047] S103, the data integration and preprocessing module initiates an outlier screening mechanism to eliminate extreme data bias caused by external non-geological factors. Specifically, it extracts the concentration set of each element across all sampling points in the study area, calculates the statistical median and absolute median difference of this concentration set, and uses these two statistical parameters to define the effective distribution range of each element.

[0048] The lower limit of the effective distribution interval is set as the statistical median minus the product of the threshold multiple and the absolute median difference, and the upper limit is set as the statistical median plus the product of the threshold multiple and the absolute median difference. This threshold multiple controls the strictness of outlier removal, and its value is set between 2 and 3; for example, a value of 2.5 obtained from the system configuration file. Element concentration values ​​outside this effective distribution interval are identified as outliers and invalid, and are cleared.

[0049] S104, the data integration and preprocessing module performs missing value imputation for missing records generated after clearing outliers or attribute fields omitted in the original data acquisition. Since geological fluid activity exhibits gradual changes in physical space, a spatial distance weighting algorithm is performed based on the planar spatial coordinate parameters of the sampling points. Using the target sampling point with missing data as the central origin, the geometric straight-line distance between this central origin and surrounding valid sampling points with complete attribute parameters is calculated.

[0050] Following the geometric straight-line distance from nearest to farthest, a preset number of valid sampling points are retrieved. This preset number defines the spatial neighborhood range for the interpolation calculation, and its value is set between 4 and 12; for example, based on the overall sampling grid density of the study area, the value is 8. Subsequently, using the reciprocal of each calculated geometric straight-line distance as a weighting coefficient, the element concentrations corresponding to the obtained valid sampling points are weighted and summed. The result of this weighted summation is then divided by the sum of all weighting coefficients involved in the calculation, and the final normalized result is used as the interpolation value and written to the missing location of the target sampling point. This interpolation logic ensures the integrity of the geochemical dataset, establishing a data foundation for the subsequent construction of a continuous and comprehensive geological anomaly spatial field.

[0051] After completing the verification and missing value imputation of the preliminary data, the data integration and preprocessing module further performs dimensionless processing and spatial attribute association operations to address the structural differences in the numerical values ​​themselves, and specifically performs the following steps.

[0052] S105, the data integration and preprocessing module performs dimensionless processing of element concentrations. Major element concentrations in geological samples are typically measured as percentages, while trace element concentrations are often on the order of parts per million. The absolute magnitude difference between the two can reach tens of thousands of times. To avoid the larger absolute values ​​of major elements masking the anomalous changes in trace elements during subsequent calculations of the comprehensive geological anomaly index, the system must eliminate the influence of inherent dimensions.

[0053] Specifically, the standard score algorithm is used to transform the data of all target elements. This process is achieved through the following formula:

[0054] ;

[0055] in, Indicates the first At the sampling point, the first Standardized values ​​of elements; This indicates the original concentration value at this location after preliminary processing; Indicates the first The arithmetic mean of the elements in the entire valid sample set of the study area; Indicates the first The standard deviation of a seed element across the entire sample set; The range of values ​​is ,in This represents the total number of valid sampling points within the study area; The range of values ​​is ,in This represents the total number of major and minor elements involved in the calculation. Based on this calculation logic, the numerical distributions of various elements are uniformly shifted and scaled to a standard space with a mean of 0 and a standard deviation of 1, effectively eliminating statistical weight bias caused by different units.

[0056] S106, the data integration and preprocessing module performs a two-dimensional discretization mapping of attribute data and geospatial coordinates. The evolution analysis of geological fluids is highly dependent on the relative positions of materials in physical space. The system reads the location field from the original record file, obtains the latitude and longitude information corresponding to each sampling point, and converts it uniformly into planar projected coordinates. For the mathematical operations of converting latitude and longitude coordinates to a specific planar projected coordinate system, those skilled in the art can use conventional mapping methods such as the Gauss-Kruger projection; the specific calculation of projection parameters is well-known in the field and will not be elaborated here.

[0057] After completing the coordinate system transformation, the system extracts the first... Planar spatial coordinate parameters of each sampling point .in Represents the horizontal axis projection coordinates. This represents the projected coordinates of the vertical axis. The system associates and binds these spatial coordinate parameters one-to-one with the standardized concentration vector generated in the previous step, thereby constructing a discrete spatial dataset for the entire region.

[0058] This discrete space dataset is mathematically defined as .in, This represents the total number of valid sampling points within the study area. Representing the The concentration vector is constructed by normalizing the values ​​of all target elements contained in each sampling point. By the At the sampling point, the first Standardized values ​​of elements Composed of elements in order, i.e., the mathematical expression is: Through the above binding operation, chemical property data is associated with spatial coordinates, and the generated discrete spatial dataset is transmitted to the index calculation engine as input data for constructing the fluid evolution model.

[0059] After receiving a discrete spatial dataset containing raw concentration values ​​and standardized values, the index calculation engine extracts scalar feature parameters indicating geothermal activity at each discrete coordinate point based on the material exchange mechanism between geological fluids and surrounding rocks. The specific calculation steps are as follows.

[0060] S201. Extract major element data and calculate the hydrothermal alteration intensity index corresponding to each spatial coordinate point. During upward migration, geothermal fluids undergo metasomatic reactions with the surrounding rocks, leading to the loss of soluble alkaline oxides and the relative enrichment of ferromagnesian components. The calculation of this index relies on the original concentration values ​​in the major element subsets, rather than dimensionless data, to maintain its inherent rock chemical and physical meaning. The calculation formula is as follows:

[0061] ;

[0062] In this formula, Indicates the first Hydrothermal alteration intensity index at each sampling point; , , and Representing the first The original concentration values ​​of potassium oxide, magnesium oxide, sodium oxide, and calcium oxide at each sampling point. This index reflects the intensity of water-rock reaction in a local space by quantitatively assessing the relative proportion of feldspar mineral destruction to alteration mineral formation.

[0063] S202 calculates an alkaline migration index used to characterize deep fluid channels. Deep, high-temperature fluids in geothermal systems are typically rich in potassium. As the fluid rises along fault lines and undergoes cooling and decompression, potassium ions in the fluid undergo extensive cation exchange reactions with sodium ions in the surrounding rocks. Therefore, the potassium-to-sodium concentration ratio constitutes a direct scalar for tracing hydrothermal migration paths.

[0064] The original concentration values ​​of potassium oxide and sodium oxide at each sampling point are extracted, and an alkaline migration index is generated by calculating their ratio. Mathematically, this index is defined as the first... At each sampling point and The direct quotient of this ratio. Areas with high spatial values ​​for this ratio typically correspond to major thermally conductive fracture zones or central channels where deep fluids surge strongly.

[0065] S203, the index calculation engine calculates the volatility and enrichment coefficients of ore-forming elements indicating shallow fluid activity. Deep geothermal fluids contain a large number of volatile or easily migratable trace elements dissolved in them. When the fluid migrates to near the surface, these elements are unloaded and precipitated in shallow tectonic fractures due to changes in the physicochemical environment. The system extracts lithium, fluorine, chlorine, and boron as target indicator elements from a subset of trace elements. For any of these target indicator elements, its relative enrichment at local spatial points is calculated.

[0066] The enrichment coefficient is calculated as the ratio of the original concentration of the target element at a specific sampling point to the regional background concentration. The system reads the statistical median of the original concentration of the target indicator element at all valid sampling points throughout the study area, sets it as the regional background concentration, and performs a division operation. When the calculated enrichment coefficient is greater than a preset enrichment threshold, it confirms the presence of superimposed injection of exogenous hydrothermal fluids at that specific location. This enrichment threshold is used to define the lower limit of elemental anomalies, and its value is set between 1.5 and 2.0. The system reads the specific value of this threshold based on the overall geochemical variability coefficient of the survey area, for example, a value of 1.5. The magnitude of the enrichment coefficient directly records the intensity of shallow hydrothermal activity.

[0067] S204, the index calculation engine calculates geothermal temperature scale parameters used to indirectly estimate the state of deep geothermal reservoirs. The specific element concentration ratios resulting from water-rock reactions are strictly controlled by the equilibrium temperature of the deep geothermal reservoir. The index calculation engine extracts the original concentration values ​​from the major element subsets and constructs an element geochemical temperature scale assemblage. Based on the lithological characteristics of the study area, the system systematically extracts the original concentration values ​​of sodium oxide, potassium oxide, and calcium oxide at each sampling point and converts them into the corresponding molar concentrations of sodium, potassium, and calcium ions according to their stoichiometric mass.

[0068] Subsequently, the system calculates the sodium-potassium-calcium geothermal temperature scale parameters using the converted molar concentrations. For the specific logarithmic function calculation equations and the selection of empirical coefficients for these parameters, those skilled in the art can utilize conventional geothermal empirical temperature scale formulas; the specific mathematical calculation process is well-known in the field and will not be elaborated upon here. These geothermal temperature scale parameters, serving as independent scalars to assist in verifying the temperature magnitude of deep heat sources, are stored in the attribute set of each discrete coordinate point for subsequent verification of comprehensive anomaly characteristics.

[0069] After completing the calculation of the basic scalar characteristic parameters, the index calculation engine extracts characteristic element pairs indicating different evolutionary phases based on the differences in geochemical behavior of geological fluids under temperature and pressure gradient changes, and then performs the following steps.

[0070] S205 extracts characteristic elements indicating the deep high-temperature phase and shallow low-temperature phase of geothermal fluids. During the migration of geothermal fluids from deep heat sources to shallower areas, the solubility of different elements varies with changes in the physicochemical environment. The index calculation engine divides the trace element subset into deep high-temperature phase indicator elements and shallow low-temperature phase indicator elements. In a specific configuration, rubidium and cesium are extracted as deep high-temperature phase indicator elements, and arsenic and antimony are extracted as shallow low-temperature phase indicator elements. Rubidium and cesium tend to exist stably in high-temperature fluids and accumulate in the center of the heat source or deep channels, while arsenic and antimony have strong volatility and migration capabilities, tending to unload and precipitate at the far end of fluid migration or in shallow low-temperature fractures. The combination of these two types of elements constitutes the basic combination for tracking the spatial migration trajectory of fluids.

[0071] S206, constructing an element ratio function reflecting the fluid evolution gradient based on extracted feature elements. The index calculation engine extracts the original concentration values ​​of feature elements at each spatial coordinate point and calculates the quotient of the indicator element in the deep high-temperature phase and the indicator element in the shallow low-temperature phase. Taking the feature element pair composed of rubidium and arsenic as an example, the calculation formula is as follows:

[0072] ;

[0073] in, Indicates the first The concentration ratio of rubidium to arsenic at each sampling point; Indicates the first The original concentration values ​​of rubidium at each sampling point; Indicates the first The original arsenic concentration values ​​at each sampling point. Near the deep fluid channel, the concentration of rubidium is relatively high while that of arsenic is relatively low, and this ratio shows a spatially high value. As the fluid migrates to shallower areas, due to the premature crystallization and precipitation of elements such as rubidium, while arsenic continues to migrate and accumulate with the fluid, this ratio shows a spatially decreasing trend. This spatial gradient change in numerical value constitutes the mathematical basis for determining the direction of fluid migration.

[0074] S207 expands the dimension of the element ratio function and calculates the comprehensive hydrodynamic zoning index. A single element ratio may be affected by lithological differences or localized mineralization inhomogeneities in a local space. The index calculation engine performs the same ratio calculation on multiple preset characteristic element pairs, forming an element ratio set. The index calculation engine performs a weighted average of multiple ratio variables with similar physical meanings to calculate the comprehensive hydrodynamic zoning index. The calculation formula is as follows:

[0075] ;

[0076] in, Indicates the first Comprehensive fluid dynamics zoning index at each sampling point; This represents the total number of extracted feature element pairs, and its value is set between 2 and 5. For example, when configuring 3 element pairs (rubidium / arsenic, rubidium / antimony, and cesium / bismuth) to participate in the calculation in the system, The value is 3; Indicates the first At the sampling point, the first The concentration ratio of a pair of characteristic elements; Indicates the first The weight coefficients of the feature element pairs.

[0077] For this weighting coefficient The specific values ​​are obtained by the indicator calculation engine through principal component analysis algorithm. For all... The first principal component is extracted from the data matrix of each ratio variable. The eigenvector elements corresponding to this first principal component are then normalized, and the resulting values ​​are set as the weight coefficients for each pair of eigenvector elements. This weight allocation logic objectively assigns values ​​based on the degree of common variation of each ratio variable in its spatial distribution, reducing the bias caused by relying on manual experience to set weights. For the specific mathematical operations of constructing the covariance matrix and solving the eigenvectors in principal component analysis, those skilled in the art can use conventional linear algebra calculation programs. The matrix dimensionality reduction process is a well-known technique in this field and will not be elaborated upon here.

[0078] After the calculations are completed, the index calculation engine packages all the generated geothermal activity scalar indicators and integrated hydrodynamic zonation indices into a structured data table. This structured data table is then sent to the integrated modeling module to provide variable inputs for constructing the integrated geothermal anomaly index scalar field and solving the spatial gradient vector field.

[0079] After receiving the structured data table sent by the index calculation engine, the integrated modeling module uses principal component analysis to construct a comprehensive geothermal anomaly index in order to eliminate multicollinearity among multi-source scalar indices and extract the core features reflecting geothermal activity. The specific steps are as follows.

[0080] S301, the integrated modeling module extracts various geological indicator scalars from the structured data table. The specific indicators involved in the calculation include the hydrothermal alteration intensity index, alkali migration index, and volatile element enrichment coefficient generated in the previous steps. The integrated modeling module integrates these indicators according to the spatial sampling point order to construct an initial indicator matrix for the study area.

[0081] The total number of scalar indices participating in dimensionality reduction calculation is The total number of valid sampling points in the study area is The dimension of the initial index matrix constructed is . To eliminate the inherent differences in physical meaning and numerical magnitude between different scalar indices, the integrated modeling module performs column-wise standardization on the initial index matrix to obtain a standardized index matrix.

[0082] S302, the integrated modeling module calculates the correlation coefficient matrix of the standardized index matrix. This correlation coefficient matrix is ​​used to quantify the degree of linear correlation between the various scalar indicators in their spatial distribution. Subsequently, the integrated modeling module performs eigenvalue decomposition on the correlation coefficient matrix to obtain the corresponding eigenvalues ​​and eigenvectors.

[0083] For the calculation of eigenvalues ​​and eigenvectors of the correlation coefficient matrix, those skilled in the art can use conventional linear algebra calculation programs. The matrix decomposition and feature extraction algorithms are well-known technologies in this field and will not be described in detail here.

[0084] S303, the integrated modeling module determines the extracted principal components based on the calculated eigenvalues. The integrated modeling module sorts the acquired eigenvalues ​​in descending order and calculates the variance contribution rate of the first principal component using the ratio of the eigenvalue of the first principal component to the sum of all eigenvalues.

[0085] To ensure that a single principal component can represent the core characteristics of the comprehensive geothermal anomaly and filter out redundant noise, the system sets a principal component validity threshold. This threshold is set to 40%. The comprehensive modeling module determines whether the variance contribution rate of the calculated first principal component exceeds this validity threshold. If the determination is yes, the system retains the first principal component and abandons the extraction of subsequent minor components, thereby establishing the spatial variation component of the dominant anomaly.

[0086] S304, the integrated modeling module calculates the comprehensive geothermal anomaly index at each discrete spatial coordinate point. After establishing that the first principal component satisfies the variance constraint, the integrated modeling module multiplies the eigenvector corresponding to the first principal component with the standardized index matrix to complete the linear projection of the basic data. The system directly extracts the score value of each sampling point on the first principal component and defines it as the comprehensive geothermal anomaly index. Its calculation formula is as follows:

[0087] ;

[0088] in, Indicates the first The comprehensive geothermal anomaly index at each sampling point; This represents the total number of categories of scalar indicators involved in the dimensionality reduction calculation; Indicates to from to All calculated items are summed. This represents the eigenvector corresponding to the first principal component. Weighting coefficients for each scalar indicator; Indicates the first At the sampling point, the first The standardized value of a scalar indicator.

[0089] This computational logic integrates multiple single geological indicators with physical and chemical correlations into a dimensionless comprehensive scalar parameter, establishing a fundamental scalar field characterizing the spatial intensity of geothermal anomalies.

[0090] After obtaining the comprehensive fluid dynamics zonation index of discrete points, the integrated modeling module calculates the spatial gradient vector field representing the direction of fluid migration in order to characterize the physical migration trajectory of geothermal fluids in space. Specifically, it performs the following steps.

[0091] S305, the integrated modeling module constructs a continuous spatial grid model of the indices. Since the collected geochemical data are spatially discrete, the system needs to convert the discrete integrated hydrodynamic zonation indices into continuous two-dimensional spatial functions. The integrated modeling module uses the ordinary kriging interpolation algorithm to perform gridding processing on the study area.

[0092] For the semivariogram fitting and spatial weight calculation in ordinary Kriging interpolation, those skilled in the art can use conventional spatial statistical methods to achieve it. The spatial interpolation process is a well-known technique in this field and will not be described in detail here.

[0093] During the meshing process, the integrated modeling module sets the interpolation grid spacing parameter based on the average sampling point spacing of the study area. This interpolation grid spacing parameter controls the resolution of the generated grid, and its value is set between 20% and 50% of the average sampling point spacing. The integrated modeling module reads the survey area parameter from the system configuration file, for example, setting this proportion to 30%. After interpolation calculation, the discrete integrated hydrodynamic zonation indexes are mapped to a continuous two-dimensional coordinate plane, resulting in grid nodes. Interpolated values ​​of the comprehensive hydrodynamic zonation index at the location .

[0094] S306, the integrated modeling module calculates the spatial gradient components of each node on a continuous grid surface. Geological fluid migration is driven by temperature and pressure gradients, which manifest in geochemical characteristics as the spatial variation rate of the comprehensive fluid dynamics zonation index. The integrated modeling module uses the central difference algorithm to solve for the partial derivatives along the horizontal and vertical projected coordinates. The calculation formula is as follows:

[0095] ;

[0096] ;

[0097] in, This represents the gradient component of the grid node along the horizontal axis. This represents the gradient component of the grid node along the vertical axis. Represents grid nodes Interpolated values ​​of the comprehensive fluid dynamics zoning index at the location; , , and These represent the interpolated values ​​of the comprehensive fluid dynamics zoning index at the corresponding offset coordinate points; and These represent the interpolation grid spacing parameters set for the horizontal and vertical axes, respectively. Through the above calculations, orthogonal directional characteristic quantities reflecting the rate of change of the index were obtained.

[0098] S307, the integrated modeling module synthesizes spatial gradient vectors and maps them to the direction of fluid migration. The integrated modeling module uses calculated orthogonal gradient components to synthesize complete gradient vectors at each grid node. During the migration of geothermal fluids from deep to shallow depths, the interpolated values ​​of their integrated hydrodynamic zonation indices exhibit a decreasing trend from high to low. Based on this trend, the physical migration direction of the fluid is consistent with the direction of the negative gradient where the interpolated values ​​decay.

[0099] The integrated modeling module extracts the negative gradient direction as the fluid transport direction vector and calculates the magnitude of this vector. The magnitude is calculated as the square root of the sum of the squares of the gradient components along the horizontal and vertical axes, i.e. .in, This represents the magnitude of the spatial gradient vector. This magnitude value is used to quantitatively assess the drastic degree of change in indicators within a local space. The larger the magnitude value, the more drastic the change in the fluid's physicochemical environment, corresponding to the development of heat conduction channels or tectonic fracture zones for fluid migration.

[0100] The integrated modeling module assigns vector data with direction and magnitude attributes to corresponding grid nodes, generating a spatial gradient vector field covering the entire study area. Mathematically, this vector field transforms the spatial variation of geochemical elements into directional vectors, providing a geometric basis for the system to track the location of deep heat sources and determine fluid upwelling channels.

[0101] After obtaining the comprehensive geothermal anomaly index at each discrete coordinate point, the spatial analysis and visualization module performs the following steps to construct a continuous scalar field and perform anomaly threshold segmentation in order to achieve continuous characterization of the spatial range and quantitative extraction of the target area boundary.

[0102] S401, the Spatial Analysis and Visualization module constructs a continuous scalar field of comprehensive geothermal anomaly indices. This module extracts the comprehensive geothermal anomaly indices at each discrete sampling point. To extend the attribute values ​​of these discrete points to the entire study area and to conform to the physical transport patterns of geothermal fluids, the module uses anisotropic kriging interpolation from geostatistics to construct a continuous scalar field of anomaly indices.

[0103] Because the migration of deep geological fluids is strictly controlled by the regional main fault system, the resistance to fluid migration along the fault strike is much smaller than the resistance perpendicular to the fault strike. The spatial analysis and visualization module imports known vector data of the main structural lines in the study area from external sources and calculates the strike azimuth parameter of these main structural lines. When performing anisotropic interpolation calculations, the spatial analysis and visualization module sets this strike azimuth parameter as the guiding parameter for spatial interpolation.

[0104] When fitting the spatial semivariogram model, the system forcibly aligns the major axis of the search ellipse to the azimuth parameter of the structural line and assigns a spatial autocorrelation weight coefficient greater for the major axis than for the minor axis. Through these anisotropic parameter settings, the generated comprehensive geothermal anomaly index continuous surface accurately reflects the geometric features distributed along the structural zone, eliminating the spherical morphology distortion caused by conventional isotropic interpolation. (Mesh node acquisition) The interpolated value of the comprehensive geothermal anomaly index at the location is denoted as... .

[0105] S402, the spatial analysis and visualization module calculates multi-level anomaly thresholds based on statistical distributions. The values ​​in the continuous scalar field reflect the relative levels of geological indicators, and the system needs to define the boundaries of different levels of anomalies using mathematical methods.

[0106] The spatial analysis and visualization module extracts the comprehensive geothermal anomaly index corresponding to all valid sampling points within the study area and calculates its arithmetic mean. with standard deviation Based on the aforementioned statistical distribution parameters, the spatial analysis and visualization module adaptively defines the lower thresholds for weak anomalies, obvious anomalies, and core anomalies. The calculation formula is as follows:

[0107] ;

[0108] ;

[0109] ;

[0110] in, , and These represent the lower threshold values ​​for dividing the weak anomaly region, the obvious anomaly region, and the core anomaly region, respectively. , and These represent the anomaly classification coefficients for the corresponding levels.

[0111] This anomaly grading coefficient controls the stringency of anomaly delineation. The spatial analysis and visualization module reads the specific values ​​of this coefficient from the system configuration file. The value range is set to be between 1.0 and 1.5, for example, a value of 1.5; The value range is set to be between 2.0 and 2.5, for example, a value of 2.0; The value range is set between 3.0 and 4.0, for example, a value of 3.0. Through this hierarchical calculation logic, the system establishes a quantitative discrimination criterion for peeling away geothermal anomaly features layer by layer from a continuous background.

[0112] S403, the Spatial Analysis and Visualization module performs geostatistical thresholding and boundary delineation of the continuous scalar field. The module utilizes calculated multi-level anomaly thresholds to perform pixel-level segmentation of the generated continuous scalar field.

[0113] The spatial analysis and visualization module traverses each grid node in the quantization field and interpolates the comprehensive geothermal anomaly index at that node. Compare the data with the various abnormal thresholds within a given range. lie in and When the node is in between, mark it as a weakly anomalous node; when lie in and When the interval is between, mark it as a clearly abnormal node; when Greater than or equal to When, mark as a core abnormal node; when Less than When that happens, mark it as a background node.

[0114] After completing the attribute labeling of all nodes, the spatial analysis and visualization module extracts the outer envelopes of adjacent nodes at the same level. By drawing closed contour polygons on a two-dimensional spatial plane, nested boundaries of geothermal anomaly zones are generated. This layered delineation operation effectively filters out wide-area geological background noise and intuitively depicts the spatial progression of geothermal anomalies from the periphery to the center.

[0115] After delineating multi-level anomaly target areas, the spatial analysis and visualization module performs spatial registration and feature analysis of scalar and vector fields to reveal the dynamic relationship between geothermal anomalies and fluid migration. Specifically, the following steps are executed.

[0116] S404, the Spatial Analysis and Visualization module performs isomorphic mapping and overlay of scalar and vector fields. The module uses the geothermal anomaly polygons at various levels, generated through threshold segmentation, as the base rendering layer to present the spatial distribution of anomaly intensity.

[0117] Above the base layer, the spatial analysis and visualization module imports the spatial gradient vector field data generated in the preceding calculations. For each spatial gradient vector at a grid node, the system converts it into a directional geometric line segment. The direction of the geometric line segment is the negative gradient direction of the vector, used to indicate the fluid transport trajectory; the length of the geometric line segment is positively correlated with the magnitude of the vector, used to characterize the relative dynamic intensity of the fluid transport.

[0118] For the specific graphics operations of vector symbolic rendering and multi-layer spatial registration, those skilled in the art can use conventional geographic information system technology to achieve them. The layer overlay and geometric rendering processes are well-known technologies in this field and will not be described in detail here.

[0119] S405, the Spatial Analysis and Visualization module calculates the divergence of the vector field to quantitatively identify fluid divergence characteristics. In a typical hydrothermal upwelling pattern, as fluid migrates from the deep main channel to the shallower region, it exhibits a divergent vector topology radiating outwards from the center on a two-dimensional plane. The Spatial Analysis and Visualization module extracts this geometric pattern by calculating the divergence of the vector field.

[0120] The mathematical basis for divergence calculation is to solve for the sum of the partial derivatives of the vector components along each coordinate axis. Based on the constructed discrete grid data, the spatial analysis and visualization module uses the central difference algorithm to calculate the divergence value at each grid node. The calculation formula is as follows:

[0121] ;

[0122] in, Represents grid nodes The divergence value at; and These represent the gradient components along the horizontal axis at the corresponding offset coordinate points; and These represent the gradient components along the vertical axis at the corresponding offset coordinate points; and These represent the interpolation grid spacing parameters in the horizontal and vertical directions, respectively.

[0123] When the divergence value A positive value indicates a net outflow of fluid vectors at that grid node, exhibiting spatial divergence characteristics. The larger the value, the more pronounced the outward radiation trend from that node.

[0124] S406, the spatial analysis and visualization module combines divergence thresholding and scalar anomalies to identify hydrothermal upwelling dynamic centers. The module extracts all divergence values ​​greater than zero within the study area and calculates the arithmetic mean and standard deviation of this positive subset.

[0125] The spatial analysis and visualization module constructs a divergence identification threshold based on this, calculated using the following formula: .in, Indicates the divergence recognition threshold; This represents the arithmetic mean of all positive divergence values; This represents the standard deviation of the above positive divergence values; This represents the divergence intensity coefficient, with a value range of 1.5 to 2.0. The spatial analysis and visualization module reads this coefficient value from the system configuration, for example, setting it to 1.5.

[0126] The spatial analysis and visualization module traverses all grid nodes and extracts divergence values ​​that are greater than or equal to the divergence identification threshold. The node is marked as a strongly divergent node.

[0127] Subsequently, the spatial analysis and visualization module performs a joint spatial topology comparison of scalars and vectors. The system extracts a set of grid nodes that simultaneously possess both core anomaly node attributes and strong divergent node attributes. These grid nodes not only meet the high value requirements of the comprehensive geothermal anomaly index in space, but also exhibit strong outward fluid radiation geometry, forming a specific topological pattern of superimposed high anomaly centers and radial vector sources.

[0128] The spatial analysis and visualization module merges these extracted overlapping nodes into a continuous spatial polygon. This polygonal region is systematically marked as the main channel for deep hydrothermal upwelling or the vertical projection area of ​​a deep heat source, providing a direct dynamic spatial positioning benchmark for final geological interpretation and borehole layout.

[0129] After receiving the multi-source analysis results, such as the hydrothermal upwelling dynamic center extracted by the spatial analysis and visualization module, the results output and target area delineation module introduces multi-source prior geological constraint information for spatial registration and overlay, providing a physical verification benchmark independent of geochemical data for the final selection of mineral exploration target areas. The specific steps are as follows.

[0130] S501, the Results Output and Target Area Delineation module imports multi-source prior geological constraint data. A single geochemical anomaly needs to be validated in conjunction with the regional geological background to exclude interference from element enrichment caused by non-geothermal factors. The Results Output and Target Area Delineation module reads known deep fault zone distribution maps, magmatic rock outcrop maps, and remote sensing surface temperature anomaly data within the study area.

[0131] In geological mineralization and thermal mechanisms, deep faults typically provide the main deep circulation and upwelling channels for geothermal fluids. Igneous rock distribution areas often correspond to potential deep local heat source centers, while remote sensing surface temperature anomaly data reflects shallow thermal radiation and weak thermal seepage. The results output and target area delineation module loads multiple source files of different data formats into memory to establish an initial geological constraint dataset.

[0132] S502, the output results and target area delineation module perform unified spatial registration and grid alignment. Because the above-mentioned multi-source geological constraint information has different spatial reference systems, data modes and scales in the original state, the system needs to eliminate spatial reference differences to achieve tight overlay between layers.

[0133] The results output and target delineation module performs geometric coordinate reprojection on vector polygon data such as deep fault surfaces and magmatic rock outcrops, mapping them to a two-dimensional planar projection coordinate system consistent with the preceding geochemical anomaly scalar field. For raster image data such as remotely sensed surface temperature anomalies, the results output and target delineation module maintains the same interpolation grid resolution as in the preceding grid generation process, performs resampling operations, and forces alignment of pixel boundaries.

[0134] For affine transformations of vector coordinates and bilinear interpolation resampling algorithms for raster pixels, those skilled in the art can use conventional geographic information system (GIS) underlying libraries to implement them. The spatial registration calculation process is a well-known technology in this field and will not be elaborated here.

[0135] S503, the Results Output and Target Area Delineation module constructs standardized geological constraint surfaces based on the registered elements. In order to input the registered geological elements into subsequent logical operations, the Results Output and Target Area Delineation module needs to convert them into standardized two-dimensional polygons that characterize the geothermal activity potential.

[0136] For deep fault data, the results output and target area delineation module extracts the central axis of the fault and sets a buffer distance parameter to generate a buffer zone of the fault influence zone. This buffer distance parameter is used to define the geological influence range of rock fracturing and increased permeability around the fault zone. Its value range is set between 500 meters and 2000 meters. The results output and target area delineation module reads the corresponding distance value according to the structural level of the fault and generates a smooth buffer polygon as a structural constraint surface.

[0137] For data on the outcropping range of magmatic rocks, the output and target area delineation modules extract the outer envelope of the reprojected polygons to generate a heat source constraint surface. For the registered remote sensing surface temperature raster data, the output and target area delineation modules extract a set of continuous pixels with temperature values ​​higher than a set background temperature threshold. This background temperature threshold is determined based on the average surface temperature of the study area during the same historical period. The system converts the extracted set of continuous pixels into a surface thermal anomaly surface with closed boundaries.

[0138] Through the above geometric topology operations, the output and target area delineation modules uniformly convert geological prior information from different sources into a geometric polygon layer with a definite spatial location, thus completing the spatial preprocessing of independent verification data.

[0139] After obtaining standardized multi-source geological constraint surfaces and geochemical analysis results, the results output and target area delineation module performs the final mineral exploration target area delineation through spatial Boolean logic operations and outputs comprehensive interpretation results, specifically executing the following steps.

[0140] S504, the results output and target area delineation module performs spatial Boolean logic intersection operations to extract the target area range. Single geochemical anomalies are often false anomalies due to local lithological differences or subsequent weathering. The results output and target area delineation module uses the extracted geological constraint surface as a priori filter condition and performs a spatial logic comparison with the previously generated hydrothermal upwelling dynamic center polygon.

[0141] The results output and target delineation module uses the spatial Boolean intersection algorithm for calculation. The formula for its spatial set operation is as follows:

[0142] ;

[0143] in, This represents the set of spatial polygons that ultimately delineate the geothermal target area; This represents the set of spatial polygons corresponding to the hydrothermal upwelling dynamic center. A set of spatial polygons representing the constrained surfaces; A set of spatial polygons representing the constraint surfaces of a heat source; A set of spatial polygons representing surface thermal anomalies; The operator represents the intersection of spaces; This represents the operator for finding the union of spaces.

[0144] Using this logical formula, the output results and the target area delineation module determine that the final target area must be located within a deep geochemical and dynamic anomaly zone, and spatially overlap with at least one of the following: a tectonic fault influence zone, a magmatic heat source zone, or a shallow thermal anomaly zone. This operation eliminates isolated geochemical false anomalies lacking geological support, improving the accuracy of geothermal exploration target area location.

[0145] S505, the results output and target area delineation module performs target area geometric attribute extraction and boundary smoothing. The initial intersection polygons generated by multi-layer Boolean logic operations typically contain irregular jagged boundaries or extremely small fragmented patches. The results output and target area delineation module calculates the actual geographic area of ​​each polygon and sets a patch area threshold.

[0146] The area threshold for this map patch is determined based on the minimum working face requirements of actual geothermal exploration projects, and its value range is set at 0.2 km. 2 Up to 1.0km 2 Between, for example, reading the set value as 0.5km 2 The output and target area delineation module will remove small polygons with an area smaller than the area threshold of the patch.

[0147] Subsequently, the results output and target area delineation module performs curve smoothing on the remaining polygon boundaries to eliminate jagged nodes and generate a final geothermal target area polygon with smoothed boundaries. For polygon area calculation and Bézier curve smoothing of polygon boundaries, those skilled in the art can use conventional computational geometry algorithms; the graphic smoothing process is a well-known technique in the field and will not be elaborated upon here.

[0148] After completing the graphic processing, the output and target area delineation module extracts and records the geometric centroid coordinates, boundary node coordinate sequence, and corresponding area values ​​of the final geothermal target area polygon.

[0149] S506, the Results Output and Target Area Delineation module generates a comprehensive geological interpretation and target area analysis report. The system needs to transform the calculated data from various dimensions into a guideline document readable by engineers. The Results Output and Target Area Delineation module uses the final geothermal target area polygon as the spatial boundary and extracts key attribute information within that area from the underlying database.

[0150] The extracted attribute information specifically includes the peak value and arithmetic mean of the comprehensive geothermal anomaly index within the target area, the divergence azimuth of the spatial gradient vector, the types of geological constraint elements superimposed within the target area, and the combination characteristics of key anomaly elements such as potassium, lithium, and fluorine within the target area.

[0151] Based on the above information, the results output and target area delineation module perform logical inference of the genetic type of the geothermal system. When the system reads geochemical parameters within the target area showing an enrichment combination of high potassium, high lithium, and high fluorine elements, and the spatial topological boundary of the target area is directly adjacent to or superimposed on the magmatic heat source constraint surface, the system infers it to be a magmatic geothermal system; when the target area lacks extremely high anomalies of the above volatile elements, and is strictly controlled spatially by tectonic fault constraint surfaces, the system infers it to be a convective geothermal system.

[0152] The results output and target area delineation module fills in the geometric coordinate details, index values, and the above-mentioned causal inference results into the corresponding data fields. Based on this information, the system automatically outputs a comprehensive geological interpretation report containing the exploration target area coordinates, anomaly intensity level assessment, heat source flow control channel location, and inference of geothermal system causal type, providing data support for field borehole deployment.

[0153] Specific application examples:

[0154] Taking the exploration of a typical concealed high-temperature geothermal field (the study area covers approximately 50 square kilometers) as an example, the steps S1 to S5 were systematically and completely executed:

[0155] 200 whole-rock geochemical sampling points were arranged in a 500m × 500m grid across the study area. When the data integration and preprocessing module read the raw data, it was found that the antimony (Sb) concentration at sampling point P32 was displayed as a null value because it was below the instrument detection limit (0.1 ppm). Following the logic of step S102, the system extracted half of the detection limit (0.05 ppm) and assigned it a value; subsequently, it applied the formula... Dimensionless processing was performed on all element data.

[0156] When calculating the thermal corrosion intensity index, taking sampling point P55 around the core anomaly zone as an example, the system extracts its original concentrations of major elements: potassium oxide 4.5%, magnesium oxide 2.0%, sodium oxide 1.5%, and calcium oxide 1.0%. Substituting these into the calculation formula in step S201:

[0157] ;

[0158] Simultaneously, the system extracts rubidium (Rb, concentration 150 ppm) as an indicator element of the deep high-temperature phase and arsenic (As, concentration 10 ppm) as an indicator element of the shallow low-temperature phase, and substitutes them into step S206 to calculate the concentration ratio of the characteristic element pair. The value is 15. Along the fluid migration path from the periphery to the center (from 4 km to 0 km from the kinetic center), the thermal alteration intensity index gradually increases from 35% to 75%. The ratio increases from 2 to 40, showing a clear spatial evolution gradient.

[0159] The integrated modeling module utilizes principal component analysis to analyze the thermal corrosion intensity index and alkali migration index (…). Dimensionality reduction was performed using parameters such as the Rb / As ratio. The calculation results showed that the variance contribution rate of the first principal component reached 68.5% (greater than the 40% principal component validity threshold). The system extracted the score of each point from the first principal component as a comprehensive geothermal anomaly index, and used the central difference algorithm in step S306 to obtain the horizontal and vertical partial derivatives of each grid node, synthesizing a spatial gradient vector field pointing towards the dynamic center.

[0160] The system calculates the arithmetic mean of the comprehensive geothermal anomaly index. Standard deviation Set the anomaly classification coefficient for the core anomaly region. Calculate the lower threshold of the core anomaly region. The system extracted core anomaly nodes with a comprehensive geothermal anomaly index ≥ 4.5, and combined them with the divergence identification formula ( Strongly divergent nodes were extracted, establishing the hydrothermal upwelling dynamic center. Finally, Boolean logic formulas were used:

[0161] ;

[0162] By performing spatial intersection calculations between the geochemical anomaly center and the prior deep fault buffer zone (buffer distance 800m), a target area of ​​0.85 square kilometers (greater than the set threshold of 0.5) was successfully delineated, and the causal inference was output as a convective geothermal system.

[0163] To verify the accuracy and engineering practical value of the technical solution of this invention, the prospecting target area delineation method based on scalar vector joint field and multidimensional geological constraints adopted in this invention was compared with the traditional single element concentration anomaly lower limit delineation method (which only uses areas where the arsenic concentration is greater than the mean plus 2 times the standard deviation as the target area) through drilling engineering verification and effect comparison in the same study area.

[0164] Traditional single-element methods, due to their failure to consider the constraints of main faults and the complex evolution characteristics of multiple elements, are affected by local lithological fluctuations, resulting in anomaly areas of up to 3.2 square kilometers, and presenting multiple isolated and scattered patches. In contrast, this invention relies on spatial Boolean intersection operations (intersection of the hydrothermal upwelling dynamic center and the buffer zone of deep and large faults) to precisely lock the core target area to 0.85 square kilometers, successfully filtering out 73.4% of weak peripheral anomalies and false surface anomalies without solid tectonic support, thus narrowing the exploration range.

[0165] The traditional single-element method recommended three verification boreholes. Ultimately, only one borehole showed heat. Although the other two boreholes had trace element enrichment at the surface, they completely lacked fluid upwelling channels at depth, which is a typical example of shallow lateral seepage.

[0166] Two verification boreholes were recommended using the method of this invention. Both boreholes accurately encountered the main thermal conductivity fracture zone at a depth of 1200 to 1500 meters, with bottom temperatures reaching 165°C and 172°C respectively, demonstrating the extremely high accuracy of this invention in locating hidden geothermal target areas.

[0167] See attached document Figure 3 The horizontal axis represents the distance (in kilometers) from the center of the geothermal anomaly dynamics. The scale in the figure gradually decreases from 4 to 0 from left to right, simulating the spatial path of the geothermal system's outer cold water zone continuously approaching the deep heat source in the system's core.

[0168] The vertical axis on the left and the solid broken line represent the hydrothermal alteration intensity index. This index corresponds to the ratio of major elements (potassium oxide, magnesium oxide, sodium oxide, and calcium oxide) in the specific implementation. As shown by the solid line marked with a square in the figure, the index is only about 35 at a distance of 4 kilometers from the center. As the space moves closer to the dynamic center (0 kilometers), the fluid temperature and the intensity of the water-rock reaction increase dramatically. Feldspar minerals are extensively destroyed and high-temperature alteration minerals are generated, causing the index to rise synchronously in an almost linear fashion, eventually reaching a peak value of 75 at the center.

[0169] The vertical axis on the right and the dashed line represent the concentration ratio of the characteristic element. This index corresponds to the ratio of the deep high-temperature phase indicator element (rubidium) to the shallow low-temperature phase indicator element (arsenic) in the specific implementation. As shown by the dashed line with circular markings in the figure, at a distance of 4 kilometers from the center, highly volatile and easily migrating arsenic is heavily enriched, resulting in a low ratio of less than 5. However, as the distance approaches the kinetic center (0 kilometers), the concentration of rubidium, representing the high-temperature phase, becomes absolutely dominant, and the ratio exhibits an exponential increase, eventually exceeding 40 at the center.

[0170] The trend of the two broken lines (the ratio of hydrothermal alteration intensity index to characteristic element concentration) in the attached figure increasing synchronously as the spatial distance decreases confirms the theoretical model that the direction of physical transport of the fluid is consistent with the direction of the negative gradient of the interpolated numerical decay. Figure 3 The horizontal axis represents the peak values ​​of various indicators at 0 km, intuitively reflecting the increasing distribution pattern of geochemical indicators from the periphery to the center of the heat source. The calculation model of this invention synthesizes a two-dimensional vector field and calculates the divergence based on this spatial gradient law, thereby achieving quantitative positioning of the hydrothermal upwelling dynamic center. This calculation process provides quantified spatial coordinate constraints for borehole layout, effectively improving the borehole heat encounter rate in engineering verification.

Claims

1. A method for identifying geothermal anomaly zones based on whole-rock geochemical data, characterized in that, include: The raw whole-rock geochemical dataset was obtained, and data cleaning and spatial mapping were performed to form a standardized discrete spatial dataset. Receive the discrete spatial dataset, calculate the geothermal activity indicator scalar and characteristic element concentration ratio variables, synthesize the comprehensive fluid dynamics zonation index and generate a structured data table; Receive the structured data table, calculate the comprehensive geothermal anomaly index through principal component analysis, and convert the comprehensive hydrodynamic zonation index into a spatial gradient vector field; Receive the comprehensive geothermal anomaly index and the spatial gradient vector field, construct a continuous scalar field and define multi-level anomaly boundaries, and extract the hydrothermal upwelling dynamic center by combining the divergence value of the spatial gradient vector field. The system receives the hydrothermal upwelling dynamic center, performs spatial intersection judgment based on the imported multi-source prior geological constraint data, delineates the target area boundary set, and outputs comprehensive geological interpretation results.

2. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 1, characterized in that, The steps of obtaining the raw whole-rock geochemical dataset, performing data cleaning and spatial mapping to form a standardized discrete spatial dataset specifically include: The instrument detection limit parameter is used to compare the elemental concentration values ​​in the original whole-rock geochemical dataset with the values ​​of the instrument detection limit parameter, and the outlier removal mechanism is activated to clear invalid values ​​and complete the data cleaning. For missing records, spatial imputation of missing values ​​is performed, and the standard fraction algorithm is used to complete the dimensionless processing; The latitude and longitude information is obtained and converted into planar projected coordinates. The planar projected coordinates are then associated and bound with the standardized concentration vector to construct the discrete spatial dataset.

3. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 1, characterized in that, The steps for calculating the scalar indicator of geothermal activity and the ratio of characteristic element concentrations specifically include: The original concentration values ​​of major elements were extracted to calculate the hydrothermal alteration intensity index, and the original concentration values ​​of potassium oxide and sodium oxide were extracted and the quotients were calculated to generate the alkali migration index. The enrichment coefficient is calculated by extracting the original concentration values ​​of the target indicator element and comparing them with the regional background concentration. The original concentration values ​​of the indicator elements in the deep high-temperature phase and the indicator elements in the shallow low-temperature phase are extracted and the quotient is calculated to obtain the concentration ratio variable of the characteristic elements.

4. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 1, characterized in that, The steps of synthesizing comprehensive fluid dynamics zonation indices and generating structured data tables specifically include: The weight coefficients of feature element pairs are obtained through principal component analysis algorithm; The comprehensive fluid dynamics zoning index is calculated by performing a weighted average on the concentration ratio variables of multiple characteristic elements with similar physical meanings. The generated geothermal activity scalar and the integrated fluid dynamics zonation index are packaged and converted to generate the structured data table.

5. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 1, characterized in that, The steps for calculating the comprehensive geothermal anomaly index through principal component analysis specifically include: The hydrothermal alteration intensity index, alkali migration index and enrichment coefficient are extracted from the structured data table to construct an initial index matrix, and standardized index matrix is ​​obtained by performing standardization processing on columns. Calculate the correlation coefficient matrix of the standardized index matrix and perform eigenvalue decomposition to obtain the corresponding eigenvalues ​​and eigenvectors; The comprehensive geothermal anomaly index is defined by multiplying the eigenvector corresponding to the first principal component with the standardized index matrix and extracting the score value of each sampling point on the first principal component.

6. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 5, characterized in that, The specific steps for converting the comprehensive fluid dynamics zonal index into a spatial gradient vector field include: The discrete integrated fluid dynamics zoning indexes are meshed using the ordinary Kriging interpolation algorithm to obtain interpolated values; The gradient components are obtained by solving the partial derivatives in the directions of the horizontal and vertical projected coordinates using the central difference algorithm. The complete gradient vectors at each grid node are synthesized, the negative gradient direction is extracted as the fluid transport direction vector and the magnitude is calculated to generate the spatial gradient vector field.

7. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 1, characterized in that, The steps of constructing a continuous scalar field and defining multi-level anomaly boundaries specifically include: The anisotropic kriging interpolation algorithm is invoked, and the strike azimuth parameter of the main structural line is set as the steering parameter to construct the continuous scalar field of the comprehensive geothermal anomaly index. Extract the comprehensive geothermal anomaly index corresponding to all valid sampling points and calculate the arithmetic mean and standard deviation, and define the lower limit thresholds for weak anomalies, obvious anomalies, and core anomalies; The continuous scalar field is segmented at the pixel level using a multi-level anomaly threshold, and the outer envelopes of adjacent nodes at the same level are extracted to generate the multi-level anomaly boundaries.

8. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 1, characterized in that, The specific steps for extracting the hydrothermal upwelling dynamic center by combining the divergence value of the spatial gradient vector field include: The divergence value at each grid node is calculated using the central difference algorithm. Extract divergence values ​​greater than zero and calculate the arithmetic mean and standard deviation of the positive subset to construct a divergence identification threshold. Nodes with divergence values ​​not less than the divergence identification threshold are marked as strongly divergent nodes. Perform a joint spatial topology comparison of scalar and vector, extract the set of mesh nodes that simultaneously possess core anomalous node attributes and strong divergent node attributes, and merge them into a continuous spatial polygon, which is marked as the hydrothermal upwelling dynamic center.

9. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 1, characterized in that, The steps for performing spatial intersection determination using imported multi-source prior geological constraint data specifically include: Perform unified spatial registration and grid alignment on the distribution map of deep and large fault zones, the map of magmatic rock outcrops, and remote sensing surface temperature anomaly data; Extract the fracture center axis to generate a buffer polygon as a structural constraint surface, extract the outer envelope of igneous rocks to generate a heat source constraint surface, and extract a continuous set of pixels with temperature values ​​higher than the background temperature threshold to convert them into a surface thermal anomaly surface. The spatial Boolean intersection algorithm is invoked to determine that the target area is located within the deep geochemical dynamic anomaly zone and overlaps spatially with the tectonic constraint surface, the heat source constraint surface, or the surface thermal anomaly surface, thereby generating an initial intersection polygon.

10. The method for determining geothermal anomaly zones based on whole-rock geochemical data according to claim 9, characterized in that, The steps of delineating the target area boundary set and outputting the comprehensive geological interpretation results specifically include: The initial intersection polygons are subjected to area threshold culling and curve smoothing to generate the final geothermal target area polygons, which constitute the target target area boundary set. Extract the attribute information within the final geothermal target area polygon to perform logical inference of the geothermal system's genetic type; The output includes a comprehensive geological interpretation report containing the polygon coordinate parameters of the final geothermal target area and the inference of its genetic type.

Images