Geothermal water multi-factor remote sensing detection method, system, equipment, medium and product

Abstract

The invention discloses a geothermal water multi-factor remote sensing detection method, system and device, a medium and a product, and the method comprises the steps: obtaining a thermal radiation anomaly index and a surface temperature anomaly index through multi-temporal thermal infrared preliminary screening and surface temperature inversion; the method comprises the following steps: extracting a vegetation stress index VIBS by utilizing the characteristic that terrestrial heat inhibits vegetation, and obtaining a comprehensive altered mineral index MI according to mineral alteration characteristics and a spectral unmixing method; based on ALOS PALSAR, a fracture-fracture zone is extracted, and a construction advantage degree SF is quantified; gRACE gravity image data are introduced to form hydrological favorable degree HF; geothermal water indication characteristics are integrated, multi-factor scoring and grading are fused to obtain a comprehensive geothermal anomaly distribution map, and finally boundary refinement and field verification are carried out through unmanned aerial vehicle thermal infrared / multispectral aerial survey. According to the method, false detection and missing detection can be remarkably reduced, and the accuracy and field efficiency of geothermal anomaly recognition are improved.

Description

Technical Field

[0001] This invention belongs to the field of geology and remote sensing technology, specifically relating to a multi-factor remote sensing method, system, equipment, medium and product for geothermal water. Background Technology

[0002] Geothermal systems on the surface often exhibit thermal anomalies, vegetation stress, acidic sulfate-type altered minerals, structural features such as fault-fracture zones acting as water and heat conduction channels, and hydrological anomalies related to groundwater recharge / flow. Traditional single data sources (such as single-scene thermal infrared or single-mineral interpretation) are easily affected by factors such as topography, season, water thermal inertia, urban heat islands, and bare rock albedo, leading to misjudgments and omissions. In recent years, multi-temporal TIR combined with temperature-emissivity separation (TES) at night has improved the accuracy of surface temperature retrieval; hyperspectral imaging can identify geothermal vegetation stress characteristics characterized by red-edge blue shift and SWIR reflectance uplift in vegetated areas, and quantitatively retrieve alteration minerals such as kaolinite and alunite in bare / sparsely vegetated areas through linear spectral unmixing (LSU); radar ALOS-1 / 2 PALSAR can provide characterization of surface tectonic fracture zones and fault density; GRACE / GRACE-FO provides monthly-scale terrestrial water storage (TWS) anomalies, reflecting changes in regional hydrological background. Integrating these multi-factor evidences within a unified spatiotemporal framework, supplemented by on-site confirmation by UAV thermal infrared sounding, can significantly reduce uncertainty in remote sensing interpretation and improve field efficiency and accuracy. Summary of the Invention

[0003] To address the technical problem of misjudgment / missed judgment in geothermal remote sensing identification in complex coverage areas caused by the influence of topography, water bodies, season, and surface type, this invention aims to provide a multi-factor remote sensing detection method, system, equipment, medium, and product for geothermal water. It constructs a remote sensing detection technology with multi-factor constraints of "geothermal-vegetation / mineral-structure-hydrology" and UAV closed-loop verification, and realizes the extraction of geothermal anomaly range, priority, and boundaries.

[0004] To solve the technical problem, the technical solution of the present invention is as follows: A multi-factor remote sensing method for geothermal water, the multi-factor remote sensing method for geothermal water includes: ASTER and Landsat TIR images were acquired, and the surface temperature was obtained through radiometric calibration, atmospheric correction and TES inversion. Combined with cloud, water body and urban heat island removal, a thermal stability anomaly map H and a mask Mask_TIR were generated. Within the masked area of ​​GF-5 hyperspectral data, the image is analyzed using the mask to divide the region into vegetated and non-vegetated areas. For vegetated areas, the red edge and SWIR indices are calculated and standardized to synthesize the vegetation stress index VIBS. For non-vegetated areas, the abundance of kaolinite and alunite is extracted using the endmember unmixing method to construct the alteration index MI. Fracture structures were extracted using ALOS PALSAR and tectonic factors (SF) were generated. These were then combined with Grace TWS data to decompose water storage trends, amplitudes, and frequencies, thus constructing hydrological favorability (HF). TWS ; Finally, the H, VIBS, MI, SF, and HF components are fused. TWS The comprehensive score is calculated using multiple factors and normalized weighting. Candidate anomalous patches are obtained through clustering and classification. The results are then verified by UAV thermal infrared and multispectral analysis to output accurate geothermal anomaly distribution results.

[0005] Furthermore, the generation of the thermal stability anomaly map H and the mask Mask_TIR includes: At least 6 ASTER TIR images were acquired, supplemented with Landsat-8 / 9 TIRS. Radiometric calibration, atmospheric correction, and geometric registration were performed on each image. The LST (Land Surface Temperature) was retrieved using the TES algorithm. During the processing of each image, cloud masks and the JRC water body dataset were applied to remove interference from cloud and water reflections. In addition, significant urban heat island effects were identified and removed using land use data. Finally, a stable thermal anomaly map H and a temperature anomaly mask Mask_TIR (without water and urban interference) were obtained.

[0006] Furthermore, the construction of the alteration index MI specifically includes: Regional division: Acquire GF-5 AHSI hyperspectral imagery, perform atmospheric correction on the hyperspectral imagery, mosaic multiple images together, and calculate the Normalized Difference Vegetation Index (NDVI) using the near-infrared and red bands. The formula is as follows:

[0007] Wherein, NIR is the reflectance in the near-infrared band, and R is the reflectance in the red band; Based on the NDVI value, areas with an NDVI greater than 0.30 are defined as vegetated areas, and other areas are defined as non-vegetated areas. The output generates Mask_TIR ∩ Vegetated Area and Mask_TIR ∩ Non-Vegetated Area zoning masks. Vegetation area treatment: The vegetation mask obtained above is used to calculate the mean of adjacent bands using the shortwave infrared (SWIR) bands in the hyperspectral image. NDVI1 is then calculated according to the definition, using the formula:

[0008] Among them, R2274 and R671 are the average reflectances of the SWIR long-wavelength band and red band, respectively; The slope of the red-edge band can be approximated by using the first derivative to calculate the reflectivity of the band near the red edge.

[0009] Where ΔR is the slope per unit wavelength; Then the red-edge blue shift index was calculated. To monitor the impact of geothermal water on vegetation;

[0010] Standardize NDVI1 and NDVI2, calculate their mean and standard deviation in the vegetation zone, and apply the Z-score standardization formula:

[0011] in, and These are the mean and standard deviation, respectively. By combining the standardized NDVI1 and NDVI2, the composite index VIBS is obtained. The stress level of vegetation is assessed by a formula, and the vegetation stress index map VIBS is output, covering the Mask_TIR∩ vegetation zone.

