A method for extracting surface geothermal anomaly information under low solar radiation

Abstract

The present invention discloses a method for extracting surface geothermal anomaly information under low solar radiation. The method selects remote sensing data from the winter solstice period to extract geothermal anomaly information, uses a single-window algorithm to invert the surface temperature to form time series data, and calculates the variance pixel by pixel to reflect the degree of dispersion of the surface temperature. When the dispersion is lower, the constant temperature effect produced by the deep heat source of the earth is more obvious, thereby extracting geothermal anomaly information.

Description

Technical Field

[0001] The present invention relates to the field of surface analysis technology, and more particularly to a method for extracting surface geothermal anomaly information under low solar radiation. Background Art

[0002] The Earth's surface temperature is primarily derived from solar radiation, followed by deep-Earth heat sources. The former, which transmits heat in the form of electromagnetic radiation, plays a dominant role in warming the Earth's surface and is present everywhere. The latter, geothermal resources, transmit heat through conduction and convection, primarily controlled by tectonic forces and influenced by the physical properties of strata and rocks. This manifests as localized warming, creating thermal anomalies on the surface and providing a direct clue to the search for geothermal resources.

[0003] With the recent advancement of remote sensing technology, many researchers have conducted research on surface temperature inversion in geothermal anomaly areas using thermal infrared remote sensing, achieving promising results. Most existing studies focus on validating surface temperature inversion algorithms, or extracting or predicting geothermal anomaly information in the study area through cross-comparison of algorithms or other auxiliary methods. However, these research methods are mostly based on a single scene and a single phase, ignoring the "homeothermal effect" of geothermal heat as an endogenous energy source. This means that the heat provided by the deep Earth to the surface remains relatively stable regardless of changes in the external environment. However, due to factors such as heat source type, depth, geological structure, heat migration, and thermal reservoir caprock, thermal anomalies on the surface are often not obvious, making it difficult to accurately extract surface temperature from a single period of thermal infrared data.

[0004] Therefore, how to accurately extract surface geothermal anomaly information is an urgent problem that needs to be solved by those skilled in the art. Summary of the Invention

[0005] In light of this, the present invention provides a method for extracting surface geothermal anomalies under low solar radiation conditions. This method uses multiple Landsat 8 data periods during the winter solstice to minimize the surface heat source from solar radiation. Surface temperature time series data are then inverted and the variance is calculated to reflect the dispersion of surface temperature. Lower dispersion indicates a more pronounced "constant temperature effect," allowing for the extraction of geothermal anomaly information. The period around the winter solstice is when the sun's altitude is lowest and the ground receives the least solar radiation. During this time, deep Earth heat sources contribute significantly to surface temperature. Selecting this period as the optimal time window for geothermal anomaly extraction improves its accuracy.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] A method for extracting surface geothermal anomaly information under low solar radiation includes the following specific steps:

[0008] Step 1: Collect Landsat8 remote sensing data of the winter solstice period in the area to be extracted according to work needs; Landsat8 remote sensing data includes OLI data and TIRS thermal infrared data;

[0009] Step 2: Perform radiation correction on the remote sensing data of each period to obtain the corrected data;

[0010] Step 3: Use the single window algorithm to perform surface temperature inversion on the corrected data to obtain the surface temperature;

[0011] Step 4: Superimpose the surface temperatures inverted in each period to obtain surface temperature time series data;

[0012] Step 5: Calculate the variance of the surface temperature time series data pixel by pixel to obtain surface temperature dispersion data;

[0013] Step 6: Select OLI data from a period for atmospheric correction and use the spectral angle (SAM) method to extract interfering objects as mask data; the interfering objects include vegetation, water bodies, clouds, and shadows. Select OLI data from a period with less interfering objects such as less clouds and snow;

[0014] Step 7: Masking and filtering the surface temperature dispersion data according to the mask data to obtain processed surface temperature dispersion data, thereby achieving the purpose of interference removal and data optimization;

[0015] Step 8: Threshold cutting is performed on the processed surface temperature dispersion data according to the set cutting threshold, and the portion of the processed surface temperature dispersion data that is less than the threshold is used as geothermal anomaly information. The geothermal anomaly area is obtained based on its inherent spatial coordinate attributes.

