Surface temperature anomaly detection and attribution analysis method based on thermal infrared data

Abstract

The invention relates to the technical field of geothermal detection, in particular to a surface temperature anomaly detection and attribution analysis method based on thermal infrared data. Comprising the following steps: S1, collecting MODIS thermal infrared brightness temperature data, and forming a sequence input sample of each city; s2, constructing a prediction model, and training the prediction model to obtain a reconstruction error and a regression residual error of the prediction model; s3, obtaining a joint anomaly score of each month of each city, setting a judgment threshold value, when the joint anomaly score is greater than the judgment threshold value, indicating that the corresponding month is an abnormal month, and otherwise, indicating that the corresponding month is a normal month; s4, for each abnormal month of each city, acquiring a plurality of core variables, performing standardization processing and multiple linear regression on the plurality of core variables, and quantifying the contribution degree of each core variable to temperature anomaly; the accuracy of surface temperature prediction and the accuracy of surface temperature anomaly screening are improved, and quantitative analysis is carried out on the cause of abnormal temperature.

Description

Technical Field

[0001] This invention relates to the field of geothermal detection technology, and in particular to a method for detecting and attributing anomalies in land surface temperature based on thermal infrared data. Background Technology

[0002] Surface temperature, as a crucial physical parameter reflecting the energy exchange at the Earth's surface, directly characterizes urban thermal environment features and is a key indicator for assessing the urban heat island effect and analyzing spatial differentiation. Related research indicates that surface temperature information obtained through remote sensing has significant application value in urban-scale thermal environment monitoring, anomaly area identification, and environmental response analysis. Particularly at regional and long-term scales, constructing stable and continuous surface temperature sequences based on thermal infrared remote sensing data is of great importance for revealing the evolution patterns of the urban thermal environment.

[0003] Currently, land surface temperature acquisition primarily relies on multi-source Earth observation satellite data. Among them, the Moderate Resolution Imaging Spectroradiometer (MODIS), by setting multiple thermal infrared observation bands, provides fundamental data support for land surface temperature retrieval and atmospheric parameter estimation. Its resulting land surface temperature products exhibit good continuity and consistency in both temporal scale and spatial coverage, making them suitable for long-term thermal environment monitoring research at regional and even global scales. Besides MODIS, data from the Landsat series and Sentinel satellites are also widely used in land surface environment research. However, differences in spatial resolution, revisit periods, and sensor design among different satellite systems result in varying applicability for long-term series analysis.

[0004] In terms of surface temperature analysis methods, existing technologies mainly focus on surface temperature inversion and anomaly identification driven by physical models or statistical methods. These methods typically rely on radiative transfer models, single-window algorithms, or empirical formulas, and are highly dependent on intermediate parameters such as surface emissivity and atmospheric water vapor content. Parameter acquisition is complex, and errors are prone to accumulate at multiple stages. Furthermore, with the continuous growth in the time span and volume of remote sensing observation data, traditional methods struggle to balance computational efficiency with the ability to represent temporal structure when processing long-term series data. In recent years, machine learning and deep learning methods have been gradually introduced into the field of remote sensing time series analysis to uncover temporal patterns and potential changes in the data. However, some existing models still have limitations in long-series modeling capabilities, parameter stability, or computational complexity, making it difficult to simultaneously ensure the representation of temporal continuity and meet the needs of large-scale applications.

[0005] Relevant patent documents retrieved:

[0006] This document, published in China (CN120654148A) on September 16, 2025, discloses a geothermal anomaly identification system based on remote sensing. The system comprises a data layer, a processing layer, an analysis layer, and an output layer. The data layer collects data including thermal infrared remote sensing data, multispectral / hyperspectral remote sensing data, radar data, geological and topographic data, meteorological data, and geophysical data. The processing layer extracts key features, captures dynamic changes in geothermal anomalies, distinguishes between persistent geothermal anomalies and transient events, and combines hyperspectral data with magnetotelluric sounding data to construct a "surface..." The deep layer is used to learn the diffusion law of geothermal anomalies along the fault zone through the graph propagation model. The physical information neural network is used to constrain the heat conduction equation of the geothermal field as the loss function, and the temperature field prediction model is jointly trained. Based on gravity and geophysical data, a deep three-dimensional geological structure is constructed. The geothermal anomaly detection model is constructed through the features provided by the processing layer to identify geothermal anomaly points. The output layer is used to generate a geothermal anomaly probability map.