[0012] Treatment of non-vegetated areas: Based on the acquired non-vegetated area mask, endmember unmixing was performed using hyperspectral images. Kaolinite and alunite were selected as endmembers. The laboratory spectral data were first resampled to the GF-5 gas spectrum band, and the above spectrum was used for unmixing. The abundance of kaolinite and alunite in each pixel was extracted. After obtaining the abundance, the mean and standard deviation of the abundance map of kaolinite and alunite in the non-vegetated area were calculated. Z-score standardization was performed to ensure that the differences in illumination and background in different areas were eliminated. Using the standardized abundance values, a comprehensive alteration index MI is constructed, and a mineral anomaly index map MI is output, covering the Mask_TIR∩ non-vegetation area;

[0013] MI stands for Overall Alteration Index, which is dimensionless. The larger the value, the stronger the acidic sulfate-type alteration.

[0014] Furthermore, the extraction of fracture structures and the formation of the tectonic factor SF using ALOS PALSAR specifically includes: High-resolution radar data of the study area was acquired using ALOS PALSAR imagery. Image processing algorithms were used to perform edge detection on the radar images and identify geological structural features, including fault zones and linear structures. The detected structural line information was converted into a raster format for unified processing with other optical resolution datasets. Based on the spatial resolution of the imagery, the structural information was rasterized to the same resolution as the optical imagery. In the generated raster data, the number of fault lines per unit area and the number of intersections of different fault lines were calculated. These two indicators can be used to assess the complexity and activity of the regional structure. Finally, a structural favorability map (SF) was generated, which was used to quantify the potential impact of structural features on geothermal anomalies within the region. The construction of hydrological favorable HF TWS Specifically, it includes: Data from the GRACE satellite is used to monitor changes in Earth's water storage, specifically acquiring GRACE TWS data. Time-series analysis is performed on the GRACE TWS data, applying the STL method to decompose the water storage data into seasonal variations, long-term trends, and random fluctuations. The slope of the long-term trend (SLOPE) is calculated; a positive value indicates an increase in water storage, while a negative value indicates a decrease. The slope is expressed in cm / year to reflect the dynamic changes of the hydrological system. The seasonal variation amplitude (AMP) is calculated, representing the difference between the maximum and minimum water storage values ​​within a year; a larger amplitude indicates a more active hydrological cycle. The frequency of positive anomalies (PFREQ) within a year is counted to assess the persistence of hydrological conditions and their environmental impact. This value is normalized to between 0 and 1, representing the frequency proportion. Based on set weights, SLOPE, AMP, and PFREQ are weighted and combined to derive a comprehensive hydrological favorability index, which reflects the supply capacity and abundance of water resources. Finally, a hydrological favorability map (HF) is generated. TWS It is used to assess the potential impact of hydrological factors on geothermal anomalies.

[0015] Furthermore, obtaining accurate geothermal anomaly distribution results specifically includes: Perform multi-factor fusion scoring and priority allocation: Multi-factor fusion calculation is performed on all valid pixels within the Mask_TIR range, using the following formula:

[0016] Wherein, Score represents the overall advantage. Represents the weight of the i-th factor. This represents the normalized value of the i-th factor. The input factor data includes the thermal stability index H, vegetation stress index VIBS, mineral anomaly index MI, tectonic factor SF, and total water storage hydrological favorableness HF. TWS The result after 0-1 linear standardization; The fusion score is thresholded, connected component clustering analysis is performed, the mean and 90th percentile of each patch are calculated, and the patches are divided into high, medium and low priorities according to the quantile threshold, and a candidate anomaly list and level distribution map are output. UAV Detailed Survey and Confirmation: Based on the high-priority patch list output above, select the highest priority patch for detailed UAV survey, and then conduct UAV aerial survey, i.e., LWIR thermal imaging + multispectral / visible light, to generate centimeter-level temperature field, vegetation stress texture, and gas anomaly direction, and overlay H, VIBS, MI, SF, and HF to output high-precision geothermal anomaly confirmation results.

[0017] A multi-factor remote sensing system for geothermal water, the system being used to perform any of the methods described above, the system comprising: Data Access and Management Module: Used for unified collection and management of remote sensing data from multiple sources, including: ASTER and Landsat TIR images, GF-5 hyperspectral data, ALOS PALSAR images, DEM, JRC water body data, GRACE water storage data, and UAV data; Preprocessing and correction module: Used to perform a series of preprocessing on the acquired images to ensure data usability, including: radiometric calibration, atmospheric correction and geometric fine registration of ASTER and Landsat TIR images, inversion of land surface temperature using TES algorithm, removal of clouds, water bodies and urban heat island effect, generation of thermal stability anomaly map H and temperature anomaly mask Mask_TIR; Feature Extraction Module: This module extracts various features from the processed image, including: Thermal Anomalies: Generating thermal anomaly maps (H) and mask (Mask_TIR); Hyperspectral Vegetation Stress / Minerals: Analyzing vegetated and non-vegetated areas within the masked region of GF-5 hyperspectral data; Calculating and standardizing red edge and SWIR indices to synthesize the vegetation stress index (VIBS); Simultaneously, endmember unmixing is performed on non-vegetated areas to extract kaolinite and alunite abundances, constructing the alteration index (MI); Tectonic Landforms: Extracting fault structures using ALOS PALSAR to form tectonic factors (SF); Hydrological Field: Performing time-series analysis using GRACE TWS data to extract water storage trends, amplitudes, and frequencies, constructing the hydrological favorability index (HF). TWS ; Multi-factor fusion and hierarchical module: used to fuse and analyze the extracted features: fusion of thermal anomaly map H, vegetation stress index VIBS, mineral anomaly index MI, tectonic factor SF, and hydrological favorableness HF. TWS A comprehensive score is calculated using normalized weighted averages, and clustering and grading are performed to identify candidate anomalous patches. The UAV-based detailed investigation and ground verification module is used to further verify candidate anomalous patches: using UAVs equipped with thermal infrared and multispectral cameras to conduct detailed investigations of high-priority patches; collecting and analyzing UAV data to generate accurate geothermal anomaly distribution results, and comparing and confirming the authenticity of the anomalous areas; The Results Management and Publishing module is used for final processing and management of all obtained results, ensuring data reusability and convenient access for users, outputting the final geothermal anomaly distribution results, publishing them together with relevant reports or documents, and sharing key data and information with the outside world.

[0018] A computer device includes: a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor executes the computer program to implement the steps of the geothermal multi-factor remote sensing detection method described in any one of the above descriptions.

[0019] A computer-readable storage medium having a computer program stored thereon, characterized in that, when executed by a processor, the computer program implements the steps of the geothermal multi-factor remote sensing detection method described above.

[0020] A computer program product includes a computer program, characterized in that, when executed by a processor, the computer program implements the steps of the geothermal multi-factor remote sensing detection method described above.

[0021] Compared with the prior art, the advantages of the present invention are as follows: 1) Multi-factor verification significantly reduces misjudgments caused by bare rock, urban heat islands and thermal inertia of water bodies; 2) Sensitive to vegetation cover areas, VIBS captures vegetation stress caused by red-edge blue shift + SWIR elevation; 3) Quantitative characterization of acidic sulfate alteration minerals using MI-LSU-GF5 supports the determination of the origin of "why heat"; 4) SAR fracture density and structural geometry such as intersections / turns constrain hydrothermal channels; 5) GRACE TWS incorporates regional hydrological background to enhance the physical consistency of anomaly stability interpretation; 6) Closed-loop verification of drone thermal imaging and on-site inspection improves final inspection rate and field efficiency. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 Flowchart of the method; Figure 2 System functional module block diagram; Figure 3 A schematic diagram of multi-factor feature extraction; Figure 4 A schematic diagram of multi-factor fusion scoring and geothermal potential classification; Figure 5 1. Schematic diagram of on-site verification. Detailed Implementation