[0016] Preferably, the formula for surface temperature inversion using a single window algorithm in step 3 is:

[0017] T s =[a(1-CD)+(b(1-CD)+C+D)T b -DT a ] / C (1)

[0018] C=ετ (2)

[0019] D=(1-τ)[1+τ(1-ε)] (3)

[0020] Among them, T s is the surface temperature; a and b are set constants; ε is the surface emissivity; τ is the atmospheric transmittance; T a is the average atmospheric temperature; Tb is the brightness temperature.

[0021] Preferably, the surface emissivity needs to be estimated and determined. When estimating, the surface structures of different regions of the earth's surface are divided into water bodies, vegetation, bare land, and built-up land. The emissivity of vegetation is 0.986, and that of water bodies is 0.995. The estimated expression is:

[0022] When NDVI < 0.05, it indicates that the water body is dominant, ε = 0.995;

[0023] When NDVI>0.72, it indicates that vegetation is dominant, ε=0.986;

[0024] When 0.05≤NDVI≤0.72, it means that the land is mainly bare land or construction land.

[0025] ε=0.985×P v +0.960×(1-P v )+0.06×P v ×(1-P v ) (4)

[0026] in,

[0027] Where NDVI is the normalized difference vegetation index; R represents the red band, which corresponds to the 4th band in the OLI data; NIR represents the near infrared band, which corresponds to the 5th band in the OLI data;

[0028]

[0029] Where, P v is the proportion of vegetation in the pixel, NDVI s and NDVI v They are the bare land NDVI threshold and the vegetation NDVI threshold, respectively.

[0030] Preferably, the variance calculation formula in step 5 is:

[0031]

[0032] Where x is the surface temperature value inverted for a period corresponding to the pixel of the time series data; M is the average surface temperature corresponding to the pixel; n is the number of periods; S 2 is the variance.

[0033] Preferably, the surface temperature dispersion data includes a surface temperature dispersion map. The surface temperature dispersion map is obtained by calculating the deviation of the surface temperature from the mean value for each period based on the variance, and summing and averaging the results using a square method. A smaller variance indicates more clustered data and lower dispersion; a larger variance indicates more discrete data and higher dispersion.

[0034] Preferably, the spectral angle method in step 6 includes:

[0035] Step 61: Calculate the spectral angle of the OLI data using the formula:

[0036]

[0037] Where α is the spectral angle; d is the reference interference object spectrum; y is the image spectrum in the OLI data;

[0038] Step 62: The OLI data corresponding to the calculated spectral angle less than the set spectral angle threshold is determined as an interfering object and used as mask data. The reference interfering object spectrum is obtained by visual interpretation to establish clouds, vegetation, water bodies, and shadows.

[0039] Preferably, the masking process is to delete the area corresponding to the interfering ground objects on the surface temperature dispersion map.

[0040] Preferably, the surface temperature dispersion data conforms to a normal distribution, the standard deviation of the surface temperature dispersion data is used as a cutting threshold, and the surface temperature area information with a temperature less than the cutting threshold is used as geothermal anomaly information.

[0041] Through the above technical solutions, it can be seen that compared with the existing technology, the present invention discloses a method for extracting surface geothermal anomaly information under low solar radiation. The period around the winter solstice is the period with the least solar radiation throughout the year. At this time, the contribution of heat sources deep in the earth to the surface temperature is relatively high. The thermal infrared data of this period is selected as the optimal time window for extracting geothermal anomaly information; Landsat8 remote sensing data with similar time periods is used to ensure that each period of data has similar climate characteristics, solar radiation and other background environments, and is comparable. Compared with single-period data, multi-period remote sensing data has higher reliability, further improving the accuracy of surface geothermal anomaly information extraction; a single-window algorithm is used to invert the surface temperature to form time series data, and the variance is calculated pixel by pixel to reflect the dispersion of the surface temperature. When the dispersion is lower, the constant temperature effect is more obvious, thereby extracting geothermal anomaly information. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are merely embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying any creative work.