[0007] The prior art represented by the aforementioned literature has at least the following unresolved technical problems or defects: The prior art represented by CN120654148A can identify the spatial distribution and diffusion patterns of geothermal anomalies along fault zones to a certain extent, and improve the accuracy of anomaly identification by using the heat conduction equation as a constraint condition through a physical information neural network. However, this type of method mainly targets short-term or transient changes in geothermal anomalies, lacking a systematic approach to the long-term temporal evolution characteristics of surface temperature. Furthermore, the system does not fully consider the quantitative correlation between surface temperature anomalies and multiple environmental factors, making it difficult to achieve in-depth attribution analysis of the anomaly formation mechanism. In addition, due to its reliance on complex deep structural data and physical constraints, the system faces challenges in data consistency and computational efficiency when processing large-scale, multi-year surface temperature time-series data, making it difficult to meet the needs of continuous monitoring and comprehensive analysis of surface temperature anomalies.

[0008] How to reliably extract anomaly information while ensuring stable modeling of land surface temperature time series data, and further reveal the influence of environmental factors on land surface temperature anomalies, remains a problem that needs further research in the current technology.

[0009] Therefore, there is an urgent need to provide a method for detecting and attributing surface temperature anomalies based on thermal infrared data, which can improve the accuracy of surface temperature prediction, enhance the accuracy of surface temperature anomaly screening, and quantitatively analyze the causes of abnormal temperatures compared to existing technologies. Summary of the Invention

[0010] This invention addresses the technical problems existing in the prior art and provides a method for detecting and attributing anomalies in land surface temperature based on thermal infrared data.

[0011] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A method for detecting and attributing anomalies in land surface temperature based on thermal infrared data includes the following steps: S1. Collect MODIS thermal infrared brightness temperature data, divide it according to the four bands of MODIS, and form a sequence input sample for each city. S2. Construct a prediction model. The constructed prediction model structure includes an LSTM encoder, a pooling module, an LSTM decoder, and a regression head. Process the sequence input samples to obtain a time series dataset. Use the time series dataset to train the prediction model to obtain the learned prediction model. Obtain the reconstruction error and regression residual of the prediction model. S3. Based on the reconstruction error and regression residual of the prediction model, obtain the joint anomaly score for each city and each month, set a judgment threshold, and compare the joint anomaly score with the judgment threshold. If the joint anomaly score is greater than the judgment threshold, it indicates that the corresponding month is an abnormal month; otherwise, it is a normal month. S4. Filter out the abnormal months and their corresponding cities identified in step S3. For each abnormal month in each city, obtain multiple core variables. Construct a local time window for the abnormal month and the k months preceding it. Then, perform standardization and multiple linear regression on the multiple core variables to quantify the contribution of each core variable to the temperature anomaly.

[0012] Furthermore, S2 specifically includes the following steps: S21. Set the hidden layer dimensions of the LSTM model. For any time step t, the LSTM model first checks the hidden state from the previous time step. The forget gate is computed together with the current input x. Input gate With output gate Forgotten Gate Input gate With output gate All values ​​are obtained by the Sigmoid activation function and are within the range of [0,1]; forget gate Used to determine the memory unit of the previous moment. Dimensions that need to be preserved, input gate Used to identify candidate memories Write to the current memory cell, output gate Used to control the degree to which the memory state after nonlinear transformation is exposed as a hidden state; S22. Update the memory units and hidden states, and output the final memory units and the final hidden states; S23, the final memory unit output by the encoder. and the final hidden state After passing through the pooling module, the data is simultaneously input into both the regression head and the decoder; the regression head then processes the data... , The monthly average temperature prediction value is obtained after processing. At the same time, the true value of the monthly average temperature y is obtained; the decoder will... , As the initial state of the decoder, the zero vector is used as the decoding input. The original 12-month input sequence is reconstructed by generating sequences step by step over time, resulting in the reconstructed sequence. To obtain the true sequence Based on the monthly average temperature forecast , Monthly average temperature true value y, Reconstructed sequence Real sequence This yields the overall optimization objective of the prediction model during the training process; S24. Using the prediction model trained in step S2, obtain the predicted monthly average temperature for the identification phase. Actual monthly average temperature Reconstructed sequence Real sequence The reconstruction error and regression residuals are obtained.

[0013] Furthermore, in step S22, the memory cells are updated using the following formula: ; In the above formula, This represents the memory unit corresponding to time step t.

[0014] Furthermore, in step S22, the hidden state is updated using the following formula: ; In the above formula, This represents the hidden state corresponding to time step t.

[0015] Furthermore, S23 specifically includes the following steps: S231, Based on the reconstructed sequence Real sequence The reconstruction loss is obtained; S232, Based on the predicted monthly average temperature value Given the true monthly average temperature value y, calculate the temperature regression loss for the corresponding month. S233. Based on the reconstruction loss and temperature regression loss, the overall optimization objective of the model during training is obtained, calculated using the following formula: ; In the above formula, This indicates the overall optimization goal. Indicates the reconstruction loss. .