[0024] The specific implementation of the present invention is described below with reference to embodiments: It should be noted that the structures, proportions, sizes, etc. shown in this specification are only used to complement the content disclosed in the specification for those skilled in the art to understand and read, and are not intended to limit the conditions under which the present invention can be implemented. Any modifications to the structure, changes in the proportions, or adjustments to the size, without affecting the effects and objectives that the present invention can produce, should still fall within the scope of the technical content disclosed in the present invention.

[0025] Furthermore, the terms such as "upper," "lower," "left," "right," "middle," and "one" used in this specification are merely for clarity of description and are not intended to limit the scope of the invention. Any changes or adjustments to their relative relationships, without substantially altering the technical content, should also be considered within the scope of the invention. Example

[0026] like Figure 1 As shown in the figure, this embodiment provides a method for geothermal anomaly detection based on multi-source remote sensing data fusion, including the following steps: S1: ASTER thermal infrared surface temperature inversion and preliminary anomaly screening (delineating candidate range) At least six ASTER TIR images (B10–B14, 90 m, covering the dry / wet seasons) were acquired, supplemented with Landsat-8 / 9 TIRS images (100 m). Each image underwent radiometric calibration, atmospheric correction, and geometric registration. The LST (Land Surface Temperature) was retrieved using the TES algorithm. Cloud masks, JRC water bodies (Joint Research Centre Global Surface Water dataset), and land use data were used to remove significant urban heat islands. The final result was a stable thermal anomaly map (H) and a temperature anomaly mask (Mask_TIR) removed from water bodies and urban interference.

[0027] S: GF-5 hyperspectral differential processing S2.1 Region Delineation: The Normalized Difference Vegetation Index (NDVI) was calculated using the near-infrared and red bands of the mosaicked and atmospherically corrected GF-5 AHSI hyperspectral image.

[0028] Where NIR is the near-infrared reflectance and R is the red reflectance. An area is defined as vegetated when NDVI > 0.30; otherwise, it is defined as non-vegetated. This step is only performed within the "Thermal Anomaly Mask (Mask_TIR)" range, employing different processing strategies for vegetated and non-vegetated areas respectively.

[0029] S2.2 Vegetated Zone Treatment (VIBS): Within the “Mask_TIR ∩ Vegetation Zone”, two hyperspectral indices related to vegetation stress monitoring were calculated, as follows: 1) SWIR Long End Lift Index (NDVI1) First, take the average of adjacent bands (the horizontal line above indicates the "mean value of this band").

[0030]

[0031] Definition of index:

[0032] 2) Red-edge blue shift index (NDVI2) The first derivative approximation is made using the red-edge adjacent band (ΔR is the slope per unit wavelength).

[0033] in: , The reflectivity of adjacent bands (or windows) The wavelength interval is in μm or nm. The center wavelength is Reflectance (units can be consistent). :around Multi-band average reflectivity (noise reduction, robustness). by The first derivative approximation centered on the curve (spectral slope).

[0034] Selected for the purpose of geothermal water to suppress the red-edge blue shift of vegetation ΔR707, ΔR737 The calculation is as follows, where the denominator "10" is an example of the difference between the two points 707nm and 737nm, each taken by ±5 nm. This can be replaced according to the actual band interval:

[0035]

[0036] The index is defined as:

[0037] 3) Vibration Inhibition Strength of Composite Vegetation (VIBS) Calculate the scene statistics (mean and standard deviation) of the two indices within the “Mask_TIR ∩ Vegetation Zone”: First, standardize any index (Z-score general expression): given an index (or band ratio). I k (k can be 1, 2, 3..., such as NDVI, etc.), calculate the mean m within a specific mask area. k and standard deviation std k The standardized formula is:

[0038] For multiple standardized indices (or wave characteristics), a general composite index can be defined:

[0039] in: The weight of the k-th feature, For the standardized value, when =1 means the composite Z-score is of equal weight. Based on the suppression vegetation index characteristics of geothermal water according to the present invention, the determination is... m1, Mean and standard deviation of NDVI1; m2, Mean and standard deviation of NDVI2.

[0040] Standardize NDVI1 and NDVI2 (Z-score) and substitute them into the aforementioned formula to obtain the composite vegetation inhibition index:

[0041] Finally, output the vegetation stress intensity map VIBS, which only covers the “Mask_TIR ∩ Vegetation Zone”.

[0042] S2.3 Non-vegetated area treatment (MI-GF5): Within the “Mask_TIR ∩ Non-vegetated area”, linear spectral unmixing (LSU-GF5) was performed, with endmembers including kaolinite and alunite. The laboratory endmember spectra were resampled to the GF-5 band and then unmixed to obtain the endmember abundance map. The overall alteration index was calculated.

[0043] Where: MI: Overall Altering Index (MI-GF5), dimensionless. The larger the value, the stronger the acidic sulfate alteration (related to geothermal fluids).

[0044] Kaolinite LSU： The kaolinite endmember abundance obtained by linear spectral unmixing (LSU) represents the proportion of kaolinite in that pixel (0–1). The endmember spectra are resampled from laboratory / field spectral libraries to the actual GF-5 AHSI band before being used for unmixing.

[0045] Alunite LSU： Abundance of alunite endmembers obtained by LSU (ibid., 0–1).

[0046] m k m a These are the mean values ​​of the kaolinite abundance map and the alunite abundance map within the statistical domain (i.e., "Mask_TIR ∩ non-vegetation area"). They are used for scene standardization to eliminate differences in lighting / background.

[0047] , These are the standard deviations of the two abundance plots mentioned above within the same statistical domain. Together with the mean, they form the denominator of the Z-score standardization.

[0048] Output the MI-map, which only covers “Mask_TIR ∩ Non-vegetation Area”. The intersection of the two is used as the working area for LSU calculation and statistics of m, σ, avoiding interference from vegetation cover and non-thermal background.

[0049] S3: Combined constraints of tectonic landforms and hydrology (pixelation to optical resolution) S3.1 Tectonic Geomorphological Extraction: Using ALOS PALSAR imagery, edge detection filtering (such as Sobel, Canny, or directional gradient methods) is employed to identify linear tectonic zones, fault fracture zones, and geomorphic boundaries. The results are then rasterized and unified to the optical image resolution. The tectonic favorability index (SF) is obtained through a combination of fault density and intersection density.

[0050] S3.2 GRACE Hydrological Analysis: The NASA Gravity Recovery and Climate Experiment (GRACE) satellite was loaded onto the GEE platform. The study area was cropped and monthly series STL decomposition was performed to extract water storage trend slope and interannual amplitude information. The hydraulic favorability (HF) of total water storage (TWS) was obtained through weighted assignment.