[0043] Figure 1 The accompanying drawing is a schematic flow chart of a method for extracting surface geothermal anomaly information under low solar radiation provided by the present invention;

[0044] Figure 2 The accompanying drawing is a schematic diagram of the surface temperature dispersion degree diagram provided by the present invention;

[0045] Figure 3 The accompanying drawing is a schematic diagram of the distribution of the extracted abnormal areas provided by the present invention. DETAILED DESCRIPTION

[0046] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0047] The embodiment of the present invention discloses a method for extracting ground thermal anomaly information under low solar radiation. Figure 1 Schematic diagram of the process shown.

[0048] Example

[0049] The method of the present invention is described by taking the Taxkorgan Valley as an example.

[0050] S1: Using Landsat8 OLI / TIRS remote sensing data, the imaging time was selected near the winter solstice and there was no cloud or snow cover to obtain image data passing through the study area. The time was December 15, 2014, December 2, 2015, December 20, 2016, December 29, 2019, and December 31, 2020, a total of five periods;

[0051] S2: Perform radiation correction on the remote sensing data collected in each period to obtain the corrected data;

[0052] S3: Surface temperature is inverted using a single window algorithm on the corrected data of each period to obtain the surface temperature;

[0053] Due to calibration instability in the Landsat 8 TIRS Band 11, the split-window algorithm is not recommended for quantitative research. The single-window algorithm directly incorporates the effects of the atmosphere and surface into the calculation formula and is simpler and easier to use than the radiative transfer equation method and the single-channel algorithm.

[0054] The formula is as follows:

[0055] T s =[a(1-CD)+(b(1-CD)+C+D)T b -DT a ] / C (1)

[0056] C=ετ (2)

[0057] D=(1-τ)[1+τ(1-ε)] (3)

[0058] Where, T s is the surface temperature (K), a and b are constants, ε is the surface emissivity, τ is the atmospheric transmittance, T a is the average atmospheric temperature (K), T b is the brightness temperature;

[0059] (1) Atmospheric transmittance τ

[0060] Atmospheric transmittance can be obtained by querying the NASA official website and entering the imaging time and the center longitude and latitude coordinates;

[0061] (2) Average atmospheric temperature T a

[0062] There are four estimation equations for the atmospheric mean effective temperature of standard atmospheres. The historical temperature T0 (converted to Kelvin) of the study area can be queried on the website and the atmospheric mean effective temperature can be estimated using Table 1. In this embodiment, the atmospheric mean effective temperature of the mid-latitude winter mode is selected for estimation.

[0063] Table 1 Estimation method of atmospheric average effective temperature

[0064] Atmospheric model Estimation equation for atmospheric mean effective temperature tropical atmosphere Ta=17.9769+0.91715T0 mid-latitude summer Ta=16.0110+0.92621T0 mid-latitude winter Ta=19.2704+0.91118T0 1976 U.S. Standard Atmosphere Ta=25.9396+0.88045T0

[0065] (3) Surface emissivity ε

[0066] The surface emissivity estimation method is selected. It is considered that although the surface structure of different regions of the earth's surface is complex, it can be roughly divided into water bodies, vegetation, bare land and built-up land. The emissivity of vegetation is 0.986 and that of water bodies is 0.995. The expressions are as follows:

[0067] ① When NDVI < 0.05, it indicates that the water body is dominant, ε = 0.995;

[0068] ② When NDVI>0.72, it indicates that vegetation is dominant, ε=0.986;

[0069] ③ When 0.05≤NDVI≤0.72, it means that the land is mainly bare land or construction land.

[0070] ε=0.985×P v +0.960×(1-P v )+0.06×P v ×(1-P v ) (4)

[0071] in,

[0072] Where NDVI is the normalized difference vegetation index, R represents the red band, and NIR represents the near-infrared band, which correspond to bands 4 and 5 of Landsat 8 OLI, respectively.

[0073]

[0074] Where, P v is the proportion of vegetation in the pixel, NDVI s and NDVI v are the NDVI thresholds for bare land and vegetation, respectively;

[0075] S4: superimpose the five phases of surface temperature inversion data to form time series data;