[0016] Furthermore, the reconstruction loss and temperature regression loss are specifically calculated using the following formula: ; ; In the above formula, Indicates the time step.

[0017] Furthermore, S24 specifically includes the following steps: S241, Based on the reconstructed sequence Real sequence The reconstruction error of the LSTM decoder is obtained, which is used to characterize the degree of structural deviation between the sequence and the historical pattern. Specifically, it is calculated using the following formula: ; In the above formula, Indicates reconstruction error; S242. The regression residuals are calculated using the following formula: ; In the above formula, This represents the regression residual.

[0018] Furthermore, S3 includes the following steps: S31. Normalize the reconstruction error and regression residuals using the following formula: ; ; In the above formula, This represents the normalized reconstruction error. express The standard deviation on the training set. This represents the normalized regression residuals. express Standard deviation on the training set; S32. Calculate the joint anomaly score based on the normalized reconstruction error and the normalized regression residual, specifically using the following formula: ; In the above formula, Indicates joint anomaly score, Indicates the first weight. This indicates the second weight, where the sum of the first and second weights equals 1. S33. Set a judgment threshold. The judgment threshold is obtained by processing the average of the joint scores corresponding to the samples that deviate most from the normal mode by 5% in the training set using the following formula: ; In the above formula, Indicates the threshold for judgment. Represents the 95th percentile. This represents the average joint score of the samples in the training set that deviate most from the normal pattern by 5%. S34. Determine whether the average monthly temperature of the city in this month is an abnormal month or a normal month.

[0019] Furthermore, in step S4, the core variables include three core variables: population density, NDVI, and evaporation. For population density, the annual population density of the corresponding city is obtained, and for NDVI and evaporation, the corresponding values ​​for the abnormal month are obtained.

[0020] Furthermore, S4 specifically includes the following steps: S41. For each month with an anomaly, construct an analysis window, represented as follows: , This represents the set of temperature values ​​for the i-th abnormal month; S42. For each abnormal month, in the analysis window Within this framework, we construct a system incorporating multiple linear regression, expressed as: ; In the above formula, Let P represent the standardized temperature and P represent the standardized population density. Indicates standardized evaporation rate. Indicates standardized NDVI, Represents a constant term. , , , , , Both represent regression coefficients. Indicates the error term; S43. Based on multiple linear regression, calculate the contribution of each core variable. The larger the contribution, the greater the contribution of that core variable to the temperature anomaly. The specific contribution is calculated using the following formula: ; In the above formula, This represents the contribution of the j-th core variable, where j ranges from 1 to 3. This represents the regression coefficient of the j-th core variable in the multiple regression model. Let represent the regression coefficient of the l-th variable, where l ranges from 1 to q, and q represents the total number of core variables included in the regression model.

[0021] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention constructs a nationwide monthly land surface temperature dataset based on MODIS thermal infrared brightness temperature data. A novel strategy model using an LSTM-based encoder-decoder framework is proposed, employing a sliding window spanning multiple months for sequence prediction and reconstruction. A semi-supervised anomaly detection mechanism is built by combining reconstruction error and prediction residuals to achieve automated identification of temperature anomalies, improving the accuracy of land surface temperature prediction and anomaly screening. Furthermore, the anomaly detection model can identify temperature abrupt change points in several typical cities across the country, analyzing the significant spatial heterogeneity of the dominant factors influencing temperature anomalies in different cities, and quantitatively analyzing the causes of anomalous temperatures. These research results provide strong technical support for the precise monitoring, scientific regulation, and sustainable environmental and climate governance of urban thermal environments. Attached Figure Description

[0022] Figure 1 This is a flowchart of the method of the present invention.

[0023] Figure 2 This is a structural block diagram of the prediction model of this invention.

[0024] Figure 3 This is a schematic diagram of the reconstruction error distribution of the encoder of the present invention on the evaluation set.

[0025] Figure 4 This is a diagram showing the predicted monthly average temperature for different cities and the anomaly detection results.

[0026] Figure 5 This is a schematic diagram illustrating the attribution analysis of abnormal temperature data in different cities. Detailed Implementation

[0027] The technical solution of the present invention will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are not all embodiments of the present invention. All other embodiments obtained by those skilled in the art without creative effort are within the protection scope of the present invention.

