A time-series temperature anomaly remote sensing extraction method for coal fire area identification

Abstract

The application discloses a time-series temperature anomaly remote sensing extraction method for coal fire area identification and belongs to the field of remote sensing algorithm development. Optical image data with a thermal infrared band are acquired, radiation calibration and atmospheric correction are performed on the optical image data, NDVI, vegetation coverage VFC and surface emissivity epsilon are calculated, atmospheric profile parameters are acquired, a temperature inversion method is selected to invert the surface temperature, an adaptive window-based temperature anomaly extraction algorithm is used to extract temperature anomaly pixels from the inverted surface temperature data, meanwhile, the time continuity of coal fire combustion is considered, a time parameter is introduced, time-series temperature anomaly pixels are extracted through time-series analysis, and a time-series temperature anomaly image is composed of the time-series temperature anomaly pixels. The algorithm based on unbiased estimation is simple to operate, fast in operation, avoids introduction of artificial subjective errors and can obviously improve the accuracy of extraction of temperature anomaly regions.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of remote sensing algorithm development, in particular to a time series temperature anomaly remote sensing extraction method for coal fire area identification. TECHNICAL BACKGROUND

[0002] At present, there are five kinds of methods for coal fire detection and monitoring, including geophysical prospecting, geochemical prospecting, thermal prospecting, drilling and remote sensing. The first four methods are commonly used and effective, but they need field research, consume a lot of manpower and financial resources, and it is difficult to achieve for large fire areas. Remote sensing technology has the advantages of large-area synchronous observation, traceability and economy, and has gradually become one of the mainstream methods for coal fire detection in recent years. The heat generated by coal fire combustion will be transferred to the ground, causing the ground temperature in the coal fire area to be higher than that in the non-coal fire area. Using thermal infrared remote sensing images to retrieve ground temperature and extract temperature anomaly areas to identify the scope of the fire area is the main means of remote sensing coal fire detection, and the key and difficulty is how to select the appropriate temperature threshold to extract the temperature anomaly area.

[0003] The current commonly used temperature anomaly extraction methods mainly include two categories: single threshold algorithm and multi-threshold algorithm; the single threshold algorithm mainly includes fixed threshold algorithm (FTA), gradient threshold method, adaptive gradient threshold algorithm and space-time temperature threshold algorithm; the multi-threshold algorithm is mainly a sliding window algorithm based on unbiased estimation (Temperature Anomaly Extraction, TAE). The single threshold algorithm selects a fixed threshold to extract the temperature anomaly region in the whole image, which will lead to the loss or surplus of the temperature anomaly region, and has subjectivity. The sliding window algorithm based on unbiased estimation can effectively improve the disadvantages of the single threshold algorithm by selecting a local threshold to extract the temperature anomaly region, but this method lacks stability, which is shown in the following aspects: (1) the temperature anomaly extraction result is greatly affected by the window size; (2) the standard for selecting the temperature anomaly pixel in the window is single, which affects the recognition accuracy. In addition, the existing algorithms all identify the temperature anomaly region by using independent images, and do not fully consider that the coal fire combustion is a continuous process in time, so that the temperature anomaly recognition accuracy of the fire area is reduced. In view of the above problems, the present application provides a time-series temperature anomaly remote sensing extraction method for coal fire area identification (Temperature Anomaly of Time-Series Extraction Algorithm based on Adaptive Windows, TATSAW), which comprises a temperature anomaly extraction method based on adaptive window (Temperature Anomaly Extraction Based On Adaptive Window, TAEAW) and a time-series processing method. The adaptive window is designed to establish multiple standards in the window to screen the temperature anomaly pixels; at the same time, the time continuity of the coal fire combustion is considered, the time parameter is introduced, and the time-series temperature anomaly information is extracted through time-series analysis. SUMMARY

[0004] In view of the problems existing in the prior art, the present application provides a time-series temperature anomaly remote sensing extraction method for coal fire area identification, which adopts an adaptive window to establish multiple standards in the window to screen the temperature anomaly pixels; at the same time, the time continuity of the coal fire combustion is considered, the time parameter is introduced, and the time-series temperature anomaly information is extracted, so as to overcome the shortcomings of the traditional temperature anomaly extraction method and more accurately extract the temperature anomaly region of the coal fire area.