[0076] S5: Calculate the variance of each pixel, count the deviation of each period of data from the mean, and use the square method to sum and average to avoid the mutual cancellation of positive and negative numbers, thereby obtaining the surface temperature dispersion map of the study area, such as Figure 2 As shown in the figure, when the variance is smaller, the data is more clustered and the dispersion is lower; when the variance is larger, the data is more discrete and the dispersion is higher.

[0077]

[0078] Where x is the surface temperature value of a certain period inverted by the pixel of the time series data; M is the average surface temperature corresponding to the pixel; n is the number of periods; S 2 is the variance; Based on this method, the geothermal anomaly information in the area of ​​Quman Village, Taxkorgan County is extracted;

[0079] S6: The ground surface is blocked by shadows, clouds, snow, vegetation, etc., and solar radiation cannot reach directly, which affects the ground heating. These interference factors need to be eliminated. The spectral angle (SAM) method is used to extract interfering ground objects and use them as masks.

[0080] The spectral angle method is used to calculate the similarity between the test spectrum and the reference spectrum. The calculation result can be regarded as the cosine angle between the two arrays. When the spectral angle a value is smaller, it means that the two sets of spectra are more matched and the similarity is higher. Based on this, the above-mentioned interference objects can be extracted;

[0081] S61: Calculate the spectral angle of OLI data using the formula:

[0082]

[0083] Where α is the spectral angle; d is the reference interference object spectrum; y is the image spectrum in the OLI data;

[0084] S62: The OLI data corresponding to the calculated spectral angle less than the set spectral angle threshold is judged as an interfering object and used as mask data; cloud bodies, vegetation, water bodies, and shadows are established through visual interpretation to obtain a reference interfering object spectrum;

[0085] S7: Mask the surface temperature dispersion data to reduce the impact of interfering objects on the extraction of geothermal anomaly information. Mask processing means deleting the area corresponding to the mask data on the surface temperature dispersion map based on the spatial coordinate attributes of the OLI data. Filter the surface temperature dispersion data. Isolated and fragmented patches often appear on the image. Filtering is used to remove isolated pixels or merge them into surrounding patches to ensure that the encircled anomaly patches have good continuity and regularity. It is recommended to use a 3×3 pixel window for mean filtering.

[0086] S8: According to the principles of geostatistics, the surface temperature dispersion map conforms to the normal distribution. The standard deviation is used as the reference value for anomaly cutting, and values ​​less than this reference value are regarded as geothermal anomaly information.

[0087] By analyzing geothermal anomaly information, a total of 5 anomaly areas with certain scales were obtained, such as Figure 3 As shown, the anomalous areas are primarily distributed along the western bank of the Tashkurgan River north of Quman Village and in the nearby piedmont alluvial plain. Field verification and exploration projects were conducted on geothermal anomaly areas 1, 2, 3, and 4, revealing the presence of geothermal bodies. Anomalies 1 and 2 can be connected as a single entity, as can anomalies 3 and 4. Area 5 has not been verified. The characteristics are shown in Table 2.

[0088] Table 2 Field verification of geothermal anomaly areas

[0089]

[0090] The distribution of geothermal anomalies in the study area is controlled by fault structures, and the spatial distribution of the two is consistent. Five geothermal anomaly areas were discovered, four of which are consistent with the field situation, and one requires further field verification. These results indicate that the proposed method has certain reliability and reference value for extracting geothermal anomaly information.

[0091] Beneficial effects of the present invention:

[0092] 1) The period around the winter solstice is the period with the least solar radiation throughout the year. During this period, the contribution of heat sources deep in the earth to the surface temperature is relatively high. Therefore, the thermal infrared data of this period is selected as the best time window for extracting geothermal anomaly information.

[0093] 2) Using multiple Landsat 8 remote sensing data sets from similar time periods ensures that each data set has similar climate characteristics, solar radiation, and other background conditions, making it comparable. Compared to single-period data, multi-period remote sensing data is more reliable and further improves the accuracy of surface geothermal anomaly information extraction.