[0028] like Figure 1 As shown, this invention provides a method for detecting and attributing surface temperature anomalies based on thermal infrared data, comprising the following steps: S1. Collect MODIS thermal infrared brightness temperature data nationwide from 2002 to 2022, perform monthly synthesis to form a sequence input sample, initially divide according to the MODIS thermal infrared brightness temperature data of each city, and then divide according to the four bands of MODIS. For each band, sort the MODIS thermal infrared brightness temperature data in chronological order to form a sequence input sample, represented as: ; In the above formula, This represents a sequence of input samples from a city. This represents the dataset for the first band of MODIS. This represents the dataset for the second band of MODIS. This represents the dataset for the third band of MODIS. This represents the dataset for the fourth band of MODIS. The feature dimension is represented by F=4 in this invention.

[0029] MODIS thermal infrared brightness temperature data can be seen in Table 1 below: Table 1

[0030] S2. For each band of the sequence input samples generated in step S1, extract the corresponding monthly average data to form the monthly average brightness temperature sequence for each city. The monthly average brightness temperature sequences for all cities form a time-series dataset. Construct a prediction model, which includes an LSTM encoder, pooling module, LSTM decoder, and regression head. Train the prediction model using the time-series dataset to obtain the learned prediction model, and obtain the reconstruction error and regression residuals. Specifically, the following steps are included: S21. Set the hidden layer dimension H of the LSTM model to 64. For any time step t, the LSTM model first checks the hidden state from the previous time step. The forget gate is computed together with the current input x. Input gate With output gate Forgotten Gate Input gate With output gate All are obtained by processing with the Sigmoid activation function, and their values ​​are all between [0,1]. They are used to control the retention ratio of old memories, the writing ratio of new information, and the output ratio of the current state, respectively. Among them, the forget gate... Used to determine the memory unit of the previous moment. Dimensions that need to be preserved, input gate Used to identify candidate memories Write to the current memory cell, output gate This is used to control the degree to which a memory state after a nonlinear transformation is exposed as a hidden state; where candidate memories... Generated by tanh(.), it represents the potential new information extracted by the model from the current input without considering gating factors.

[0031] S22. Update memory units and hidden states; for memory units, the current long-term memory representation is formed by integrating forgotten old memories with filtered new information; the hidden state represents the effective sequence features exposed at time step t and serves as the input for the next time step; after 12 time steps of recursive updates to memory units and hidden states, the final memory units and the final hidden state are output. The final hidden state is regarded as a high-dimensional semantic representation of the entire input sequence and is used for subsequent temperature prediction and sequence reconstruction tasks.

[0032] The memory unit is updated using the following formula: ; In the above formula, This represents the memory unit corresponding to time step t.

[0033] The hidden state is updated using the following formula: ; In the above formula, This represents the hidden state corresponding to time step t.

[0034] S23, the final memory unit output by the encoder. and the final hidden state After passing through the pooling module, the data is simultaneously input into both the regression head and the decoder; the regression head then processes the data... , The monthly average temperature prediction value is obtained after processing. At the same time, the true value of the monthly average temperature y is obtained; the decoder will... , As the initial state of the decoder, the zero vector is used as the decoding input. The original 12-month input sequence is reconstructed by generating sequences step by step over time, resulting in the reconstructed sequence. To obtain the true sequence Based on the monthly average temperature forecast , Monthly average temperature true value y, Reconstructed sequence Real sequence The overall optimization objective of the model during training is obtained, which includes the following steps: S231. Reconstruct the sequence based on the training. Real sequence The reconstruction loss is obtained, and the prediction model is constrained by the reconstruction loss to extract more robust temporal structure features. The specific reconstruction loss is calculated by the following formula: ; In the above formula, Indicates the reconstruction loss. Indicates the time step. .

[0035] S232, Based on the predicted monthly average temperature value The actual monthly average temperature y is used to calculate the temperature regression loss for the corresponding month, specifically calculated using the following formula: ; In the above formula, This represents the temperature regression loss.

[0036] S233. Based on the reconstruction loss and temperature regression loss, the overall optimization objective of the model during training is obtained, specifically expressed as: ; In the above formula, This indicates the overall optimization goal.

[0037] By jointly minimizing the aforementioned losses, the model reconstructs the original temporal structure as much as possible while maintaining prediction accuracy, thereby obtaining a more stable temporal representation and residual signal with clearer physical meaning for subsequent anomaly detection. The final training objective consists of both reconstruction loss and regression loss, enabling the model to simultaneously learn regression patterns that depend on long-term structural and temperature changes, providing a more reliable representation and residual signal for the subsequent anomaly detection module.