[0051]

[0052] in A comprehensive index of "hydrological favorableness" constructed based on GRACE Total Water Storage (TWS); , , The weights are respectively the interannual (seasonal) amplitude, trend slope, and frequency of positive anomalies.

[0053] AMP: Interannual (Seasonal) Amplitude, the amplitude of the seasonal component (reflecting the intensity of replenishment-depletion within the year). The larger the amplitude, the stronger the replenishment / discharge cycle and the greater the active water volume; it is crucial for identifying aquifers in monsoon regions.

[0054] SLOPE: Trend slope, the linear slope of the total water storage trend over time. Unit / Value: cm / year; only positive slopes are considered (indicating an increase in storage over a multi-year timescale). A signal of continuous recharge / uplift; negative slopes often indicate depletion or the impact of extraction and are not included in the score.

[0055] PFREQ: Positive anomaly frequency, unit / value: dimensionless 0–1 (e.g., 0.67 if positive for 8 out of 12 months in a year). Measures the persistence of "more frequently in a wet / full state".

[0056] S4: Multi-factor fusion scoring and priority allocation, such as Figure 3 and 4 As shown.

[0057] S4.1 Pixel-level scoring: Perform multi-factor fusion calculation on all valid pixels within the “Mask_TIR” range. Calculate using the following formula:

[0058] in, Score: Represents the overall advantage. This represents the weight of the i-th factor.

[0059] This represents the normalized value of the i-th factor (standardized to the range of 0 to 1). The input factor data includes the thermal stability index H, vegetation stress index VIBS, mineral anomaly index MI, tectonic factor SF, and total water storage hydrological favorableness HF. TWS The result after 0-1 linear standardization.

[0060] S4.2 Patch Clustering and Grading: Threshold the fusion score, perform connected component clustering analysis, calculate the mean and 90th percentile (P90) of each patch, and classify the patches into high, medium and low priorities according to the quantile threshold (e.g. 80% / 50%), and output a candidate anomaly list and grade distribution map.

[0061] S5: Drone Inspection and Verification S5.1 Aerial Survey Planning: Cover only high-priority patches; flight altitude 120–300 m, side / heading overlap ≥70%; payload includes LWIR thermal imager (8–14 μm, blackbody calibrated) and visible / multispectral phase.

[0062] S5.2 Detailed Product Analysis: Generates centimeter-level temperature fields, vegetation stress textures, and gas anomaly trends; overlays with H, VIBS, MI, SF, and HF to refine boundary and causal explanations.

[0063] S5.3 Confirmation Rule: Under the constraints, it is confirmed as a genuine geothermal anomaly; otherwise, it is marked as a false anomaly or pending verification.

[0064] Example 2: In one exemplary embodiment, the Yuanjiang Hot Spring Pond (Yuanjiang County, Yunnan Province) experimental area was selected as an example. The study area was the Yuanjiang Hot Spring Pond area (where hot springs are concentrated, vegetation coverage is high, and the influence of topography and river valleys is significant) to verify the effectiveness of the present invention under complex cover conditions in subtropical mountainous areas.

[0065] Step 1: Initial screening of spaceborne thermal infrared temperature (delineating the candidate range) 1.1 Data and Regional Settings 1.1.1 Study Area and Coordinate System 1) Study Area: The experimental area is formed by setting out a 5-10 km buffer zone around the Yuanjiang Hot Spring Group as the core; 2) Coordinate system: The WGS 84 / UTM projection is used uniformly (the appropriate zone number is selected according to the longitude of the study area), and the elevation data adopts the SRTM 30 m digital elevation model (Shuttle Radar Topography Mission, SRTM).

[0066] 1.1.2 Image List and Time Period 1) Thermal infrared data: Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) night scene TIR (B10–B14, spatial resolution approximately 90m), selecting ≥6 scenes, covering both dry and wet seasons; 2) Alternative data: When ASTER lacks scenery or has excessive cloud cover in a certain season, supplement with Landsat-8 / 9 TIRS (Thermal Infrared Sensor, 100 m) night scenes; 3) Supporting data: Joint Research Centre Global Surface Water (JRC GSW), land use / cover (available from ESA WorldCover / domestic equivalent products), SRTM DEM.

[0067] 1.2 Pretreatment and Surface Temperature Inversion 1.2.1 Geometric and Radiation Preprocessing 1) Geometric correction: Using the highest quality ASTER scene as the master control, perform sub-pixel registration (RMSE ≤ 0.3 pixels). 2) Radiance calibration: Convert DN to radiance according to ASTER's official gain / bias settings; 3) Atmospheric correction: The FLAASH atmospheric correction model is parameterized using the radiative transfer model.

[0068] 1.2.2 Temperature-emissivity separation (TES) and land surface temperature (LST) 1) TES (Temperature Emissivity Separation, TES): Using the Temperature Emissivity Separation (TES) module of ENVI software, temperature-emissivity separation is performed on the ASTER TIR band. The wavelength unit is set to Micrometers, the minimum emissivity is 0.96, the number of iterations is 3, and the output is the surface temperature. LST, Kelvin ) and emissivity (ε) band; convert LST to Celsius temperature: LST(°C) = LST(K) – 273.15 2) Single-scene image thermal anomaly detection: for LST The background mean μ and standard deviation σ were calculated using a 3×3 sliding window.

[0069] Generate a preliminary thermal anomaly mask: ΔT ≥ + nσ Cells with n = 2.5 are marked as 1, and the rest are marked as 0.

[0070] in: This is a thermal anomaly. It is a temperature band, average value It provides a baseline value for the background temperature (the average temperature of the current pixel and its surrounding environment); the standard deviation σ provides a reference for the magnitude of temperature fluctuations, used to determine whether the current pixel is significantly higher than the background.

[0071] 3) Calculation of temporal stability index: The temporal stability index is obtained by performing pixel-by-pixel statistics on multi-scene thermal anomaly masks.

[0072]

[0073] in H Temporal stability index ΔT n This represents the number of times a pixel is in an abnormal state. N Represents the total number of images.

[0074] 4) Thermal anomaly mask output: extraction H ≥ 0.4 The pixels are used as stable thermal anomaly regions, and the thermal anomaly regions are output as binary masks (1 = thermal anomaly, 0 = non-anomaly).

[0075] 1.3 Suppression of water body and urban effects 1.3.1 Water Body Mask Since high temperatures in surface water may originate from thermal inertia and heat storage effects rather than geothermal activity, the initial infrared temperature screening stage shields the area affected by water bodies, prioritizing the identification of stable thermal anomalies on land. Hotspots near water bodies are then verified using subsequent ground-based methods.

[0076] Load the JRC GSW “Permanent Water” or other water body product (GeoTIFF). Use the ENVI-Toolbox-Raster-Reproject tool to unify the surface temperature (LST) mask projection and resolution. Generate the water body mask using Band Math, with pixels of value 1 representing water bodies.

[0077] 1.3.2 Urban Mask Load land use / cover data (e.g., GlobeLand30, ESA CCI Land Cover). Filter by "Urban Built-up Areas". Rasterize to a binary mask (Land Cover-mask) with consistent projection with the water body mask.