[0005] To achieve the above object, the present application provides a time-series temperature anomaly remote sensing extraction method for coal fire area identification, characterized in that the steps are as follows:

[0006] Select as many q amplitude Landsat 8 image data as possible in the research time range with less than 10% cloud cover; Extract the Landsat 8 image data with the thermal infrared band in the measured area, perform radiation calibration and atmospheric correction on the original optical image data, calculate the normalized vegetation index NDVI, vegetation cover VFC and surface emissivity ε, and obtain the atmospheric profile parameters;

[0007] Use the single window algorithm to perform temperature inversion on the obtained Landsat 8 image data to obtain temperature inversion information;

[0008] Use adaptive window and multi-threshold to extract abnormal temperature information in the temperature inversion information;

[0009] Considering the time continuity of coal fire combustion, introduce the time parameter of coal fire combustion, determine the time frequency threshold of coal fire pixels, and extract the time series temperature anomaly information according to the time frequency threshold, so as to determine the coal fire area.

[0010] Further, import the image into ENVI software to correct the multispectral band and thermal infrared band in turn, perform atmospheric correction on the multispectral band to obtain the reflectivity image of the visible band, and obtain the atmospheric transmittance τ of the image on the research area through the NASA atmospheric profile parameter calculation website.

[0011] Further, the temperature inversion process of the single window algorithm is as follows:

[0012] Calculate the surface emissivity of the q amplitude selected Landsat 8 image data in turn:

[0013] Calculate NDVI: input the calculation formula of normalized vegetation index (float(b4-b3)) / (b4+b3) in ENVI software bandmath, where b3 and b4 are the red band and near-infrared band after atmospheric correction, respectively. Due to atmospheric correction, the value of NDVI will exceed [-1, 1], and the method of removing outliers is used to remove the corrected NDVI;

[0014] Calculate the vegetation cover FVC using the formula Calculate, where NDVIS and NDVIV represent the minimum and maximum values of NDVI in a certain confidence interval, usually 5% and 95%;

[0015] Calculate the surface emissivity LSE. According to different ground objects, their emissivity is different. Generally, remote sensing ground objects are divided into water, town and natural surface; the emissivity calculation formula of three different types of ground object pixels is:

[0016] E(water) = 0.995, E(building) = 0.9589 + 0.086*FVC - 0.0671*FVC*FVC, E(building) = 0.9589 + 0.086*FVC - 0.0671*FVC*FVC;

[0017] Input the surface emissivity calculation formula in ENVI software bandmath:

[0018] (b1le NDVIS)*0.995+(b1gt NDVIS and b1lt NDVIV)*(0.9589+0.086*b2-0.0671*b2*b2)(b1ge NDVIV)*(0.9625+0.0614*b2-0.0461*b2*b2),

[0019] Wherein, b1 is NDVI, b2 is vegetation coverage FVC, and the surface emissivity LSE is obtained.

[0020] Inversion of the study area surface temperature: first, use the formula C = LSE * tau, D = (1-t) * (1+(1-LSE) * tau) to calculate the intermediate variables C, D; LSE is the surface emissivity, tau is the atmospheric transmittance;

[0021] Calculate the on-board brightness temperature T6, input the calculation formula in ENVI software bandmath:

[0022]

[0023] b1 selects the radiometrically corrected thermal infrared band, and the result is T6;

[0024] Inversion of the study area surface temperature,

[0025] Input the formula (-67.355351*(1-b1-b2)+(0.458606*(1-b1-b2)+b1+b2)*b3-b2*T a ) / b1 in bandmath;

[0026] In the formula, b1, b2, b3 represent intermediate parameters C, D, T6, T a is the Kelvin temperature of the study area on the day, and the calculation result obtained through the atmospheric profile parameters is the surface temperature inverted by the single-window algorithm.

[0027] Further, the steps of using adaptive window and multi-threshold to extract abnormal temperature information in Landsat 8 image temperature inversion information are as follows:

[0028] Initial window selection: the initial adaptive window size should meet the minimum resolution of remote sensing image, the resolution of thermal infrared band of Landsat 8 image data is 100m, the pixel size is 30m x 30m, therefore the initial window size is selected as 4 x 4;