[0038] S24. Using the prediction model trained in step S2, obtain the predicted monthly average temperature for the identification phase. Actual monthly average temperature Reconstructed sequence Real sequence The reconstruction error and regression residuals are obtained.

[0039] S241, Based on the reconstructed sequence Real sequence The reconstruction error of the LSTM decoder is obtained, which is used to characterize the degree of structural deviation between the sequence and the historical pattern. It belongs to unsupervised anomalous signals and is specifically calculated by the following formula: ; In the above formula, This indicates the reconstruction error.

[0040] S242. Due to the abnormal deviation in the relationship between real temperature and remote sensing features, the regression head of the prediction model will exhibit significant prediction bias in regression tasks, resulting in regression residuals. These residuals can reflect abnormal supervised signals at the target level. The regression residuals are calculated using the following formula: ; In the above formula, This represents the regression residual.

[0041] S3. Based on the reconstruction error and regression residual, obtain the joint anomaly score for the temperature of each city for each month. Set a judgment threshold and compare the joint anomaly score with the threshold to determine whether the average monthly temperature of that city in that month is an abnormal month or a normal month. This process identifies months with abnormal temperatures and includes the following steps: S31. Normalize the reconstruction error and regression residuals using the following formula: ; ; In the above formula, This represents the normalized reconstruction error. express The standard deviation on the training set is used to normalize the reconstruction error. This represents the normalized regression residuals. express The standard deviation on the training set is used to normalize the regression residuals; S32. Calculate the joint anomaly score based on the normalized reconstruction error and the normalized regression residual, specifically using the following formula: ; In the above formula, Indicates joint anomaly score, Indicates the first weight. The second weight is represented by the sum of the first and second weights, which is 1. The first and second weights are preferably both 0.5.

[0042] This design allows for the simultaneous capture of structural and data anomalies, thereby improving the robustness and sensitivity of anomaly detection. When certain months exhibit both significant structural and predictive deviations, their joint score will be significantly higher than that of normal areas, enabling the model to identify potential extreme weather or observational problems with high confidence.

[0043] S33. Set a judgment threshold. The judgment threshold is obtained by processing the average of the joint scores corresponding to the samples that deviate most from the normal mode by 5% in the training set using the following formula: ; In the above formula, Indicates the threshold for judgment. Represents the 95th percentile. This represents the average joint score of the samples in the training set that deviate most from the normal pattern by 5%.

[0044] S34. Determine whether the average monthly temperature of a city in a given month is an abnormal month or a normal month. If the combined abnormal score of a city in a given month is greater than the judgment threshold, it indicates that the month is an abnormal month; otherwise, it is a normal month.

[0045] like Figure 3 As shown, the histogram displays the reconstruction error calculated by the prediction model for each sequence segment. The larger the error, the more the pattern of that time segment deviates from the normal temporal structure learned during training.

[0046] S4. Filter out the abnormal months identified in step S3 and their corresponding cities. For each abnormal month in each city, obtain three core variables: population density, NDVI, and evaporation. For population density, obtain the city's annual population density; for NDVI and evaporation, obtain the corresponding values ​​for that abnormal month. A local time window is formed by the abnormal month and the k months preceding it. By standardizing these three core variables and performing multiple linear regression, quantify the contribution and direction of each core variable to temperature anomalies. Specifically, this includes the following steps: S41. For each month with an anomaly, construct an analysis window, represented as follows: ; In the above formula, This represents the set of temperature values ​​for the i-th abnormal month. Indicates the first three months, Indicates the first two months, Indicates the current month.

[0047] S42. For each abnormal month, in the analysis window Within this framework, a multiple linear regression including main effects and interaction terms is constructed, expressed as: ; In the above formula, Let P represent the standardized temperature and P represent the standardized population density. Indicates standardized evaporation rate. Indicates standardized NDVI, Represents a constant term. , , , , , Both represent regression coefficients. This indicates the error term.

[0048] S43. Based on multiple linear regression, calculate the contribution of each core variable. The larger the contribution, the greater the contribution of that core variable to the temperature anomaly. The specific contribution is calculated using the following formula: ; In the above formula, This represents the contribution of the j-th core variable, where j ranges from 1 to 3; This represents the regression coefficient of the j-th core variable (such as population density, NDVI, evaporation) in the multiple regression model, which reflects the direction and intensity of the variable's influence on temperature anomalies. Let represent the regression coefficient of the l-th variable, where l ranges from 1 to q, and q represents the total number of core variables included in the regression model. In this study, q = 3 (i.e., population density, NDVI, and evaporation).