[0078] 1.4 Mask Blending and Output After unifying the water mask (water_mask) and the city mask (Land Cover-mask), the disturbance region (disturb_mask) is assigned a value of 1.

[0079] Using the Band Math tool, we calculated Mask_TIR = thermal_mask × (disturb_mask eq 0) to obtain a temperature mask that has removed interference from water bodies and cities.

[0080] The final output Mask_TIR is used as thermal anomaly mask data.

[0081] Step 2: Synergistic detection of vegetation and alteration mineral anomalies based on GF-5 hyperspectral imaging 2.1 Region Division and Preprocessing 2.1.1 Data and Correction 1) Data: Advanced Hyperspectral Imager (AHSI) / GF-5 (Gaofen-5); 2) Processing: Geometric correction → Radiometric calibration → Atmospheric correction (FLAASH / QUAC can be used) → Fine registration to ASTER. Bad band removal → FLAASH → MNF / inverse transform. 2.1.2 Vegetation / Non-vegetation Zoning 1) Calculate the Normalized Difference Vegetation Index (NDVI).

[0082] 2) Threshold: NDVI > 0.30 is defined as a vegetated area; otherwise, it is a non-vegetated area. 3) This step is only performed within the “Mask_TIR” scope.

[0083] 2.2 Vegetation Zone (VIBS: Vegetation Inhibition Index) 2.2.1 Index Calculation To improve robustness, an algorithm of "averaging three bands + first-order difference" is adopted. First, band preparation is performed: Red light neighborhood: R661, R671, R681 SWIR-2.27 µm: R2264, R2274, R2285 First-order difference with red edges: R70², R71², R73², R74² Calculate the average band using Band Math: = (R661 + R671 + R681) / 3 = (R2264 + R2274 + R2285) / 3 Red-edge first-order difference approximation: ΔR707 = R702 - R712 ΔR737 = R732 - R742 The geothermal suppression vegetation index is: 1) NDVI1 (SWIR long end elevation): NDVI 1 =( - ) / ( + ) NDVI2 (Red-edge blue shift): NDVI 2 =(ΔR707 - ΔR737) / (ΔR707 + ΔR737) Increased NDVI1: SWIR is brighter relative to red light – commonly seen in cases of decreased water content / changes in leaf structure / acidification stress. Increased NDVI2: The red edge peak shifts to shorter wavelengths (“blue shift”) – corresponding to decreased chlorophyll / growth inhibition.

[0084] 2.2.2 Standardize and synthesize VIBS (stress strength) 1) Scene-by-scene standardization (avoiding scene differences) Within the “vegetation mask”, calculate the mean (m) and standard deviation (std) of NDVI1 and NDVI2: m1 / std1, m2 / std2.

[0085]

[0086] VIBS is a dimensionless Z-fraction sum. Since Z1 and Z2 reflect independent stress information from different channels, the larger the value, the stronger the geothermal stress.

[0087] 3) Output: Vegetation stress intensity map (covering only "Mask_TIR ∩ Vegetation Zone").

[0088] 2.3 Non-vegetated areas (MI: Composite Alteration Index) 2.3.1 Linear Spectral Unmixing (LSU) 1) End-ingredients: Kaolinite and Alunite; 2) Spectral library: Laboratory endmembers resampled to the GF-5 band; 3) Output: Fraction (abundance) plots for Kaolinite and Alunite, along with RMS Error for quality control and removal of abnormal negative values / overflows.

[0089] 2.3.2 Index Synthesis

[0090] Where: MI: Overall Altering Index (MI-GF5), dimensionless. The larger the value, the stronger the acidic sulfate alteration (related to geothermal fluids).

[0091] Kaolinite LSU： The kaolinite endmember abundance obtained by linear spectral unmixing (LSU) represents the proportion of kaolinite in that pixel (0–1). The endmember spectra are resampled from laboratory / field spectral libraries to the actual GF-5 AHSI band before being used for unmixing.

[0092] Alunite LSU： Abundance of alunite endmembers obtained by LSU (ibid., 0–1).

[0093] m k m a These are the mean values ​​of the kaolinite abundance map and the alunite abundance map within the statistical domain (i.e., "Mask_TIR ∩ non-vegetation area"). They are used for scene standardization to eliminate differences in lighting / background.

[0094] , These are the standard deviations of the two abundance plots mentioned above within the same statistical domain. Together with the mean, they form the denominator of the Z-score standardization.

[0095] Step 3: Constructing combined geomorphological and hydrological constraints 3.1 Extraction of structural fracture zones (ALOS PALSAR) 3.1.1 Data and Augmentation 1) Data: Advanced Land Observing Satellite (ALOS-1 / 2) PALSAR (L-band, intensity map / coherence map); 2) Enhancement: Directional gradient / anisotropic filtering improves linear boundary and cliff response.

[0096] 3.1.2 Edge Detection and Linear Volume Recognition 1) Edge detection: Sobel / Canny / directional gradient (grouped by 0° / 45° / 90° / 135° direction); 2) Obtaining the construction line via the minimum cost path Constructing the Cost Raster: Based on edge enhancement results (such as Sobel) derived from NDVI, the Raster Calculator assigns high costs to non-line areas and low costs to line areas. Regions of interest are selected to construct raster cells for the fracture zone. Endpoint Determination: The raster is converted to polygons. Using vertex-to-point conversion, each vertex position of the line or polygon is extracted as a point feature, extracting all endpoints. Minimum Cost Path Calculation: The "Path Distance" tool takes endpoint data as input raster or feature source data, simulating fracture line connections along the strike. The cost can be set to 1 for cells on the line. Considering the impact of elevation undulation on path distance, a DEM is used as input surface raster, connecting nearby endpoints (500m). Finally, a backlink raster is output to record the direction encoding of the return from each cell to the source, resulting in a distance raster (CostDistance) and a backlink raster. The minimum cost path tool inputs endpoint data, the distance raster, and the backlink raster, and outputs the path raster as the construction feature.

[0097] 3) Geomorphological evidence: Verification of steep cliff / foothill lines and valley densities derived from superimposed DEM.

[0098] 3.1.3 Constructing Favorable Quantification (SF) 1) Fragment density FD (grid method): The neighborhood is statistically determined by the focal point (radius 1000 / 2000 m, statistical type is SUM) to obtain the number of pixels in the window × pixel length, which is then converted to line length / area (divided by window area) to obtain FD. FD is then normalized to 0-1.

[0099] 2) Intersection Density: Construct intersection kernel density (number of intersection points per unit area). Perform a direction calculation on the construction feature raster to obtain the direction value of each cell. Within a sliding window of 10*10, count the proportion of cells with large direction changes, which serve as the intersection / complex construction indicator layer.

[0100] 3) Structure Factor (SF): SF=w by weighting D × fracture density + w J × Intersection Density Among them W D and w J The weights are the fracture density and the intersection density, respectively.