[0029] Adaptive window expansion criterion: extract Landsat 8 image, calculate the mean value m of the whole image;

[0030] The adaptive window slides from the top left corner of the image, and the mean value m0 of the initial window is calculated,

[0031] If m0≤m, the window is no longer expanded;

[0032] If m0>m, the window is expanded to the four directions with a step of 1, if the image boundary is encountered, the window is expanded to the non-boundary of the image with a step of 1, and the mean value m of the newly expanded window is calculated n (n = 1, 2, 3, …),

[0033] If m n ≤m0, the window is no longer expanded;

[0034] If m n >m0, the window is continuously expanded with a step of 1 until m n ≤m0, the window is no longer expanded;

[0035] Multi-threshold segmentation in adaptive window: calculate the mean value m n (n = 0, 1, 2, …) and the standard deviation std in the window, take m n +std as the threshold t1, and take m n +2*std as the threshold t2;

[0036] When the pixel value value in the window satisfies value>m n +2*std, the pixel is selected as a thermal anomaly pixel, and the intermediate variable M of the total number of times that a certain pixel is selected as a thermal anomaly pixel is increased by 1 based on the original value;

[0037] When the pixel value value in the window satisfies m n +std<value≤m n +2*std, the pixel is selected as a low-frequency thermal anomaly pixel, and the intermediate variable M of the total number of times that a certain pixel is selected as a thermal anomaly pixel is increased by 0.5 based on the original value;

[0038] When the pixel value value in the window satisfies value≤m n +std, the pixel is not considered as a thermal anomaly pixel, and the intermediate variable M of the total number of times that a certain pixel is selected as a thermal anomaly pixel is increased by 0 based on the original value.

[0039] Further, the temperature anomaly information extraction method is as follows:

[0040] The temperature anomaly information is extracted by extracting the temperature anomaly pixels in the Landsat 8 image: each pixel in the Landsat 8 image has intermediate variables N and M, the total number of times that the pixel is selected into the adaptive window N and the total number of times that the pixel is selected as a thermal anomaly pixel M, the initial N = 0, M = 0,

[0041] The temperature anomaly pixel is extracted by relying on two intermediate variables, if Then the pixel is marked as a thermal anomaly pixel; the calculation of the intermediate variable N is based on the adaptive window expansion criterion, and the pixel is selected by the initial window or the expanded window, which will be added to the original intermediate variable N; the calculation of the intermediate variable M is based on the multi-threshold segmentation criterion in the adaptive window;

[0042] Run from the upper left of the Landsat 8 image, calculate the intermediate variables N and M, calculate the intermediate variables N and M again with a step of 4, and so on, until the image is run completely, then calculate When The pixel is selected as a temperature anomaly pixel, and the temperature anomaly pixels form temperature anomaly information.

[0043] The above operation is performed on q temperature inversion results one by one to obtain q temperature anomaly information.

[0044] Further, the time frequency threshold is determined, and the time series temperature anomaly information is extracted, and the specific steps are as follows:

[0045] Determine the time frequency threshold: perform time frequency statistics on the temperature anomaly information in the q Landsat 8 images pixel by pixel to generate an image corresponding to the Landsat 8 image only with geographic coordinate information, assign each geographic coordinate to its corresponding time frequency value (maximum q) to obtain a time frequency graph composed of time frequency values; calculate the time frequency mean mean t and standard deviation std t of each pixel in the time frequency graph according to the statistical method, and select the time frequency threshold TF = mean t + std t ;

[0046] Extract the time series temperature anomaly information: compare the marked time frequency of all temperature anomaly information pixel by pixel with the time frequency threshold TF, consider the pixel position with a marked time frequency greater than the time frequency threshold as belonging to the coal fire area, and the pixel with a time frequency lower than TF does not meet the coal fire burning time continuity rule, and is considered not to belong to the coal fire area.

[0047] Advantages

[0048] The present application provides a time-series temperature anomaly remote sensing extraction method for coal fire area identification, which is convenient to implement and has remarkable effect, mainly embodied in the following three points:

[0049] 1) The single threshold algorithm for extracting temperature anomaly information may cause the loss or redundancy of temperature anomaly area and has subjectivity. Compared with the single threshold algorithm, the present method uses a multi-threshold algorithm based on unbiased estimation to improve the drawbacks of the single threshold algorithm.

[0050] 2) The present method uses an adaptive window and a multi-threshold segmentation standard in the window to extract temperature anomaly information, improves the sliding window algorithm based on unbiased estimation using a sliding window and a single threshold, and improves the stability of the algorithm.

[0051] 3) For the extraction of temperature anomaly information of coal fire area, the present method fully considers the time continuity of coal fire combustion, introduces a time parameter related to coal fire combustion, determines a time frequency threshold of coal fire combustion, and extracts time frequency information through the time frequency threshold to obtain time-series temperature anomaly information. BRIEF DESCRIPTION OF DRAWINGS

[0052] Figure 1 The present application is a flowchart of the time-series temperature anomaly remote sensing extraction method for coal fire area identification.

[0053] Figure 2 The present application is a schematic diagram of temperature inversion results in different periods of the study area, in which the temperature size increases from dark color to light color.

[0054] Figure 3 The present application is a schematic diagram of temperature anomaly extraction results of the embodiment; the dates selected are 20180324, 20180605 and 20180916, the deepest color in the figure is the extracted temperature anomaly area, and the black triangle is the coal fire area point.

[0055] Figure 4 The present application is a time frequency threshold determination diagram of the embodiment; Figure 4 (a) is the time frequency of the extracted 29 scenes of temperature thermal anomaly data, and the time frequency increases from dark color to light color; Figure 4 (b) is the number of pixels corresponding to each time frequency, the average value, the standard deviation and the time frequency threshold are calculated; Figure 4 (c) is the time-series temperature anomaly diagram after time frequency threshold screening.

[0056] Figure 5 The present application is a schematic diagram of time-series temperature anomaly extraction results of the embodiment; the dates selected are 20180324, 20180605 and 20180916, the deepest color is the extracted temperature anomaly area, and the black triangle is the coal fire area point. DETAILED DESCRIPTION

[0057] In order to make the technical problems solved by the present application, technical solutions and beneficial effects clearer, the present application will be further described in detail below in combination with the drawings and examples. It should be understood that the specific examples described herein are only used to explain the present application and not used to limit the present application.

[0058] 1. Study area and experimental data profile

[0059] This experiment selects a certain place as the study area; 29 scenes of Landsat 8 satellite data with cloud coverage less than 10% from October 2017 to December 2019 are selected, and the spatial resolutions of the visible light (bands 1-7 and 9) and thermal infrared (bands 10-11) (TIRS) are 30 m and 100 m, respectively.

[0060] 2. Temperature inversion of the study area using single-window algorithm Figure 1 The specific steps are as follows:

[0061] Radiometric calibration and atmospheric correction: Import the _MTL file of 29 scenes of Landsat 8 data into ENVI 5.3 in turn; first, perform radiometric correction, select the Radiometric Calibration tool, and perform radiometric correction on the multi-spectral bands and the thermal infrared band in turn, and select band 10 for the thermal infrared band; select the FLAASH Atmospheric Correction tool to perform atmospheric correction on the multi-spectral bands, and select Atmospheric Model parameters according to different months, and select Rural for Aerosol Model parameters, to obtain reflectivity images of visible light bands. Obtain the atmospheric transmittance τ of 29 scenes of images through the atmospheric profile parameter calculation website provided by NASA.

[0062] Calculate the surface emissivity: Select the NDVI tool to calculate the NDVI, select the multi-spectral data after atmospheric correction, and the band will be automatically selected for the infrared band and the near-infrared band. Correct the NDVI value that exceeds the range of [-1, 1] due to atmospheric correction, input (b1lt-1)*(-1)+(b1gt 1)*1+(b1ge-1and b1le 1)*b1 in the bandmath tool, b1 is the NDVI data. Input the calculation formula in the bandmath tool Calculate the fractional vegetation cover (FVC) where NDVIS, NDVIV represent the minimum and maximum values of NDVI within a certain confidence interval (5% and 95% confidence interval). Input the land surface emissivity calculation formula in ENVI software bandmath: (b1 le NDVIS)*0.995+(b1 gt NDVIS and b1 lt NDVIV)*(0.9589+0.086*b2-0.0671*b2*b2)(b1 ge NDVIV)*(0.9625+0.0614*b2-0.0461*b2*b2) where b1 is the NDVI data, b2 is the fractional vegetation cover FVC, to obtain the land surface emissivity LSE.