[0094] 3) Five phases of Landsat8 OLI / TIRS data were acquired and surface temperature inversion was performed using a single-window algorithm to form time series data. The variance was calculated pixel by pixel to reflect the dispersion of surface temperature. The lower the dispersion, the more obvious the "isothermal effect", thus extracting geothermal anomaly information.

[0095] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. Reference can be made to the common and similar parts between the various embodiments. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and the relevant parts can be referred to the method description.

[0096] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be readily apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not limited to the embodiments shown herein but is intended to conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for extracting surface geothermal anomaly information under low solar radiation, characterized in that: The specific steps include: Step 1: Collect Landsat8 remote sensing data of the winter solstice period in the area to be extracted. Landsat8 remote sensing data includes OLI data and TIRS thermal infrared data; Step 2: Perform radiation correction on the remote sensing data of each period to obtain the corrected data; Step 3: Use the single window algorithm to perform surface temperature inversion on the corrected data to obtain the surface temperature; Step 4: Superimpose the surface temperatures inverted in each period to obtain surface temperature time series data; Step 5: Calculate the variance of the surface temperature time series data pixel by pixel to obtain surface temperature dispersion data; calculate the deviation of the surface temperature from the mean value in each period based on the variance, and sum and average them using a square method to obtain a surface temperature dispersion degree map; Step 6: Select a period of OLI data for atmospheric correction and use the spectral angle method to extract interfering objects as mask data; the spectral angle method includes: Step 61: Calculate the spectral angle of the OLI data using the formula: Where α is the spectral angle; d is the reference interference object spectrum; y is the image spectrum in the OLI data; Step 62: Determine the OLI data corresponding to the spectral angle being smaller than the set spectral angle threshold as interference objects and use the data as mask data; Step 7: performing mask processing and filtering processing on the surface temperature dispersion data according to the mask data to obtain processed surface temperature dispersion data; Step 8: Perform threshold cutting on the processed surface temperature dispersion data according to the set cutting threshold, and take the part smaller than the threshold as geothermal anomaly information.

2. The method for extracting surface geothermal anomaly information under low solar radiation according to claim 1, characterized in that: The formula for surface temperature inversion using the single window algorithm in step 3 is: T s =[a(1-C-D)+(b(1-C-D)+C+D)T b -DT a ] / C (1) C=ετ (2) D=(1-τ)[1+τ(1-ε)] (3) Among them, T s is the surface temperature; a and b are set constants; ε is the surface emissivity; τ is the atmospheric transmittance; T a is the average atmospheric temperature; T b is the brightness temperature.

3. The method for extracting surface geothermal anomaly information under low solar radiation according to claim 2, characterized in that: The surface emissivity needs to be estimated and determined. When estimating, the surface structure of different areas of the earth's surface is divided into water bodies, vegetation, bare land and building land. The emissivity of vegetation is 0.986 and that of water bodies is 0.

995. The estimated expression is: When NDVI < 0.05, it indicates that the water body is dominant, ε = 0.995; When NDVI>0.72, it indicates that vegetation is dominant, ε=0.986; When 0.05≤NDVI≤0.72, it means that the land is mainly bare land or construction land. ε=0.985×P v +0.960×(1-P v )+0.06×P v ×(1-P v ) (4) in, Where NDVI is the normalized difference vegetation index; R represents the red band, which corresponds to the 4th band in the OLI data; NIR represents the near infrared band, which corresponds to the 5th band in the OLI data; Where, P v is the proportion of vegetation in the pixel, NDVI s and NDVI v They are the bare land NDVI threshold and the vegetation NDVI threshold, respectively.

4. The method for extracting surface geothermal anomaly information under low solar radiation according to claim 1, characterized in that: The variance calculation formula in step 5 is: Where x is the surface temperature value inverted for a period corresponding to the pixel of the time series data; M is the average surface temperature corresponding to the pixel; n is the number of periods; S 2 is the variance.

5. The method for extracting surface geothermal anomaly information under low solar radiation according to claim 1, characterized in that: The surface temperature dispersion data conforms to the normal distribution, and the standard deviation corresponding to the surface temperature dispersion data is used as the cutting threshold.

Images