[0101] 3.2 Regional hydrological background (GRACE / GRACE-FO) 3.2.1 Data and Decomposition 1) Data: Monthly raster of Terrestrial Water Storage (TWS) from the Gravity Recovery and Climate Experiment (GRACE / GRACE-FO) satellite; using the NASA / GRACE / MASS_GRIDS_V04 / LAND dataset, selecting the lwe_thickness_csr band (CSR solution result), where lwe_thickness is the equivalent liquid water thickness (cm), representing changes in water storage. This data is on a monthly scale with a resolution of approximately 1° (~111 km), reflecting changes in the entire water column (surface water + groundwater + soil water + snow and ice) relative to the reference period.

[0102] 2) Time-Series Decomposition: STL (Trend-Season-Residual) decomposition was performed on the TWS series of the study area and buffer zone. Multi-year variation curves were plotted, with the X-axis representing time and the Y-axis representing the equivalent liquid water thickness (cm). The curves reflect the long-term trend and seasonal fluctuations of water storage within the ROI: Peak period: High-water season (abundant precipitation and snowmelt). Low-water season: Low-water season (less precipitation and higher evaporation).

[0103] 3.2.2 Hydrological Favorability Hydrological favorableness serves only as a regional hydrological constraint layer in this study, used to characterize whether the water quantity conditions are favorable in the study area on a multi-year scale. The specific approach is as follows: (1) Indicator Construction Seasonal amplitude AMP (cm): Calculated as the monthly climatological average (12 amplitudes) of GRACE / GRACE-FO, representing the maximum to minimum of the year, indicating the seasonal strength of annual replenishment-consumption; a large amplitude indicates strong seasonal replenishment capacity such as rainy season / snowmelt.

[0104] Positive anomaly frequency PFREQ(0-1): The proportion of months with anomaly values ​​> 0 is calculated using the sequence of "monthly anomalies (current month value - climatological mean of the same month)", indicating a tendency for the weather to be frequently wet.

[0105] Multi-year trend SLOPE (cm / year): The Theil-Sen slope is calculated for anomalous sequences, reflecting long-term depletion / recovery (>0 indicates recovery, <0 indicates decline). The trend does not directly determine "favorable," but it can serve as a soft penalty to avoid giving high scores in long-term deficit areas.

[0106] (2) Normalization and Combination Linear normalize AMP to [0,1] (reference range 0-30 cm; can be relaxed to 40-50 cm in extreme monsoon / snowmelt areas), keep PFREQ to [0,1], and normalize positive slope 0-2 cm / year to [0,1].

[0107] Regional hydrological favorableness is defined as:

[0108] SLOPE only scores positive slopes; a trend soft penalty function f(S) is introduced: it takes the value 1 when S≥0; and linearly decreases to 0 when −2≤S<0, resulting in...

[0109] It emphasizes both seasonal replenishment and frequent periods of moisture, while avoiding misjudging high-favorability conditions in long-term arid areas. Note: This layer only represents the regional water volume background and is used as a constraint for subsequent superposition layers; it is not used for well-level location.

[0110] (3) Classification right Quantile classification is used: for example, Q20 / 60 / 80 divides the region into four levels: low, medium, relatively high, and high.

[0111] High-value areas (relatively high / high): These areas are characterized by strong seasonal replenishment and frequent humidity (with a non-negative trend), and can be used as the upper limit of the hydrological conditions for favorable geothermal areas. Low-value zone (medium / low): This zone has weak seasonality or is often dry, or has a long-term decline, and its weight should be reduced.

[0112] Step 4: Multi-factor fusion scoring and priority assignment 4.1 Pixel-level comprehensive score 4.1.1 Input and Alignment Input raster (continuous type, value range not uniform) H: Thermal anomaly stability index (0–1 or 0–100; the larger the value, the more stable). VIBS: Vegetation stress intensity MI: Composite Alteration Index SF: Structural Advantage Hydrological favorableness Mask: Mask_TIR (participates in fusion only within the time-stability region of thermal anomalies) 4.1.2 Unified scoring (standardization and consistency of direction) To avoid dimensional / scale differences, robust quantile standardization is used to the 0–1 interval. Quantiles are calculated for each indicator X: Q² = P²(X), Q₉₈ = P₉₈(X). The indicators are uniformly oriented as "higher = more favorable", namely H, VIBS, MI, SF, and HF. TWS The bigger it is, the better.

[0113] 4.1.3 Weight Adaptation

[0114] Score: Represents the overall advantage. Represents the weight of the i-th factor. This represents the normalized value of the i-th factor (standardized to the 0-1 range). The input factor data are the results of thermal stability index H, mineral anomaly index MI, etc., after 0-1 linear standardization.

[0115] 4.2 Patch Clustering and Grading (from Pixels to "Workable Target Areas") 1) Thresholding → Connected component clustering (minimum area ≥ 9 pixels) → Calculate patch mean and P90; 2) Divide into high / medium / low priority according to quantile threshold (e.g. 80% / 50%); 3) Perform morphological opening operation on small patches to denoise, forming a candidate list and rank map.

[0116] 4.3 Parameter Recommendations and Alternative Calibers Default weights (can be revised by region): H: 0.35 (thermal stability is the trigger condition); VIBS: 0.20 (strong forest indication); MI: 0.20 (mineral indication in bare land); SF: 0.15 (channel / converging structures); HF: 0.10 (recharge background). If the study area is mainly bare land, the weights of MI and SF can be increased; if the forest cover is heavy and acid gases are obvious, the weight of VIBS can be increased.

[0117] Step 5: Detailed inspection and ground verification of unmanned aerial vehicles (UAVs) 5.1 Aerial Survey Scheme (Covering only priority patches) Coverage: Only fly over High-level patches; if necessary, fly over adjacent (≤200 m) Medium patch boundary segments. Extend a buffer of 50–100 m outward from each patch to capture heat dissipation channels outside the boundary.

[0118] Time window and weather: 2–3 hours before dawn (to avoid solar heating and thermal hysteresis), no precipitation; wind speed ≤5m / s; relative humidity <85%.

[0119] For the same plaque, it is recommended to retest at least twice on different dates for stability assessment.

[0120] Flight and Resolution: Flight altitude is primarily 120–180 m (achieving a thermal resolution of approximately 10–20 cm GSD); large patches can be coarsely scanned at 250–300 m first, followed by low-altitude densification of hotspots. Overlap: Forward ≥80%, Lateral ≥70%; Ground speed 4–8 m / s (depending on camera integration time and texture). Flight Path: Prioritize contour lines / parallel paths, and ideally place a densification zone perpendicular to the main structural trend.

[0121] Load and calibration: LWIR thermal camera (8–14 μm, radiometric calibration type); calibration once before and after the factory double blackbody (or blackbody + high emissivity reference plate); high emissivity reference plate (ε≈0.95–0.98) and thermometer are placed on site as temperature scale constraints within the scene; 5.2 Data Processing and Thermal Anomaly Extraction (Temperature Evidence Only) Radiation-Temperature Conversion: Reads manufacturer-calibrated parameters to generate a temperature raster (°C) for each frame; outputs a thermal orthophoto mosaic (GeoTIFF, °C) and a DSM. Local background and difference temperature map are calculated on the thermal orthophoto using a circular neighborhood of r = 15–30 m. Interference elimination: Overlay road / building / aquaculture / industrial heat source layers; Water thermal inertia: Linear high temperatures ≤10–20 m immediately adjacent to the water surface and ditches are prioritized as "candidates for water body impact" and require ground verification. Shadows: Combine visible orthophotos and solar altitude angle to mark suspected shadow areas and reduce their credibility.