[0063] Invert the land surface temperature: calculate the intermediate variables C, D, C = LSE * τ, D = (1 - t) * (1 + (1 - LSE) * τ). LSE is the land surface emissivity, τ is the atmospheric transmittance; calculate the on-board brightness temperature T6, input the calculation formula in ENVI software bandmath b1 selects the radiometrically corrected thermal infrared band; invert the land surface temperature of the study area, input the formula in bandmath (-67.355351 * (1 - b1 - b2) + (0.458606 * (1 - b1 - b2) + b1 + b2) * b3 - b2 * T a ) / b1. Where b1, b2, b3 represent the intermediate parameters C, D, T6, T a is the study area's Kelvin temperature on the day (obtained from the atmospheric profile parameters), and the calculation result is the land surface temperature (shown in Figure 2 ) inverted by the single-window algorithm.

[0064] 3. Extract the time-series temperature anomaly information using the time-series temperature anomaly remote sensing extraction method for coal fire area identification, the specific steps are as follows:

[0065] Extract the temperature anomaly information:

[0066] Calculate the mean temperature m of the entire image.

[0067] Select the starting window (4x4), calculate the mean value m0 of the pixels in the window, the threshold values t1 (one standard deviation of the window mean) and t2 (two standard deviations of the window mean), count the number of pixels sampled by the window N (each time a pixel is sampled, N is incremented by 1), and determine the relationship between the window mean m0 and the image mean m. If m0≤m, the window will not be expanded, and Step 3( Figure 1 ) will be run, if m0>m, Step 2_2( Figure 1 ) will be run, and the window will spread outwards by a step of 1, calculating the mean value m n, threshold t1 and t2, the number of times N that the pixel is window-sampled, running Step3 Figure 1 ).

[0068] The number of times M that the pixel is selected as a thermal anomaly pixel is counted Figure 1 , Step3). When the pixel value γ value ≤t1, the pixel is not a thermal anomaly pixel, and M is incremented by 0 based on the original value; when the pixel value t1<γ value ≤t2, the pixel is selected as a low-frequency thermal anomaly pixel, and M is incremented by 0.5 based on the original value; when the pixel value γ value >t2, the pixel is selected as a high-frequency thermal anomaly pixel, and M is incremented by 1 based on the original value.

[0069] It is determined whether m n is 0. If m n =0, the window is not expanded, and it is determined whether the entire image is run to completion. If the entire image is not run to completion, the window is slid, and Step2_1 is run. If m n is not 0, the window is expanded, and it is determined whether m n and m0 are related. If m n >m0, Step2_2 is run. If m n ≤m0, it is determined whether the entire image is run to completion. If the entire image is not run to completion, the window is slid, and Step2_1 is run. Until the entire image is run to completion, the pixel-by-pixel calculation of is performed. When the pixel value γ is selected as a thermal anomaly pixel. The final thermal anomaly image is composed of thermal anomaly pixels. The above operations are repeated for 29 images to obtain temperature anomaly information of the 29 images. The dates are selected as 20180324, 20180605, and 20180916. The deepest color in the figure is the extracted temperature anomaly region, and the black triangle is the coal fire region point Figure 3 .

[0070] Extracting time-series temperature anomaly information: time-frequency statistics are performed on the 29 images of temperature anomaly information pixel by pixel to obtain a time-frequency graph. The time-frequency mean mean t and the standard deviation std t of the time-frequency graph are calculated, and the time-frequency threshold (TF) is TF=mean t +std t . Figure 4 The pixel position where the time frequency is greater than the time frequency threshold is marked by the time frequency threshold TF, wherein Figure 4 (a) is the time-frequency of the extracted 29 images of temperature thermal anomaly data, and the time frequency increases from dark to light; Figure 4(b) Calculate the time-frequency average, standard deviation and time-frequency threshold value for counting the number of pixels corresponding to each time-frequency; Figure 4 (c) Time-series temperature anomaly map screened by time-frequency threshold value. The pixels in the marked pixel positions are retained, and the pixels in the unmarked pixel positions are removed. The time-series temperature pixels constitute a time-series temperature image, the date is selected as 20180324, 20180605 and 20180916, the color of the deepest temperature anomaly region is extracted, and the black triangle is the coal fire region point. Figure 5 As shown).