[0049] This invention verifies the stability and generalization performance of the prediction model's temperature prediction capability across different time spans in multi-step prediction tasks: First, quantitative analysis was conducted. Multi-step prediction (Step 3, Step 6, and Step 12) can reflect the model's ability to capture local disturbances, seasonal structures, and long-term trends at different scales. Therefore, comparing the prediction error indices under different step lengths helps to clarify the optimal configuration of the model in actual climate prediction tasks. Table 2 summarizes the MAE, RMSE, and R2 results for the three prediction step lengths to reveal the impact of step length changes on model performance. As shown in Table 2, although all step lengths achieved high accuracy (R2 all exceeded 0.95), Step 12 performed best in terms of comprehensive indices, indicating that the model has stronger stability and generalization ability in climate series prediction over long periods.

[0050] Table 2

[0051] Specifically, Step 12 achieved a MAE of 1.1503, an RMSE of 1.6638, and an R² of 0.9765, representing the lowest error and highest fit among the three settings. This performance indicates that the model can effectively extract long-term dependency information from the 12-month input sequence, enabling temperature predictions to have both low bias and accurately reproduce the overall temperature change trend of the real sequence. In contrast, Step 6 showed a slight decrease in prediction accuracy, with MAE and RMSE of 1.3855 and 2.0653, respectively, and R² dropping to 0.9634. Step 3 further declined to MAE 1.6057, RMSE 2.3774, and R² 0.9512, indicating that while shorter time windows focus more on local perturbation features, their ability to capture the complete seasonal structure is relatively insufficient.

[0052] Overall, while all three step lengths achieved high levels of fit, Step 12 significantly outperformed Step 6 and Step 3, maximizing the utilization of long-term structural features of the time series without increasing model complexity. Therefore, choosing 12 months as the window length is more suitable for the urban temperature evolution modeling task in this study, not only maintaining lower prediction errors but also providing a more stable and reliable temporal representation for subsequent anomaly detection.

[0053] To further verify the model's applicability and anomaly detection capabilities in different climate regions, Figure 4 The actual and predicted temperature sequences of eight typical cities were visualized, and the months of anomalies were marked. The broken line with triangles represents the actual temperature, the broken line with circles represents the model predicted temperature, the gray shaded area is the abnormal time window identified by the model, and the black circle marks the abnormal months where the joint anomaly score exceeds the threshold. Figure 4 In the figures, (a) Beijing; (b) Chengdu; (c) Hong Kong Special Administrative Region; (d) Shanghai; (e) Shigatse; (f) Qujing; (g) Jinzhong; and (h) Shenzhen. As shown in the figure, all detected anomalies correspond to significant temperature jumps or trend shifts in the actual temperature series, indicating that the model can make a stable response to structural deviations in actual climate change.