[0122] 5.3 Verification of Surface Gases and Surface Temperature (Targeting Aerial Survey Hotspots) Objective: To conduct rapid and low-cost on-site confirmation and interference source investigation of UAV thermal anomalies (hot spots); without conducting a comprehensive geochemical survey, only verifying "whether it is consistent with geothermal activity".

[0123] Principle: Only place data points in the hot spots (aerial survey P90 and hot spot center) within the High patch; record background-hot spot-comparison paired data.

[0124] Temperature: Contact thermometer / thermocouple (surface probe, ±0.5 °C); handheld thermal imager (e.g., 7.5–14 μm, with radiation thermometry function) for point scanning and confirmation; Gases: Portable H2S detector (range ≥ 50 ppm, resolution ≤ 0.1 ppm); Portable CO2 detector (0–50,000 ppm, resolution ≤ 10 ppm).

[0125] Example 3: Based on the same inventive concept, this application also provides a system for implementing the aforementioned multi-factor remote sensing detection system for geothermal water. The solution provided by this system is similar to the solution described in the above method. Therefore, the specific limitations of one or more embodiments of the multi-factor remote sensing detection system for geothermal water provided below can be found in the above-described limitations of a multi-factor remote sensing detection method for geothermal water, and will not be repeated here.

[0126] In one exemplary embodiment, such as Figure 2 As shown, a multi-factor remote sensing system for geothermal water is provided, comprising: Module 1: Data Access and Management (Unified collection of ASTER / Landsat, GF-5, ALOS, DEM, JRC, GRACE, UAV); Module 2: Preprocessing and Correction (Radiation / Atmosphere / Geometrics, TES, Effect Suppression); Module 3: Feature Extraction (Parallel by Element): Thermal Anomaly (H and Mask_TIR), Hyperspectral Vegetation Stress / Minerals (VIBS, MI), Tectonic Landforms (SF), Hydrological Field (HF) TWS ); Module Four: Multi-factor Fusion and Hierarchy; Module 5: Unmanned Aerial Vehicle (UAV) Inspection and Ground Verification; Module Six: Results Management and Publication.

[0127] In one exemplary embodiment, a computer device is also provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above-described method embodiments.

[0128] In one exemplary embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0129] In one exemplary embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0130] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments described above. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM).

[0131] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0132] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A geothermal water multi-factor remote sensing detection method, characterized in that, The geothermal water multi-factor remote sensing detection method comprises: ASTER and Landsat TIR images are acquired, and surface temperature is obtained through radiation calibration, atmospheric correction and TES inversion, combined with cloud, water body and urban heat island removal to generate a thermal stability anomaly map H and a mask Mask_TIR; In the mask area of the GF-5 hyperspectral data, the image is analyzed using the mask, the area is divided into a vegetation area and a non-vegetation area, for the vegetation area, red edge and SWIR index are calculated and standardized, and a vegetation stress index VIBS is synthesized; for the non-vegetation area, an end member unmixing method is used to extract the abundance of kaolinite and alunite, so as to construct an alteration index MI; The fault structure is extracted by ALOS PALSAR and the structure factor SF is formed, and the water storage trend, amplitude and frequency are obtained by combining with the GRACE TWS data to build the hydrological favorability HF TWS ; Final fusion of the H, VIBS, MI, SF, HF TWS Multi-factor, normalized weighted calculation of comprehensive score, clustering classification of candidate abnormal patches, verification by unmanned aerial vehicle thermal infrared and multispectral precision inspection, and output of accurate geothermal anomaly distribution results.

2. The method according to claim 1, wherein, The generation of the thermal stability anomaly map H and the mask Mask_TIR comprises: More than 6 scenes of ASTER TIR are acquired, and Landsat-8 / 9 TIRS is supplemented, radiation calibration, atmospheric correction and geometric precise registration are performed on each scene, and TES algorithm is used to invert the surface temperature LST; in the process of each scene image processing, cloud mask and JRC water body data set are applied to remove the interference caused by clouds and water body reflection in the image, in addition, through land use data identification and removal of obvious urban heat island effect, finally, the stable thermal anomaly map H and the temperature anomaly mask Mask_TIR without water body and urban interference are obtained.

3. The method according to claim 1, wherein, The construction of the alteration index MI specifically comprises: Regional division: ASTER and Landsat TIR images are acquired, and surface temperature is obtained through radiation calibration, atmospheric correction and TES inversion, combined with cloud, water body and urban heat island removal to generate a thermal stability anomaly map H and a mask Mask_TIR; ; In the mask area of the GF-5 hyperspectral data, the image is analyzed using the mask, the area is divided into a vegetation area and a non-vegetation area, for the vegetation area, red edge and SWIR index are calculated and standardized, and a vegetation stress index VIBS is synthesized; for the non-vegetation area, an end member unmixing method is used to extract the abundance of kaolinite and alunite, so as to construct an alteration index MI; The generation of the thermal stability anomaly map H and the mask Mask_TIR comprises: More than 6 scenes of ASTER TIR are acquired, and Landsat-8 / 9 TIRS is supplemented, radiation calibration, atmospheric correction and geometric precise registration are performed on each scene, and TES algorithm is used to invert the surface temperature LST; in the process of each scene image processing, cloud mask and JRC water body data set are applied to remove the interference caused by clouds and water body reflection in the image, in addition, through land use data identification and removal of obvious urban heat island effect, finally, the stable thermal anomaly map H and the temperature anomaly mask Mask_TIR without water body and urban interference are obtained. The construction of the alteration index MI specifically comprises: ; Regional division: ASTER and Landsat TIR images are acquired, and surface temperature is obtained through radiation calibration, atmospheric correction and TES inversion, combined with cloud, water body and urban heat island removal to generate a thermal stability anomaly map H and a mask Mask_TIR; ; In the mask area of the GF-5 hyperspectral data, the image is analyzed using the mask, the area is divided into a vegetation area and a non-vegetation area, for the vegetation area, red edge and SWIR index are calculated and standardized, and a vegetation stress index VIBS is synthesized; for the non-vegetation area, an end member unmixing method is used to extract the abundance of kaolinite and alunite, so as to construct an alteration index MI; The red edge blue shift index is then calculated monitoring the effect of geothermal water on vegetation; ； The generation of the thermal stability anomaly map H and the mask Mask_TIR comprises: ； wherein and are the mean and standard deviation, respectively; More than 6 scenes of ASTER TIR are acquired, and Landsat-8 / 9 TIRS is supplemented, radiation calibration, atmospheric correction and geometric precise registration are performed on each scene, and TES algorithm is used to invert the surface temperature LST; in the process of each scene image processing, cloud mask and JRC water body data set are applied to remove the interference caused by clouds and water body reflection in the image, in addition, through land use data identification and removal of obvious urban heat island effect, finally, the stable thermal anomaly map H and the temperature anomaly mask Mask_TIR without water body and urban interference are obtained. ； The construction of the alteration index MI specifically comprises: Regional division: Based on the obtained non-vegetation area mask, end member unmixing is performed on the hyperspectral image, kaolinite and alunite are selected as end members, the laboratory spectral data is resampled to the GF-5 spectral band, and the above-mentioned spectrum is used for unmixing processing, the abundance of kaolinite and alunite in each pixel is extracted, after the abundance is obtained, the mean value and standard deviation of the kaolinite and alunite abundance map in the non-vegetation area are calculated, and the Z-score standardization processing is performed to ensure that the differences in illumination and background in different regions are eliminated; The standardized abundance value is used to construct a comprehensive alteration index MI, and a mineral anomaly index map MI is output and overlaid on Mask_TIR intersection non-vegetation; ; Wherein, MI represents the comprehensive alteration index, dimensionless, the larger the value is, the stronger the acid sulfate type alteration is.