Claims

1. A time-series temperature anomaly remote sensing extraction method for coal fire area identification, characterized in that, The steps are as follows: Select as many q Landsat 8 image data as possible in the research time range with cloud cover less than 10%; extract the Landsat 8 image data of the measured area with the thermal infrared band, perform radiation calibration and atmospheric correction on the original optical image data, calculate the normalized vegetation index NDVI, vegetation coverage VFC and surface emissivity ε, and obtain the atmospheric profile parameters; Use the single-window algorithm to perform temperature inversion on the obtained Landsat 8 image data to obtain temperature inversion information; Use adaptive windows and multiple thresholds to extract abnormal temperature information from the temperature inversion information; Considering the time continuity of coal fire combustion, introduce the time parameter related to coal fire combustion, determine the time frequency threshold of coal fire pixels, and extract the time series temperature abnormal information according to the time frequency threshold to determine the coal fire area; The steps of using adaptive windows and multiple thresholds to extract abnormal temperature information from Landsat 8 image temperature inversion information are as follows: Initial window selection: the initial adaptive window size should meet the minimum resolution of remote sensing image, the thermal infrared band resolution of Landsat 8 image data is 100m, and the pixel size is 30m*30m, so the initial window size is selected as 4*4; Adaptive window expansion criterion: extract Landsat 8 image and calculate the mean value m of the whole image; The adaptive window starts from the top left corner of the image and calculates the initial window mean value m0, If m0≤m, the window will not be expanded; If m0>m, the window expands by 1 step in all directions, if it meets the image boundary, it expands by 1 step to the non-boundary of the image, and the mean value m of the new expanded window is calculated n (n = 1, 2, 3,...), If m n ≤ m0, the window is no longer expanded; If m n > m0, the window is expanded by 1 step, until m n ≤ m0, the window is not expanded anymore; Adaptive windowed multi-threshold segmentation: compute mean m n (n = 0, 1, 2,...) and standard deviation std within the window, with m n + std as threshold t1, with m n + 2 * std as threshold t2; When the pixel value value in the window satisfies value > m n + 2 * std, the pixel is selected as a hot anomaly pixel, and the intermediate variable M of the total number of times that a certain pixel is selected as a hot anomaly pixel is increased by 1 based on the original value. When the pixel value value in the window satisfies m n + std< value≤m n + 2 * std, the pixel is selected as a low-frequency thermal anomaly pixel, and the intermediate variable M of the total number of times that a certain pixel is selected as a thermal anomaly pixel is increased by 0.5 on the basis of the original value. When the pixel value value in the window satisfies value≤m n +std, the pixel is not considered as a thermal anomaly pixel, and the intermediate variable M of the total number of times that a certain pixel is selected as a thermal anomaly pixel is increased by 0 on the basis of the original value. The method of extracting temperature abnormal information is as follows: Extract temperature abnormal pixels from Landsat 8 image to reflect the extraction of temperature abnormal information: each pixel in Landsat 8 image has intermediate variables N and M, the total number of times that the pixel is selected into the adaptive window N and the total number of times that the pixel is selected as a thermal abnormal pixel M, initial N=0, M=0, The temperature abnormal pixel is extracted by means of two intermediate variables, if then the pixel is marked as a thermal abnormal pixel; The calculation of intermediate variable N is based on the adaptive window expansion criterion, and the selection of pixel by the initial window or the expanded window will increase 1 based on the original intermediate variable N; The calculation of intermediate variable M is based on the multi-threshold segmentation criterion in the adaptive window; From the upper left of Landsat 8 image, the intermediate variables N and M are calculated, after the calculation, the sliding is carried out with a step of 4, the intermediate variables N and M are calculated again, and the process is repeated until the whole image is calculated, and then the pixel-by-pixel calculation is carried out When The pixel is selected as a temperature anomaly pixel, and the temperature anomaly pixels constitute temperature anomaly information; Perform the above operations on the temperature inversion results of q Landsat 8 images one by one to obtain the temperature abnormal information of q Landsat 8 images.