[0054] Beijing ( Figure 4 Taking (a) as an example, the observed temperature plummeted at the end of 2014 and then rebounded sharply by about 10°C at the beginning of 2015, forming a short-term surge. The magnitude of the change deviated significantly from the normal seasonal curve. The model showed a significant prediction bias at this point, and the joint anomaly score exceeded the threshold, thus accurately marking the anomaly point. Chengdu ( Figure 4The sequence in (b) is generally smooth, but between June and July 2015, observations showed a local peak increase of approximately 3–4°C. This peak was higher than the preceding and following months, causing the prediction curve to fail to fully follow this sudden increase, thus raising the anomaly score. The gray shaded area in the figure precisely covers this short-term disturbance. Hong Kong Special Administrative Region ( Figure 4 In the temperature series (c) of the model, the observed values ​​fluctuate significantly. A sharp temperature drop of approximately 8–10°C occurred in the winter of 2015. During this period, the predicted values ​​were significantly higher than the actual values, resulting in a large regression bias. The model outputs continuous anomalous windows for the month in which this extreme temperature drop occurred and its adjacent periods. This result indicates that the model has a clear response capability to such drastic seasonal disturbances. For Shanghai (… Figure 4 In (d) of the study, a slight but persistent temperature rise was observed in the spring of 2015, with values ​​approximately 2–3°C higher than the normal trend. Although the absolute magnitude was not extreme, its persistent nature caused continuous bias in the LSTM predictions over multiple steps, resulting in outlier scores crossing thresholds and forming independent outliers. Shigatse City, located in the plateau region ( Figure 4 (e) shows a relatively larger temperature fluctuation. In the winter of 2015, the actual temperature rapidly dropped to near a local minimum within a short period, forming a sharp structural abrupt change with a drop exceeding 10°C. Model prediction lag caused both reconstruction error and regression residuals to increase simultaneously, thus enabling the stable detection of this anomaly. In summary, the model constructed in this study demonstrates robust regression performance. It can effectively identify temperature anomalies in temperature sequences from different geographical regions and cities at different levels of development across the country, laying a reliable data analysis foundation for further in-depth research into the meteorological and socioeconomic causes behind these anomaly detection points.

[0055] like Figure 5The study selected eight typical cities (Qujing, Jinchang, Shenzhen, Chengdu, Hong Kong SAR, Beijing, Shanghai, and Shigatse) and systematically compared the main effects, interaction effects, and policy simulation results of three factors: population, evaporation, and normalized vegetation index (NDVI). From the perspective of city level and geographical location, population density remains the primary cause of warming in most cities (such as Qujing, Jinchang, and Shanghai), contributing 14% to 39%. The direction of the effects of vegetation and evaporation shows significant regional differences. For example, in Shenzhen and Chengdu, NDVI exhibits a cooling effect; while in Beijing and Shanghai, vegetation cover shows a warming effect. This difference may be related to vegetation type, climate background, and urbanization stage. Southern cities (such as Shenzhen) have predominantly evergreen broad-leaved forests, resulting in a strong transpiration-cooling effect; while in Beijing and Shanghai, greening is mostly artificial turf and buffer zones. In arid and semi-arid climates or under high summer temperatures, water stress may weaken transpiration, or even cause a warming effect due to albedo being lower than the urban underlying surface. Regarding the differences, the impact of the interaction effect is particularly prominent: the synergistic effect of population and evaporation manifests as enhanced warming in Qujing and Chengdu, while exhibiting an inhibitory effect in Jinchang and Shenzhen. This vividly reflects the "policy transplantation" dilemma in comparative public policy research. A noteworthy anomaly is that in high-density cities such as Hong Kong and Beijing, the main effect of population density on temperature is weak, or even negative. This phenomenon may stem from the shadow effect created by high-rise buildings, wind field changes caused by complex structures, or structural differences in evaporation conditions within high-density areas.

[0056] This invention constructs a nationwide monthly land surface temperature dataset based on MODIS thermal infrared brightness temperature data from 2002 to 2022. A novel anomaly detection model using an LSTM-based encoder-decoder framework is proposed, employing a 12-month sliding window for sequence prediction and reconstruction. A semi-supervised anomaly detection mechanism is built by combining reconstruction error and prediction residuals to achieve automated identification of temperature anomalies. Furthermore, the anomaly detection model can identify temperature abrupt change points in multiple typical cities across the country, and the analysis reveals significant spatial heterogeneity in the dominant factors of temperature anomalies in different cities. These findings provide strong technical support for the accurate monitoring, scientific regulation, and sustainable environmental and climate governance of urban thermal environments.

[0057] Finally, it should be noted that the above content is only used to illustrate the technical solution of the present invention, and is not intended to limit the scope of protection of the present invention. Simple modifications or equivalent substitutions made by those skilled in the art to the technical solution of the present invention do not depart from the essence and scope of the technical solution of the present invention.

Claims

1. A method for detecting and attributing surface temperature anomalies based on thermal infrared data, characterized in that, Includes the following steps: S1. Collect MODIS thermal infrared brightness temperature data, divide it according to the four bands of MODIS, and form a sequence input sample for each city. S2. Construct a prediction model. The constructed prediction model structure includes an LSTM encoder, a pooling module, an LSTM decoder, and a regression head. Process the sequence input samples to obtain a time series dataset. Use the time series dataset to train the prediction model to obtain the learned prediction model. Obtain the reconstruction error and regression residual of the prediction model. S3. Based on the reconstruction error and regression residual of the prediction model, obtain the joint anomaly score for each city and each month, set a judgment threshold, and compare the joint anomaly score with the judgment threshold. If the joint anomaly score is greater than the judgment threshold, it indicates that the corresponding month is an abnormal month; otherwise, it is a normal month. S4. Filter out the abnormal months and their corresponding cities identified in step S3. For each abnormal month in each city, obtain multiple core variables. Construct a local time window for the abnormal month and the k months preceding it. Then, perform standardization and multiple linear regression on the multiple core variables to quantify the contribution of each core variable to the temperature anomaly.

2. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 1, characterized in that, S2 specifically includes the following steps: S21. Set the hidden layer dimensions of the LSTM model. For any time step t, the LSTM model first checks the hidden state from the previous time step. The forget gate is computed together with the current input x. Input gate With output gate Forgotten Gate Input gate With output gate All values ​​are obtained by the Sigmoid activation function and are within the range of [0,1]; forget gate Used to determine the memory unit of the previous moment. Dimensions that need to be preserved, input gate Used to identify candidate memories Write to the current memory cell, output gate Used to control the degree to which the memory state after nonlinear transformation is exposed as a hidden state; S22. Update the memory units and hidden states, and output the final memory units and the final hidden states; S23, the final memory unit output by the encoder. and the final hidden state After passing through the pooling module, the data is simultaneously input into both the regression head and the decoder; the regression head then processes the data... , The monthly average temperature prediction value is obtained after processing. At the same time, the true value of the monthly average temperature y is obtained; the decoder will... , As the initial state of the decoder, the zero vector is used as the decoding input. The original 12-month input sequence is reconstructed by generating sequences step by step over time, resulting in the reconstructed sequence. To obtain the true sequence Based on the monthly average temperature forecast , Monthly average temperature true value y, Reconstructed sequence Real sequence This yields the overall optimization objective of the prediction model during the training process; S24. Using the prediction model trained in step S2, obtain the predicted monthly average temperature for the identification phase. Actual monthly average temperature Reconstructed sequence Real sequence The reconstruction error and regression residuals are obtained.

3. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 2, characterized in that, In step S22, the memory cells are updated using the following formula: ； In the above formula, This represents the memory unit corresponding to time step t.

4. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 3, characterized in that, In step S22, the hidden state is updated using the following formula: ； In the above formula, This represents the hidden state corresponding to time step t.

5. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 2, characterized in that, S23 specifically includes the following steps: S231, Based on the reconstructed sequence Real sequence The reconstruction loss is obtained; S232, Based on the predicted monthly average temperature value Given the true monthly average temperature value y, calculate the temperature regression loss for the corresponding month. S233. Based on the reconstruction loss and temperature regression loss, the overall optimization objective of the model during training is obtained, calculated using the following formula: ； In the above formula, This indicates the overall optimization goal. Indicates the reconstruction loss. .

6. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 5, characterized in that, Reconstruction loss and temperature regression loss are specifically calculated using the following formulas: ； ； In the above formula, Indicates the time step.

7. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 2, characterized in that, S24 specifically includes the following steps: S241, Based on the reconstructed sequence Real sequence The reconstruction error of the LSTM decoder is obtained, which is used to characterize the degree of structural deviation between the sequence and the historical pattern. Specifically, it is calculated using the following formula: ； In the above formula, Indicates reconstruction error; S242. The regression residuals are calculated using the following formula: ； In the above formula, This represents the regression residual.

8. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 6, characterized in that, S3 includes the following steps: S31. Normalize the reconstruction error and regression residuals using the following formula: ； ； In the above formula, This represents the normalized reconstruction error. express The standard deviation on the training set. This represents the normalized regression residuals. express Standard deviation on the training set; S32. Calculate the joint anomaly score based on the normalized reconstruction error and the normalized regression residual, specifically using the following formula: ； In the above formula, Indicates joint anomaly score, Indicates the first weight. This indicates the second weight, where the sum of the first and second weights equals 1. S33. Set a judgment threshold. The judgment threshold is obtained by processing the average of the joint scores corresponding to the samples that deviate most from the normal mode by 5% in the training set using the following formula: ； In the above formula, Indicates the threshold for judgment. Represents the 95th percentile. This represents the average joint score of the samples in the training set that deviate most from the normal pattern by 5%. S34. Determine whether the average monthly temperature of the city in this month is an abnormal month or a normal month.

9. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 1, characterized in that, In step S4, the core variables include three core variables: population density, NDVI, and evaporation. For population density, the annual population density of the corresponding city is obtained, and for NDVI and evaporation, the corresponding values ​​of the abnormal month are obtained.

10. The method for detecting and attributing surface temperature anomalies based on thermal infrared data according to claim 9, characterized in that, S4 specifically includes the following steps: S41. For each month with an anomaly, construct an analysis window, represented as follows: , This represents the set of temperature values ​​for the i-th abnormal month; S42. For each abnormal month, in the analysis window Within this framework, we construct a system incorporating multiple linear regression, expressed as: ； In the above formula, Let P represent the standardized temperature and P represent the standardized population density. Indicates standardized evaporation rate. Indicates standardized NDVI, Represents a constant term. , , , , , Both represent regression coefficients. Indicates the error term; S43. Based on multiple linear regression, calculate the contribution of each core variable. The larger the contribution, the greater the contribution of that core variable to the temperature anomaly. The specific contribution is calculated using the following formula: ； In the above formula, This represents the contribution of the j-th core variable, where j ranges from 1 to 3. This represents the regression coefficient of the j-th core variable in the multiple regression model. Let represent the regression coefficient of the l-th variable, where l ranges from 1 to q, and q represents the total number of core variables included in the regression model.

Images