4. The method according to claim 1, wherein, The use of ALOS PALSAR extracts fault structure and forms structure factor SF, specifically includes: Using synthetic aperture radar ALOS PALSAR image, high-resolution radar data of the study area is obtained, image processing algorithm is used for edge detection of the radar image, geological structure features including fault zone and linear structure are identified; the structure line information detected is converted into a raster format for unified processing with other optical resolution data sets, and the structure information is rasterized to the same resolution as the optical image according to the spatial resolution of the image; in the generated raster data, the number of fault lines per unit area and the number of intersection points of different fault lines are calculated, which can be used to evaluate the complexity and activity of the regional structure; finally, a structure favorability map SF is generated, which is used to quantify the potential influence of the structure characteristics in the region on the geothermal anomaly.

5. The method according to claim 4, wherein, The constructing hydrological favorability HF TWS , specifically comprising: Data from the GRACE satellite is used to monitor changes in Earth's water storage, specifically acquiring GRACE TWS data. Time-series analysis is performed on the GRACE TWS data, applying the STL method to decompose the water storage data into seasonal variations, long-term trends, and random fluctuations. The slope of the long-term trend (SLOPE) is calculated; a positive value indicates an increase in water storage, while a negative value indicates a decrease. The slope is expressed in cm / year to reflect the dynamic changes in the hydrological system. The seasonal variation amplitude (AMP) is calculated, representing the difference between the maximum and minimum water storage values ​​within a year; a larger amplitude indicates a more active hydrological cycle. The frequency of positive anomalies (PFREQ) within a year is counted to assess the persistence of hydrological conditions and their environmental impact. This value is normalized to between 0 and 1, representing the frequency proportion. Based on set weights, SLOPE, AMP, and PFREQ are weighted and combined to derive a comprehensive hydrological favorability index, which reflects the supply capacity and abundance of water resources. Finally, a hydrological favorability map (HF) is generated. TWS It is used to assess the potential impact of hydrological factors on geothermal anomalies.

6. The method according to claim 1, wherein, The precise geothermal anomaly distribution result is obtained, specifically including: Multi-factor fusion scoring and priority division: Multi-factor fusion calculation is performed on all valid pixels within the Mask_TIR range, and the calculation is performed according to the following formula: ; wherein Score represents the comprehensive favorable degree, represents the weight of the i-th factor, represents the normalized value of the i-th factor, and the input factor data are thermal stability index H, vegetation stress index VIBS, mineral anomaly index MI, structural factor SF, and total water storage hydrological favorable degree HF TWS the result after 0-1 linear standardization processing; Threshold processing is performed on the fusion score Score, connected domain clustering analysis is performed, the mean value and 90th percentile of each patch are calculated, and the patches are divided into high, medium and low priority according to the percentile threshold, and a candidate anomaly list and grade distribution map are output; Unmanned aerial vehicle precision check and verification: based on the output high-priority patch list, the highest-priority patch is selected for unmanned aerial vehicle precision check, then UAV aerial survey is performed, that is, LWIR thermal image + multispectral / visible light, a centimeter-level temperature field, vegetation stress texture and gas anomaly trend are generated, H, VIBS, MI, SF and HF are superimposed, and high-precision geothermal anomaly verification results are output.

7. A geothermal water multi-factor remote sensing detection system, characterized in that, The system is used to execute the method of any one of claims 1-6, and the system comprises: Data access and management module: used for unified collection and management of multi-source remote sensing data, including ASTER and Landsat TIR images, GF-5 hyperspectral data, ALOS PALSAR images, DEM, JRC water body data, GRACE water storage data and UAV data; Pre-processing and correction module: used for a series of pre-processing on the acquired images, to ensure data availability, including: radiometric calibration, atmospheric correction and geometric rectification of ASTER and Landsat TIR images, using TES algorithm to retrieve land surface temperature, removing cloud, water body and urban heat island effect, generating thermal stability anomaly map H and temperature anomaly mask Mask_TIR; Feature extraction module: used to extract a variety of feature information from the processed image, including: thermal anomaly: generate thermal anomaly map H and mask Mask_TIR; hyperspectral vegetation stress / mineral: in the mask area of GF-5 hyperspectral data, analyze the vegetation area and non-vegetation area; calculate the red edge and SWIR index and standardize, synthesize the vegetation stress index VIBS, and extract kaolinite and alunite abundance in the non-vegetation area by end member unmixing to construct the alteration index MI; construct topography: use ALOS PALSAR to extract fault structure to form the structure factor SF; hydrological field: combine GRACE TWS data for time series analysis to extract water storage trend, amplitude and frequency to construct hydrological advantage degree HF TWS ; Multi-factor fusion and hierarchical module: used to fuse and analyze the extracted individual features: fusion of thermal anomaly map H, vegetation stress index VIBS, mineral anomaly index MI, structural factor SF and hydrological advantage degree HF TWS ; calculate the comprehensive score by normalized weighting, and perform clustering and grading to identify candidate anomaly patches; Unmanned aerial vehicle investigation and ground verification module, used for further verification of candidate anomaly patches: using unmanned aerial vehicle to carry thermal infrared and multispectral cameras to investigate high-priority patches; collecting and analyzing unmanned aerial vehicle data to generate accurate geothermal anomaly distribution results, comparing and confirming the authenticity of anomaly areas; Achievement management and release module, used for final processing and management of all obtained results, to ensure data reusability and user convenient access, outputting final geothermal anomaly distribution results, releasing together with related reports or documents, and sharing key data and information with the outside world.

8. A computer device comprising: A memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that the processor executes the computer program to implement the steps of the geothermal water multi-factor remote sensing detection method of any one of claims 1-5.

9. A computer readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to implement the steps of the geothermal water multi-factor remote sensing detection method of any one of claims 1-5.

10. A computer program product comprising a computer program, characterized in that, The computer program is executed by the processor to implement the steps of the geothermal water multi-factor remote sensing detection method of any one of claims 1-5.

Images