2. The method of claim 1, wherein the method comprises: Import the image into ENVI software to correct the multispectral band and thermal infrared band in turn, perform atmospheric correction on the multispectral band to obtain the reflectivity image of the visible band, and obtain the atmospheric transmittance τ of the image on the research day through the NASA atmospheric profile parameter calculation website.

3. The method of claim 1, wherein the method comprises: The temperature inversion process of the single-window algorithm is as follows: Calculate the surface emissivity of q Landsat 8 image data selected in turn: Calculate NDVI: input the calculation formula of normalized vegetation index in ENVI software bandmath: float(b4-b3)) / (b4+b3); where b3 and b4 are the red band and near-infrared band after atmospheric correction, respectively. Due to atmospheric correction, the value of NDVI will exceed [-1, 1], so use the method of removing outliers to remove the corrected NDVI; The vegetation cover FVC is calculated using the formula FVC = NDVIV - NDVIS where NDVIS, NDVIV represent the minimum and maximum values of NDVI at a certain confidence interval, typically 5% and 95%. The surface emissivity LSE is calculated, and the emissivity of different ground objects is different. Generally, remote sensing ground objects are divided into water, town and natural surface three kinds; the emissivity calculation formula of three different types of ground object pixels is: E(water) = 0.995, E(building) = 0.9589 + 0.086*FVC-0.0671*FVC*FVC, E(building) = 0.9589 + 0.086*FVC-0.0671*FVC*FVC; The surface emissivity calculation formula is input in the bandmath of ENVI software: (b1 le NDVIS)*0.995+(b1 gt NDVIS and b1 lt NDVIV)*(0.9589+0.086*b2-0.0671*b2*b2)(b1 ge NDVIV)*(0.9625+0.0614*b2-0.0461*b2*b2), Wherein, b1 is NDVI, b2 is vegetation coverage FVC, and the surface emissivity LSE is obtained; The surface temperature is retrieved: first, the intermediate variables C and D are calculated by the formula C = LSE*tau, D = (1-t)*(1+(1-LSE)*tau); LSE is the surface emissivity, and tau is the atmospheric transmittance; The on-board brightness temperature T6 is calculated, and the calculation formula is input in the bandmath of ENVI software: b1 selects the thermal infrared band after radiation correction, and the result is T6; Inversion of land surface temperature in the study area, input formula in bandmath (-67.355351*(1-b1-b2)+(0.458606*(1-b1-b2)+b1+b2)*b3-b2*T a ) / b1; In the formula, b1, b2, b3 respectively represent intermediate parameters C, D, T6, T a In order to study the Kelvin temperature of the area on the day, the calculation result obtained through the atmospheric profile parameter is the surface temperature inverted by the mono-window algorithm.

4. The method of claim 1, wherein the method comprises: The time frequency threshold is determined, and the time sequence temperature anomaly information is extracted, and the specific steps are as follows: Determine the time frequency threshold: the temperature anomaly information in q amplitude Landsat 8 image is carried out time frequency statistics pixel by pixel, generates an image corresponding to Landsat 8 image only with geographic coordinate information, assigns each geographic coordinate to its corresponding time frequency value, the maximum is the total number of images q, obtains the time frequency graph composed of time frequency value; according to the statistical method, the time frequency mean value mean t And standard deviation std t Of each pixel in the time frequency graph are calculated, and the time frequency threshold TF = mean t +std t Is selected. Extracting time sequence temperature anomaly information: compare all the time frequency of pixel marked temperature anomaly information with the time frequency threshold TF one by one, and consider that the pixel position with the marked time frequency greater than the time frequency threshold belongs to the coal fire area, and the pixel with the time frequency lower than TF does not conform to the time continuity rule of coal fire combustion, and is considered not to belong to the coal fire area.

Images