Land coverage change and thermal environment influence research method based on multi-temporal data

Abstract

The invention discloses a land cover change and thermal environment influence research method based on multi-temporal data, and the method comprises the following implementation steps: 1), obtaining a multi-temporal remote sensing image of a to-be-monitored region, carrying out the land cover classification mapping research of the remote sensing data of different time phases of the to-be-monitored region, and carrying out the statistics of the area of each land cover type; 2) performing surface temperature inversion on the remote sensing data of different time phases of the to-be-monitored area, grading the temperature, and analyzing the temporal and spatial change conditions of different temperature grades; 3) in the aspect of exploring the influence of the land coverage and the change thereof on the surface temperature, firstly, the temperature response of the surface temperature of different land coverage types in different seasons is analyzed; secondly, semi-quantitatively researching the relationship between the surface temperature and the surface temperature by introducing surface characteristic parameters; then quantitatively analyzing the contribution ratio of the land cover type to the surface temperature after the land cover type is changed; and finally, by introducing a distribution function, quantitatively analyzing heat source area information of different land coverage types in the research area.

Description

Technical Field

[0001] The present invention relates to the field of land cover change and thermal environment impact assessment, and in particular to a method for studying the impact of land use / cover change on the surface thermal environment based on LUCC. Background Art

[0002] Land use and cover change (LUCC) is an important research topic in global environmental change and sustainable development. Under the joint promotion of the International Geosphere-Biosphere Programme (IGBP) and the International Human Dimensions Programme (IHDP), LUCC research has become one of the core contents of current global environmental change research. Differences in the thermodynamic properties of land use / cover types have always been considered to be the fundamental cause of land surface temperature. Therefore, with the change of land use / cover types, the land surface temperature will also change, which will in turn affect the surface thermal environment. Therefore, understanding the impact of land use / cover changes on the surface thermal environment is very important for evaluating the thermal effect configuration of land landscape planning. In recent years, with the acceleration of my country's urbanization process and the continuous enhancement of human ability to transform the surface, the land use / cover characteristics have changed greatly in time and space. Among them, the accelerated urbanization process has caused a significant increase in the area of ​​urban construction land, a decrease in the area of ​​urban green space, an intensification of the urban heat island effect, and a decline in urban environmental quality and livability. These issues have attracted widespread attention from scholars. Academician Ye Duzheng proposed at the 21st Council of the International Science and Technology Union that equal attention should be paid to global environmental problems caused by land use and environmental problems caused by greenhouse gases. This proposal was accepted by many countries, and LUCC research has become an international research topic. Early LUCC research focused on manual surveys, classification and mapping of land use / cover types. With the introduction of agricultural location theory, land is more likely to be regarded as an economic resource, with the goal of maximizing benefits and optimal analysis of static land use / cover patterns. By the first half of the 20th century, the rapid development of industrialization and urbanization led to drastic changes in land use / cover. The main characteristics of LUCC research during this period also shifted to the establishment of LUCC research theories and models in various disciplines, and LUCC research tended to be systematic and scientific. By the end of the 20th century, LUCC research covered all fields of geography and all levels of spatial scale. At present, with the continuous deepening of LUCC research, new technical means such as remote sensing, geographic information systems, deep learning, and data mining have also played an important role in it.

[0003] Remote sensing technology must be combined with socio-economic activities to study LUCC. Remote sensing mainly obtains land cover classification results through images, and then uses geographic information system technology to combine them with natural or human characteristics for analysis. Geoscientists engaged in social sciences often analyze the mutual influence and effect of LUCC from aspects such as human behavior and policy orientation. It can be imagined that the combination of the two will form a special advantage. However, it is difficult to realize spatialization of the indispensable socio-economic factors (policies, population, income, output, consumption, trade and population migration, etc.) in some LUCC studies. Moreover, there are great difficulties in matching socio-economic data and remote sensing data in terms of time and spatial resolution.

[0004] Land surface temperature (LST) is a key parameter of the land surface system process at regional and global scales. It integrates the interaction between the land surface and the atmosphere and the energy exchange between the atmosphere and the land. As a key parameter in many basic disciplines and application fields, LST can provide information on the spatiotemporal variation of the surface energy balance state and is widely used in research fields such as numerical forecasting, global circulation models, and regional climate models. Accurate surface temperature information is not only helpful in evaluating the surface energy and hydrological balance, thermal inertia, and soil moisture, but also helps to obtain the global surface temperature and grasp its long-term changes. Because the surface temperature has significant spatial heterogeneity, it is difficult to obtain the surface temperature of the entire region using traditional observation methods based on meteorological stations, and it is difficult to grasp it on a large scale. At present, with the continuous development of remote sensing technology, the use of thermal infrared remote sensing inversion to obtain surface temperature has become a hot topic in geoscience research. The use of thermal infrared remote sensing to obtain surface temperature has the characteristics of large range, rapidity, convenience, large amount of information, and high accuracy. Therefore, it is widely used in the monitoring and research of surface temperature. Land use / cover changes can cause changes in surface temperature. Therefore, it is of great significance to explore the relationship between land use / cover changes and surface temperature.

[0005] Land use affects the regional and global environment to a large extent through land cover changes. According to the mode of energy transfer in nature, different land cover types and patterns will form different heat characteristics, radiation characteristics and cultural characteristics on the surface, which will greatly affect the regional thermal environment. In the study of the relationship between regional LUCC and thermal environment, multi-phase remote sensing images are usually used to extract land cover type information, and land surface temperature inversion is performed at the same time, and then superposition and comparison analysis is performed to explore the impact of LUCC on the thermal environment. Most researchers choose large and medium-sized cities as research areas to explore the relationship between LUCC and surface temperature, but do not consider the adaptability of typical areas. Based on this, the present invention uses remote sensing images of different phases as data sources, obtains multi-phase land cover data through remote sensing interpretation and combines it with surface temperature inversion data, and comprehensively uses RS and GIS spatial analysis technology to explore the relationship between land cover type and surface thermal environment. Response and contribution of land cover change to urban thermal environment, quantitatively reveal the relationship between land cover change and urban thermal environment, provide a basis for further studying the impact of land cover change on thermal environment, and provide a reference for optimizing urban ecological environment. Summary of the invention

[0006] The technical problem to be solved by the present invention is to provide a method for studying land cover changes and thermal environment impacts based on multi-temporal data, which has the characteristics of high and accurate classification accuracy, good versatility, reliable and practical inversion results, wide application range and strong intuitiveness.

[0007] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0008] A method for studying land cover change and thermal environment impact based on multi-temporal data, the implementation steps are as follows:

[0009] 1) Obtain multi-temporal remote sensing images of the area to be monitored, conduct land cover classification mapping research on remote sensing data of different phases of the monitored area, and count the area of ​​each land cover type. Study the change speed of land cover types in different phases of the monitored area, and set up research nodes to compare the changes in two different dynamic degrees of single land cover dynamic degree and comprehensive land cover dynamic degree of different land cover types before and after the establishment of the monitored area;

[0010] 2) Invert the surface temperature of remote sensing data of different phases in the monitored area, classify the temperature, and analyze the spatiotemporal changes of different temperature levels;

[0011] 3) In exploring the impact of land cover and its changes on surface temperature; first, analyze the temperature response of surface temperature of different land cover types in different seasons. Secondly, introduce the relationship between surface characteristic parameters and surface temperature semi-quantitatively. Then, after partially eliminating the seasonal impact of surface temperature by using impervious surface information, quantitatively analyze the contribution of land cover type changes to surface temperature. Finally, by introducing distribution functions, quantitatively analyze the heat source area information of different land cover types in the study area, and further explore the impact of land cover type changes on heat source areas through the trajectory of center of gravity movement.

[0012] As a further improvement of the above technical solution of the present invention:

[0013] The detailed steps of conducting land cover classification mapping research on remote sensing data of different phases of the monitored area in step 1) and counting the area of ​​each land cover type include:

[0014] A1) Based on the two classification systems, the USGS Land Use LandCover Classification System and the Chinese Academy of Sciences Land Resources Classification System, and combined with the image features used for interpretation and the actual land cover situation of the monitored area, the land use / land cover classification standard of the monitored area is constructed, and the land cover classification system of the study area is formulated by comprehensively considering factors such as the spatial resolution of the monitored area data and seasonal differences;

[0015] A2) Select training samples. The selection of training samples plays a decisive role in the image classification results. To ensure that the selected training samples are representative, the following requirements should be noted: the number of training samples for each type of land object must be appropriate. Too few training samples are not representative and it is difficult to ensure accuracy, while too many samples will increase the amount of calculation; the area of ​​the training samples should be large enough, and try not to choose the pixels near the boundaries between different types of land objects to ensure that the selected pixels are as pure as possible; the selected training samples should be evenly distributed in the study area;

[0016] A3) Establish interpretation marks, based on the image features of the target and its surroundings on the remote sensing image - texture, size, tone, shape, pattern and the spatial arrangement and combination rules of the target on the image, and through the relationship between the objects, through comprehensive reasoning and in-depth analysis to identify the target and obtain the interpretation marks of the area to be monitored;

[0017] A4) Add characteristic indexes. The Normalized Difference Vegetation Index (NDVI), the Modified Normalized Difference Water Index (MNDWI), the Residential Ratio Index (RRI), the Enhanced Vegetation Index (EVI) and the Digital Elevation Model (DEM) are selected as characteristic indexes and added to the classified images. These indexes can highlight the gaps between objects, increase the spectral evidence during classification, and effectively improve the classification accuracy.

[0018] A5) Select the classification method. After comparing the minimum distance method, maximum likelihood method, neural network classification method, Mahalanobis distance method, support vector machine classification method, and random forest classification method, the random forest method with better classification performance in the study area is selected for classification. This method is a machine learning method that selects training samples, randomly constructs a decision tree model, conducts a comprehensive analysis of the decision tree classification results, and votes to obtain the final classification results;

[0019] A6) Evaluate the classification accuracy. Only by analyzing the accuracy of the image classification results can the classifier determine whether the classification method is effective; the user can correctly and effectively obtain the available information in the classification results based on the accuracy of the classification results. There are generally two forms of accuracy verification: one is to make a qualitative evaluation of the classification results based on visual inspection, which is highly subjective; the other is to use the confusion matrix (Confusion matrix), also known as the error matrix (Error matrix) or ROC curve, to quantitatively evaluate the classification results;

[0020] The detailed steps of using the confusion matrix (also known as the error matrix) or the ROC curve in step A6) to quantitatively evaluate the classification results include:

[0021] A61) Calculate the overall accuracy. Overall accuracy refers to the sum of correctly classified pixels divided by the total number of pixels. The sum of correctly classified pixels is the sum of the main diagonal values ​​of the confusion matrix. The total number of pixels is equal to the sum of the pixels of all real reference sources. Get the overall accuracy;

[0022] A62) Calculate the producer precision. Producer precision refers to the ratio of the total number of pixels that are accurately classified to the total number of pixels of this type of features as independent training sample areas in each category. Get producer precision;

[0023] A63) Calculate user precision. User precision means the probability that a pixel is classified into the correct object category. Get user precision;

[0024] A64) Calculate the Kappa coefficient. The Kappa coefficient is an important indicator for evaluating accuracy. It comprehensively uses all the parameters of the confusion matrix. Get the Kappa coefficient, where X i+ and X +i represents the sum of the i column and the i row respectively, X ii represents the value on row i and column i, n represents the number of rows in the error matrix after classification, and M represents the total number of all ground object classification samples that need to be detected during classification;

[0025] In the step 1), the change speed of the land cover types in different phases of the monitored area is studied, and a research node is set up to compare the changes of two different dynamic degrees of single land cover dynamic degree and comprehensive land cover dynamic degree of different land cover types before and after the establishment of the monitored area. The detailed steps include:

[0026] B1) Calculate the dynamics of single land cover type change. The significance of single land cover dynamics is that it can clearly reflect the change range and speed of each land cover type, and can also reflect the type difference of change through comparison between types, providing a basis for subsequent driving force analysis. The dynamic degree of a certain land cover type during the monitoring period is obtained, where K is the dynamic degree of a certain land cover type during the monitoring period; u a To monitor the amount of a certain land cover type in the initial stage; U b is the number of a certain land cover type at the end of the monitoring period; T is the length of the research period. When T is set to years, the value of K is the annual change rate of a certain land cover type in the research area;

[0027] B2) Calculate the comprehensive land cover dynamics of the monitored area. The comprehensive land cover dynamics is used to describe the change process of land cover types in the monitored area within a certain time range. The comprehensive land cover dynamics is obtained, where ΔLU i is the area of ​​the i-th land cover type at the start time of monitoring; ΔLU i-jIt is the sum of the absolute values ​​of the area of ​​land cover type i converted to non-land cover type i and the area of ​​non-land cover type i converted to land cover type i during the monitoring period; T is the length of the monitoring period. When T is set to years, the value of LC is the annual change rate of comprehensive land cover in the study area;

[0028] B3) Introduce time series analysis - Dynamic Time Warping (DTW) and Autoregressive Integrated Moving Average (ARIMA) to study the changing trends and periodicity of land cover types. DTW can be used to compare the spatial distribution and change patterns of land cover types at different time points. By analyzing the similarity matrix, the main patterns and trends of land cover type changes can be identified. The time series data of each land cover type is tested for stationarity according to the ARIMA model. If it is not stationary, differential processing is performed until the series is stationary. Then, the determined ARIMA model is used to estimate the parameters and test the model of the time series data of each land cover type to analyze the periodicity of land cover type changes.

[0029] B4) Analyze the land conversion in the monitored area, including the conversion of construction land into and spatial change analysis, the conversion of cultivated land out and spatial change analysis, and the conversion of water area out and spatial change analysis;

[0030] The detailed steps of performing surface temperature inversion on remote sensing data of different phases in the monitored area in step 2), grading the temperature, and analyzing the spatiotemporal changes of different temperature levels include:

[0031] C1) Calculate the average temperature of the pixel. The principle and algorithm of remote sensing inversion of the average temperature of the pixel are mainly based on Planck's law and the thermal infrared radiation transmission equation. Under ideal conditions, the spectral radiance of a black body is a function of its thermodynamic temperature and wavelength. The radiation information received by the satellite sensor not only includes the radiation energy transmitted from the surface through the atmosphere, but also includes the energy radiated outward by the atmosphere itself to reach the sensor.

[0032] C2) to calculate the brightness temperature, and use the calibration coefficient to calibrate the DN value of the thermal infrared band. The result is that the DN value of the pixel is converted into the spectral radiation brightness of the upper boundary of the atmosphere. Get the brightness temperature of the area to be monitored, where K1 and K2 are constants, and L i is the spectral radiance of the upper boundary of the atmosphere;

[0033] C3) Calculate the surface emissivity. This paper chooses the NDVI threshold improvement algorithm proposed by the first author Qin Zhihao to invert the surface emissivity. It not only considers the natural surface, but also the two surface types of water bodies and towns. The water body is extracted by the improved Normalized Difference Water Index (MNDWI) mask and assigned a value of 0.995. Among them, the town pixels are mixed pixels composed of buildings and green vegetation. Natural surface pixels are mixed pixels composed of vegetation canopies and bare soil in different proportions;

[0034] C4) Remote sensing inversion of atmospheric water vapor content. The physical basis for inverting atmospheric water vapor content using MODIS data is to use reflected solar radiation to detect the magnitude of radiation absorption by water vapor. The total amount of water vapor can be obtained by the difference in solar radiation reflected between the absorption channel and the non-absorption channel of water vapor. According to w = (α - ln w / β) 2 With τ w =ρ(19) / β(2) to obtain the water vapor content of the atmosphere;

[0035] C5) to invert the surface temperature. This paper inverts the surface temperature based on the Jiménez-Munoz single-channel algorithm. Since this method was proposed, it has achieved good results in different research areas. When the atmospheric water vapor content in the air is 0.5 g cm -2 ~2g·cm -2 When the temperature is between 0.04 and 0.17°, the error of inverting the surface temperature using the single-channel method should be between 1K and 2K.

[0036] C6) To classify the surface temperature, the present invention uses remote sensing images at different times. To increase comparability, the surface temperature inversion data of different time series are used more accurately to analyze the spatiotemporal changes of the surface temperature in the study area. The present invention uses the mean value and standard deviation to classify the surface temperature. By calculating the temperature mean and standard deviation of the entire study area, the temperature is classified according to the mean value ± standard deviation. For example, the area where the LST is between the mean value ± 0.5 standard deviation and the mean value ± 0.5 standard deviation is classified as 3. This method can eliminate the influence of geographical factors on temperature by comparing the relative highs and lows in regional temperatures, thereby eliminating interference factors for subsequent studies on the relationship between land cover changes and temperature.

[0037] The step 3) specifically refers to using the previous land cover classification results, combining the surface temperature map of the study area at the same time, and using the ArcGIS spatial analysis module to overlay the surface temperature map with the land type map to obtain the maximum, minimum, average, and standard deviation of the surface temperature of the land cover types in different seasons in the monitored area, and performing a single sample T test to test whether the average surface temperature between different land use types reaches a significant level, wherein the significance level is set to 0.05;

[0038] In step 3), introducing the semi-quantitative study of the relationship between the surface characteristic parameters and the surface temperature specifically refers to introducing some specific indexes that can represent the surface information to explore the relationship between them and the surface temperature. Among them, the Normalized Difference Vegetation Index (NDVI) is selected as a representative for vegetation research, the Modified Normalized Difference Water Index (MNDWI) is selected as a representative for water research, and the Normalized Building Index (NDBI) is selected as a representative for building research;

[0039] After partially eliminating the seasonal impact of the surface temperature using the impervious surface information in step 3), the contribution of the land cover type change to the surface temperature is quantitatively analyzed. Specifically, in order to further quantify the urban land cover type and its relationship with the thermal environment, the present invention divides the impervious surface coverage into categories representing different levels of urban development. The pixels are grouped using the equidistant classification method, and each pixel is divided into 10 equal levels (0-10%, 10-20%, 30%, ..., 90-100%) according to its abundance value. Grouping the land surface temperature within each impervious surface grouping in this way helps to analyze the land surface temperature in different groups in detail, and can analyze and compare the spatial distribution of urban thermal environments in areas with different urban densities;

[0040] In step 3), the distribution function is introduced to quantitatively analyze the heat source area information of different land cover types in the study area, and the impact of land cover type changes on the heat source area is further explored through the trajectory of center of gravity movement. Specifically, the concept of distribution index (DI) is introduced to quantify the contribution of each land cover type to the thermal effect and determine its nature as a heat source or sink. Get, where S i represents the area of ​​the i-th land cover type; S ih represents the area of ​​high temperature (including high temperature and extremely high temperature) in the i-th land cover type; S represents the total area of ​​the study area; S hIt indicates the area of ​​high temperature (including high temperature and extremely high temperature) in the study area. If DI>1, it means that this land cover type has a great contribution to the thermal environment of the study area and should be defined as a heat source; if DI<1, it means that this land cover type has a small contribution to the thermal environment of the study area.

[0041] The present invention has the following advantages:

[0042] 1. After comparing the minimum distance method, the maximum likelihood method, the neural network classification method, the Mahalanobis distance method, the support vector machine classification method, and the random forest classification method, the present invention selects the random forest method with better classification performance for classification, and the obtained classification result has high classification accuracy and accuracy.

[0043] 2. The present invention uses single land cover dynamics and comprehensive land cover dynamics to describe the changes of a certain land cover type in the study area within a certain time range. It can clearly reflect the change amplitude and speed of each land cover type, and can also reflect the difference in the types of changes through comparison between types, providing a basis for the subsequent driving force analysis.

[0044] 3. When performing remote sensing inversion of atmospheric water vapor content, traditionally, the atmospheric water vapor content is obtained based on empirical formulas or atmospheric radiation transmission model simulations, but the results are limited to point data, which makes it difficult to reflect the spatial distribution differences of atmospheric water vapor content. The present invention uses remote sensing to invert atmospheric water vapor content to make up for this shortcoming well, and can accurately invert the atmospheric water vapor content to the pixel scale, which can greatly improve the accuracy of surface temperature inversion.

[0045] 4. The present invention uses the Jiménez-Munoz single-channel algorithm to invert the surface temperature. Since this method was proposed, it has achieved good results in different research areas, and the inversion results are reliable and practical.

[0046] 5. The present invention uses remote sensing images from different times and uses the surface temperature normalization method to unify the surface temperature distribution range to between 0 and 1, which increases comparability and more accurately uses the surface temperature inversion data of different time series to analyze the spatiotemporal changes of the surface temperature in the study area. It is more intuitive and conducive to subsequent analysis.

[0047] 6. In order to quantify the relationship between vegetation, water area, impervious surface and surface temperature and avoid too general and broad analysis, this invention introduces some specific indexes that can represent surface information to explore their relationship with surface temperature. Among them, the Normalized Difference Vegetation Index (NDVI) is selected as a representative for vegetation research, the Modified Normalized Difference Water Index (MNDWI) is selected as a representative for water area research, and the Normalized Building Index (NDBI) is selected as a representative for building research.

[0048] 7. In order to eliminate the impact of year and seasonal differences on surface temperature when using data from different years and seasons in the study, the present invention uses surface temperature normalization, and only uses changes from summer to summer and winter to winter in the study period, and uses average temperature to participate in the calculation, which greatly reduces the impact of data year and seasonal differences on surface temperature. By introducing impermeable surface information, the impact of year and seasonal differences on surface temperature is simulated and analyzed, which eliminates the impact of seasonal differences in the analysis of the impact of land cover changes on surface temperature. It can be widely used in agriculture, environmental protection, water conservancy, land and other departments. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 Schematic diagram of the basic process of an embodiment of the present invention.

[0050] Figure 2 The land cover classification map of Pingtan Island obtained by the embodiment of the present invention, wherein: A) is the land cover classification map of 2001-08-12; B) is the land cover classification map of 2010-07-20; C) is the land cover classification map of 2013-07-12; D) is the land cover classification map of 2000-12-13; E) is the land cover classification map of 2010-12-27; F) is the land cover classification map of 2013-12-03.

[0051] Figure 3 The surface temperature inversion result map of Pingtan Island from 2000 to 2013 obtained by the embodiment of the present invention, wherein: A) is the surface temperature inversion result map of 2001-08-12; B) is the surface temperature inversion result map of 2010-07-20; C) is the surface temperature inversion result map of 2013-07-12; D) is the surface temperature inversion result map of 2000-12-13; E) is the surface temperature inversion result map of 2010-12-27; F) is the surface temperature inversion result map of 2013-12-03.

[0052] Figure 4The surface temperature level map of Pingtan Island from 2000 to 2013 according to an embodiment of the present invention, wherein: A) is the surface temperature level map of 2001-08-12; B) is the surface temperature level map of 2010-07-20; C) is the surface temperature level map of 2013-07-12; D) is the surface temperature level map of 2000-12-13; E) is the surface temperature level map of 2010-12-27; F) is the surface temperature level map of 2013-12-03.

[0053] Figure 5 The scatter plots between NDVI and MNDWI and the surface temperature obtained in the embodiments of the present invention, wherein: A) is a scatter plot between NDVI and MNDWI and the surface temperature on 2000-12-13; B) is a scatter plot between NDVI and MNDWI and the surface temperature on 2010-12-27; C) is a scatter plot between NDVI and MNDWI and the surface temperature on 2013-12-03; D) is a scatter plot between NDVI and MNDWI and the surface temperature on 2001-08-12; E) is a scatter plot between NDVI and MNDWI and the surface temperature on 2010-07-20; F) is a scatter plot between NDVI and MNDWI and the surface temperature on 2013-07-12.

[0054] Figure 6 The scatter plot between MNDWI and the surface temperature obtained in the embodiment of the present invention, wherein: A) is a scatter plot between MNDWI and the surface temperature on 2000-12-13; B) is a scatter plot between MNDWI and the surface temperature on 2010-12-27; C) is a scatter plot between MNDWI and the surface temperature on 2013-12-03; D) is a scatter plot between MNDWI and the surface temperature on 2001-08-12; E) is a scatter plot between MNDWI and the surface temperature on 2010-07-20; F) is a scatter plot between MNDWI and the surface temperature on 2013-07-12.

[0055] Figure 7 The scatter plot between NDBI and land surface temperature obtained in the embodiment of the present invention. Among them: A) is a scatter plot between NDBI and land surface temperature on 2000-12-13; B) is a scatter plot between NDBI and land surface temperature on 2010-12-27; C) is a scatter plot between NDBI and land surface temperature on 2013-12-03; D) is a scatter plot between NDBI and land surface temperature on 2001-08-12; E) is a scatter plot between NDBI and land surface temperature on 2010-07-20; F) is a scatter plot between NDBI and land surface temperature on 2013-07-12.

[0056] Figure 8The impervious surface coverage map of Pingtan Island from 2000 to 2013 obtained according to the embodiments of the present invention, wherein: A) is the impervious surface coverage map of 2001-08-12; B) is the impervious surface coverage map of 2010-07-20; C) is the impervious surface coverage map of 2013-07-12; D) is the impervious surface coverage map of 2000-12-13; E) is the impervious surface coverage map of 2010-12-27; F) is the impervious surface coverage map of 2013-12-03.

[0057] Fig. 9 This is a relationship diagram between different average surface temperatures and average impervious surface proportions obtained in an embodiment of the present invention.

[0058] Fig.10 DI values ​​and high temperature area maps of different land types on Pingtan Island obtained in an embodiment of the present invention, wherein: A) is the DI value and high temperature area map of different land types on 2000-12-13; B) is the DI value and high temperature area map of different land types on 2010-12-27; C) is the DI value and high temperature area map of different land types on 2013-12-03; D) is the DI value and high temperature area map of different land types on 2001-08-12; E) is the DI value and high temperature area map of different land types on 2010-07-20; F) is the DI value and high temperature area map of different land types on 2013-07-12.

[0059] Fig.11 The spatial distribution map of the heat source area of ​​Pingtan Island obtained according to the embodiments of the present invention, wherein: A) is the spatial distribution map of the heat source area on 2001-08-12; B) is the spatial distribution map of the heat source area on 2010-07-20; C) is the spatial distribution map of the heat source area on 2013-07-12; D) is the spatial distribution map of the heat source area on 2000-12-13; E) is the spatial distribution map of the heat source area on 2010-12-27; F) is the spatial distribution map of the heat source area on 2013-12-03. DETAILED DESCRIPTION

[0060] like Figure 1 As shown, the implementation steps of the land cover change and thermal environment impact research method based on multi-temporal data in this embodiment are as follows:

[0061] 1) Obtain multi-temporal remote sensing images of the area to be monitored, conduct land cover classification mapping research on remote sensing data of different phases of the monitored area, and count the area of ​​each land cover type. Study the change speed of land cover types in different phases of the monitored area, and set up research nodes to compare the changes in two different dynamic degrees of single land cover dynamic degree and comprehensive land cover dynamic degree of different land cover types before and after the establishment of the monitored area;

[0062] 2) Invert the surface temperature of remote sensing data of different phases in the monitored area, classify the temperature, and analyze the spatiotemporal changes of different temperature levels;

[0063] 3) In exploring the impact of land cover and its changes on surface temperature; first, analyze the temperature response of surface temperature of different land cover types in different seasons. Secondly, introduce surface characteristic parameters to semi-quantitatively study the relationship with surface temperature. Then, after using impervious surface information to partially eliminate the seasonal impact of surface temperature, quantitatively analyze the contribution of land cover type changes to surface temperature. Finally, by introducing distribution functions, quantitatively analyze the heat source area information of different land cover types in the study area, and further explore the impact of land cover type changes on heat source areas through the trajectory of center of gravity movement;

[0064] As a further improvement of the above technical solution of the present invention:

[0065] The detailed steps of conducting land cover classification mapping research on remote sensing data of different phases of the monitored area in step 1) and counting the area of ​​each land cover type include:

[0066] A1) Based on the two classification systems, the USGS Land Use LandCover Classification System and the Chinese Academy of Sciences Land Resources Classification System, and combined with the image features used for interpretation and the actual land cover situation of the monitored area, the land use / land cover classification standard of the monitored area is constructed, and the land cover classification system of the study area is formulated by comprehensively considering factors such as the spatial resolution of the monitored area data and seasonal differences;

[0067] A2) Select training samples. The selection of training samples plays a decisive role in the image classification results. To ensure that the selected training samples are representative, the following requirements should be noted: the number of training samples for each type of land object must be appropriate. Too few training samples are not representative and it is difficult to ensure accuracy, while too many samples will increase the amount of calculation; the area of ​​the training samples should be large enough, and try not to choose the pixels near the boundaries between different types of land objects to ensure that the selected pixels are as pure as possible; the selected training samples should be evenly distributed in the study area;

[0068] A3) Establish interpretation marks, based on the image features of the target and its surroundings on the remote sensing image - texture, size, tone, shape, pattern and the spatial arrangement and combination rules of the target on the image, and through the relationship between the objects, through comprehensive reasoning and in-depth analysis to identify the target and obtain the interpretation marks of the area to be monitored;

[0069] A4) Add characteristic indexes. The Normalized Difference Vegetation Index (NDVI), the Modified Normalized Difference Water Index (MNDWI), the Residential Ratio Index (RRI), the Enhanced Vegetation Index (EVI) and the Digital Elevation Model (DEM) are selected as characteristic indexes and added to the classified images. These indexes can highlight the gaps between objects, increase the spectral evidence during classification, and effectively improve the classification accuracy.

[0070] A5) Select the classification method. After comparing the minimum distance method, maximum likelihood method, neural network classification method, Mahalanobis distance method, support vector machine classification method, and random forest classification method, the random forest method with better classification performance in the study area is selected for classification. This method is a machine learning method that selects training samples, randomly constructs a decision tree model, conducts a comprehensive analysis of the decision tree classification results, and votes to obtain the final classification results;

[0071] A6) Evaluate the classification accuracy. Only by analyzing the accuracy of the image classification results can the classifier determine whether the classification method is effective; the user can correctly and effectively obtain the available information in the classification results based on the accuracy of the classification results. There are generally two forms of accuracy verification: one is to make a qualitative evaluation of the classification results based on visual inspection, which is highly subjective; the other is to use the confusion matrix (Confusion matrix), also known as the error matrix (Error matrix) or ROC curve, to quantitatively evaluate the classification results;

[0072] The detailed steps of using the confusion matrix (also known as the error matrix) or the ROC curve in step A6) to quantitatively evaluate the classification results include:

[0073] A61) Calculate the overall accuracy. Overall accuracy refers to the sum of correctly classified pixels divided by the total number of pixels. The sum of correctly classified pixels is the sum of the main diagonal values ​​of the confusion matrix. The total number of pixels is equal to the sum of the pixels of all real reference sources. Get the overall accuracy;

[0074] A62) Calculate the producer precision. Producer precision refers to the ratio of the total number of pixels that are accurately classified to the total number of pixels of this type of features as independent training sample areas in each category. Get producer precision;

[0075] A63) Calculate user precision. User precision means the probability that a pixel is classified into the correct object category. Get user precision;

[0076] A64) Calculate the Kappa coefficient. The Kappa coefficient is an important indicator for evaluating accuracy. It comprehensively uses all the parameters of the confusion matrix. Get the Kappa coefficient, where X i+ and X +i represents the sum of the i column and the i row respectively, X ii represents the value on row i and column i, n represents the number of rows in the error matrix after classification, and M represents the total number of all ground object classification samples that need to be detected during classification;

[0077] In the step 1), the change speed of the land cover types in different phases of the monitored area is studied, and a research node is set up to compare the changes of two different dynamic degrees of single land cover dynamic degree and comprehensive land cover dynamic degree of different land cover types before and after the establishment of the monitored area. The detailed steps include:

[0078] B1) Calculate the dynamics of single land cover type change. The significance of single land cover dynamics is that it can clearly reflect the change range and speed of each land cover type, and can also reflect the type difference of change through comparison between types, providing a basis for subsequent driving force analysis. The dynamic degree of a certain land cover type during the monitoring period is obtained, where K is the dynamic degree of a certain land cover type during the monitoring period; U a To monitor the amount of a certain land cover type in the initial stage; U b is the number of a certain land cover type at the end of the monitoring period; T is the length of the research period. When T is set to years, the value of K is the annual change rate of a certain land cover type in the research area;

[0079] B2) Calculate the comprehensive land cover dynamics of the monitored area. The comprehensive land cover dynamics is used to describe the change process of land cover types in the monitored area within a certain time range. The comprehensive land cover dynamics is obtained, where ΔLU i is the area of ​​the i-th land cover type at the start time of monitoring; ΔLU i-jIt is the sum of the absolute values ​​of the area of ​​land cover type i converted to non-land cover type i and the area of ​​non-land cover type i converted to land cover type i during the monitoring period; T is the length of the monitoring period. When T is set to years, the value of LC is the annual change rate of comprehensive land cover in the study area;

[0080] B3) Introduce time series analysis - Dynamic Time Warping (DTW) and Autoregressive Integrated Moving Average (ARIMA) to study the changing trends and periodicity of land cover types. DTW can be used to compare the spatial distribution and change patterns of land cover types at different time points. By analyzing the similarity matrix, the main patterns and trends of land cover type changes can be identified. The time series data of each land cover type is tested for stationarity according to the ARIMA model. If it is not stationary, differential processing is performed until the series is stationary. Then, the determined ARIMA model is used to estimate the parameters and test the model of the time series data of each land cover type to analyze the periodicity of land cover type changes.

[0081] B4) Analyze the land conversion in the monitored area, including the conversion of construction land into and spatial change analysis, the conversion of cultivated land out and spatial change analysis, and the conversion of water area out and spatial change analysis;

[0082] The detailed steps of performing surface temperature inversion on remote sensing data of different phases in the monitored area in step 2), grading the temperature, and analyzing the spatiotemporal changes of different temperature levels include:

[0083] C1) Calculate the average temperature of the pixel. The principle and algorithm of remote sensing inversion of the average temperature of the pixel are mainly based on Planck's law and the thermal infrared radiation transmission equation. Under ideal conditions, the spectral radiance of a black body is a function of its thermodynamic temperature and wavelength. The radiation information received by the satellite sensor not only includes the radiation energy transmitted from the surface through the atmosphere, but also includes the energy radiated outward by the atmosphere itself to reach the sensor.

[0084] C2) to calculate the brightness temperature, and use the calibration coefficient to calibrate the DN value of the thermal infrared band. The result is that the DN value of the pixel is converted into the spectral radiation brightness of the upper boundary of the atmosphere. Get the brightness temperature of the area to be monitored, where K1 and K2 are constants, and L i is the spectral radiance of the upper boundary of the atmosphere;

[0085] C3) Calculate the surface emissivity. This paper chooses the NDVI threshold improvement algorithm proposed by the first author Qin Zhihao to invert the surface emissivity. It not only considers the natural surface, but also the two surface types of water bodies and towns. The water body is extracted by the improved Normalized Difference Water Index (MNDWI) mask and assigned a value of 0.995. Among them, the town pixels are mixed pixels composed of buildings and green vegetation. Natural surface pixels are mixed pixels composed of vegetation canopies and bare soil in different proportions;

[0086] C4) Remote sensing inversion of atmospheric water vapor content. The physical basis for inverting atmospheric water vapor content using MODIS data is to use reflected solar radiation to detect the magnitude of radiation absorption by water vapor. The total amount of water vapor can be obtained by the difference in solar radiation reflected between the absorption channel and the non-absorption channel of water vapor. According to w = (α - lnτ w / β) 2 With τ w =ρ(19) / β(2) to obtain the water vapor content of the atmosphere;

[0087] C5) to invert the surface temperature. This paper inverts the surface temperature based on the Jiménez-Munoz single-channel algorithm. Since this method was proposed, it has achieved good results in different research areas. When the atmospheric water vapor content in the air is 0.5 g cm -2 ~2g·cm -2 When the temperature is between 0.04 and 0.17°, the error of inverting the surface temperature using the single-channel method should be between 1K and 2K.

[0088] C6) To classify the surface temperature, the present invention uses remote sensing images at different times. To increase comparability, the surface temperature inversion data of different time series are used more accurately to analyze the spatiotemporal changes of the surface temperature in the study area. The present invention uses the mean value and standard deviation to classify the surface temperature. By calculating the temperature mean and standard deviation of the entire study area, the temperature is classified according to the mean value ± standard deviation. For example, the area where the LST is between the mean value ± 0.5 standard deviation and the mean value ± 0.5 standard deviation is classified as 3. This method can eliminate the influence of geographical factors on temperature by comparing the relative highs and lows in regional temperatures, thereby eliminating interference factors for subsequent studies on the relationship between land cover changes and temperature.

[0089] The step 3) specifically refers to using the previous land cover classification results, combining the surface temperature map of the study area at the same time, and using the ArcGIS spatial analysis module to overlay the surface temperature map with the land type map to obtain the maximum, minimum, average, and standard deviation of the surface temperature of the land cover types in different seasons in the monitored area, and performing a single sample T test to test whether the average surface temperature between different land use types reaches a significant level, wherein the significance level is set to 0.05;

[0090] In step 3), introducing the semi-quantitative study of the relationship between the surface characteristic parameters and the surface temperature specifically refers to introducing some specific indexes that can represent the surface information to explore the relationship between them and the surface temperature. Among them, the Normalized Difference Vegetation Index (NDVI) is selected as a representative for vegetation research, the Modified Normalized Difference Water Index (MNDWI) is selected as a representative for water research, and the Normalized Building Index (NDBI) is selected as a representative for building research;

[0091] After partially eliminating the seasonal impact of the surface temperature using the impervious surface information in step 3), the contribution of the land cover type change to the surface temperature is quantitatively analyzed. Specifically, in order to further quantify the urban land cover type and its relationship with the thermal environment, the present invention divides the impervious surface coverage into categories representing different levels of urban development. The pixels are grouped using the equidistant classification method, and each pixel is divided into 10 equal levels (0-10%, 10-20%, 30%, ..., 90-100%) according to its abundance value. Grouping the land surface temperature within each impervious surface grouping in this way helps to analyze the land surface temperature in different groups in detail, and can analyze and compare the spatial distribution of urban thermal environments in areas with different urban densities;

[0092] In step 3), the distribution function is introduced to quantitatively analyze the heat source area information of different land cover types in the study area, and the impact of land cover type changes on the heat source area is further explored through the trajectory of center of gravity movement. Specifically, the concept of distribution index (DI) is introduced to quantify the contribution of each land cover type to the thermal effect and determine its nature as a heat source or sink. Get, where S i represents the area of ​​the i-th land cover type; S ih represents the area of ​​high temperature (including high temperature and extremely high temperature) in the i-th land cover type; S represents the total area of ​​the study area; S hIt indicates the area of ​​high temperature (including high temperature and extremely high temperature) in the study area. If DI>1, it means that this land cover type has a great contribution to the thermal environment of the study area and should be defined as a heat source; if DI<1, it means that this land cover type has a small contribution to the thermal environment of the study area.

[0093] The present invention uses Pingtan as a specific embodiment to prove that the method of the present invention is feasible:

[0094] 1. Research on land cover types from remote sensing data at different time phases:

[0095] The main remote sensing data used in the embodiment of the present invention include: 6 summer and winter remote sensing image data of the study area from 2000 to 2013, including Landsat-5, Landsat-7 and Landsat-8 data, of which 2 remote sensing image data in 2010 use Landsat-5 data, and 2 remote sensing images around 2000 use Landsat-5 and Landsat-7 data. The 2 scenes in 2013 use Landsat-8 data, and the remote sensing images are from the United States Geological Survey (https: / / earthexplorer.usgs.gov / ). The MODIS data used in the embodiment of the present invention for inversion of atmospheric water vapor content is from the website (https: / / ladsweb.modaps.eosdis.nasa.gov / ), and the imaging time is selected to be closest to the imaging time of the above-mentioned Landsat images. The DEM data with a spatial resolution of 30 meters used in the embodiment of the present invention is from the Geospatial Data Cloud (http: / / www.gscloud.cn / ).

[0096] 1 Pingtan Island Remote Sensing Image Classification

[0097] At present, the most commonly used land use classification systems in research are the USGS Land Use Land Cover Classification System and the Chinese Academy of Sciences Land Resources Classification System. This paper mainly constructs the land use / cover classification standard of the study area based on these two classification systems, and combines the image features used for interpretation and the actual land cover of the study area, and divides the land cover types of Pingtan Island into 4 primary categories, namely: construction land, forest land, water area, and cultivated land. 3 secondary categories, respectively: sandy land, bare land, and tidal flat.

[0098] 1.1 Pingtan Island training sample selection and interpretation mark establishment

[0099] The selection of training samples plays a decisive role in the image classification results. To ensure that the selected training samples are representative, the following requirements should be noted: the number of training samples for each type of land object must be appropriate. Too few training samples are not representative and it is difficult to ensure accuracy, while too many samples will increase the amount of calculation; the area of ​​the training samples should be large enough, and try not to select the pixels near the boundaries between different land object categories to ensure that the selected pixels are as pure as possible; the selected training samples should be evenly distributed in the study area. The present invention identifies the target and its surrounding image features on the remote sensing image - texture, size, tone, shape, pattern, etc., as well as the spatial arrangement and combination rules of the target on the image, and through the relationship between the land objects, through comprehensive reasoning and in-depth analysis to establish the interpretation mark.

[0100] 1.2 Selection of Pingtan Island classification method

[0101] After comparing the minimum distance method (Minimum Distance), the maximum likelihood method (MaximumLikelihood), the neural network classification method (Neural Net), the Mahalanobis distance method (Mahalanobis Distance), the support vector machine classification method (Support Vector Machine), and the random forest classification method (Random Forest), the present invention selects the random forest method with better classification performance in the study area for classification. This method is a machine learning method, which selects training samples, randomly constructs a decision tree model, conducts a comprehensive analysis of the decision tree classification results, and obtains the final classification result by voting. The present invention uses the ENVI 5.3 random forest remote sensing classification algorithm plug-in (Random Forest Classification) to classify land cover, and classifies remote sensing images of Pingtan Island in different phases from 2000 to 2013. The classification results are as follows Figure 2 Table 1 shows the land cover area results of six periods of data from 2000 to 2013.

[0102] Table 1 Statistics of land cover types and areas in Pingtan Island

[0103]

[0104] 1.3 Classification Accuracy Evaluation

[0105] Accuracy evaluation refers to the verification of image classification results with standard data or measured data to determine the accuracy of the classification process and classification results. Only by performing accuracy analysis on image classification results can the classifier determine whether the classification method is effective; only then can the user correctly and effectively obtain the available information in the classification results based on the accuracy of the classification results. Accuracy verification generally takes two forms: one is mainly based on visual inspection, which has a qualitative evaluation of the classification results and is highly subjective; the other is to use the confusion matrix (Confusion matrix), also known as the error matrix (Error matrix) or ROC curve, to quantitatively evaluate the classification results, among which the confusion matrix is ​​more commonly used. The ROC curve represents the classification accuracy in the form of a graph, and its advantage is that it is easy to understand. The confusion matrix is ​​usually a matrix with N rows and N columns, where the rows represent the category data obtained by remote sensing data classification, and the columns represent the reference data. This paper selects overall accuracy, producer accuracy, user accuracy, and Kappa coefficient as the basic statistical parameters of the confusion matrix to test the remote sensing classification accuracy.

[0106] (1) Overall accuracy: refers to the sum of correctly classified pixels divided by the total number of pixels, where the sum of correctly classified pixels is the sum of the main diagonal values ​​of the confusion matrix, and the total number of pixels is equal to the sum of the pixels of all real reference sources. Its expression is:

[0107]

[0108] (2) Producer precision: refers to the ratio of the total number of pixels that are accurately classified to the total number of pixels of this type of features as independent training samples in each category. Its expression is:

[0109]

[0110] (3) User precision: It indicates the probability that a pixel is classified into the correct object category, and its expression is:

[0111]

[0112] (4) Kappa coefficient: As an important indicator for evaluating accuracy, it comprehensively uses all parameters of the confusion matrix and its expression is:

[0113]

[0114] Where X i+ and X +i represents the sum of the i column and the i row respectively, X iirepresents the value on the i-th row and i-th column, n represents the number of rows of the error matrix after classification, and M represents the total number of all classified samples that need to be detected during classification.

[0115] The present invention selects classification samples based on the stratified random sampling method and evaluates the classification results of the random forest method. The accuracy evaluation of the classification results mainly includes: the original land use status data; high spatial resolution remote sensing data. Among them, the classification results of remote sensing images in 2000 refer to the 2000 Pingtan Island land use status map, and use the visual interpretation method and the original image for accuracy verification. For the classification results of remote sensing images in 2010, refer to the 2010 Pingtan Island land use status map, and combine with high spatial resolution images (2009 Quickbird data) for accuracy verification. The classification results of remote sensing images in 2013 are combined with Google images to verify accuracy using the visual interpretation method. It can be seen from Tables 2, 3, and 4 that the user accuracy and producer accuracy of water and forest land are relatively high, and the user accuracy of bare land and tidal flats is low. The overall classification accuracy of the three phases of data is greater than 83.88%, and the Kappa coefficient is greater than 0.8099.

[0116] Table 2 Land cover classification error matrix on December 23, 2000

[0117]

[0118] Table 3 Land cover classification error matrix on December 20, 2010

[0119]

[0120] Table 4 Land cover classification error matrix on December 30, 2013

[0121]

[0122] 1.4 Speed ​​of land cover change in Pingtan Island

[0123] The dynamics of land cover change can quantitatively reflect the speed and degree of regional land cover change as a whole, and plays a positive role in comparing regional differences in land cover change and predicting future land cover change trends. It includes single land cover dynamics and comprehensive land cover dynamics.

[0124] 1.4.1 Dynamics of changes in single land cover types on Pingtan Island

[0125] The single land cover dynamics is used to describe the changes of a certain land cover type in the study area within a certain time range. The significance of the single land cover dynamics is that it can clearly reflect the change range and speed of each land cover type, and can also reflect the type differences of change through comparison between types, providing a basis for subsequent driving force analysis. The expression is:

[0126]

[0127] Where K is the dynamic degree of a certain land cover type during the study period; U a is the number of a certain land cover class at the beginning of the study; U b is the number of a certain land cover type at the end of the study period; T is the length of the study period. When the period of T is set to , the value of K is the annual change rate of a certain land cover type in the study area. According to the above formula, the dynamic degree of change of a single land cover type in the study area is calculated. The calculation shows that the interannual change rate of construction land from 2000 to 2013 is greater than 4.35%, and the change rate from 2010 to 2013 is 23.48%, indicating that the construction land has increased significantly in the three years after the establishment of the Pingtan Comprehensive Experimental Zone. The interannual change rate of bare land from 2000 to 2013 is greater than 54%, indicating that bare land has continued to increase in these 13 years. The interannual change rate of water area from 2000 to 2013 is greater than 5.18%, indicating that the area of ​​water area has continued to decrease in these 13 years. Among them, the interannual change rate of cultivated land from 2000 to 2013 is greater than 0.66%, indicating that the area of ​​cultivated land has continued to decrease in these 13 years. The interannual change rate of forest land from 2000 to 2013 was not significant, and the change rate was less than 1.76%, which shows that the change of forest land in the past 13 years is very small, which is closely related to the local policy of protecting forest land.

[0128] 1.4.2 Pingtan Island's comprehensive land cover dynamics

[0129] The comprehensive land cover dynamics is used to describe the change process of land cover types in a certain study area within a certain time range. Its expression is:

[0130]

[0131] Where ΔLU i-j ΔLU is the sum of the absolute values ​​of the area of ​​land cover type i converted to non-i land cover type and the area of ​​non-i land cover type converted to i land cover type during the monitoring period. i is the area of ​​the i-th land cover type at the monitoring start time; T is the length of the monitoring period. When the period of T is set to years, the value of LC is the annual change rate of the comprehensive land cover in the study area.

[0132] 2. Remote sensing inversion of Pingtan Island surface temperature:

[0133] Based on Planck's law, under ideal conditions, the spectral radiance of a black body is a function of its thermodynamic temperature and wavelength. The radiation information received by the satellite sensor not only includes the radiation energy transmitted from the surface through the atmosphere, but also includes the energy radiated outward by the atmosphere itself to the sensor. Planck's blackbody radiation law and the thermal infrared radiation transfer equation jointly describe the transmission of thermal infrared radiation, which is also the basis of the surface temperature inversion algorithm.

[0134] 2.1 Single-channel algorithm for inversion of surface temperature

[0135] The Jiménez-Munoz single-channel algorithm is proposed for retrieving land surface temperature from a single thermal infrared band. The general equation for its calculation is:

[0136] T s =γ[(ψ1L sen +ψ2) / ε+ψ3]+δ

[0137] Where ε represents the surface emissivity, T S Indicates the brightness temperature corresponding to the radiant brightness of the star, L sen is the radiation intensity measured by the remote sensor at the satellite altitude, ψ1, ψ2 and ψ3 are atmospheric function parameters, and δ is a constant used to correct system deviations. There are different calculation methods for different types of Landsat data.

[0138] 2.2 Brightness temperature estimation

[0139] The DN value of the thermal infrared band is calibrated by the calibration coefficient. The result is that the DN value of the pixel is converted into the spectral radiation brightness of the upper boundary of the atmosphere. Then the brightness temperature on the star is inferred by the Planck equation. The calculation formula is as follows:

[0140]

[0141] In the formula, K1 and K2 are constants, L i is the spectral radiation brightness of the upper boundary of the atmosphere.

[0142] 2.3 Pingtan Island surface emissivity

[0143] Surface emissivity is a necessary parameter in the process of surface temperature inversion, but it is quite difficult to obtain real-time surface emissivity data. Many scholars have completed remote sensing estimation of surface emissivity through the Normalized Different Vegetation Index (NDVI) and the empirical model of surface emissivity. This paper selects the NDVI threshold improvement algorithm proposed by the first author Qin Zhihao to invert the surface emissivity. It not only considers the natural surface, but also considers the two types of surface: water bodies and towns. The water body is extracted by the modified Normalized Difference Water Index (MNDWI) mask and assigned a value of 0.995. Among them, the town pixels are mixed pixels composed of buildings and green vegetation. The natural surface pixels are mixed pixels composed of different proportions of vegetation canopies and bare soil. Furthermore, we can define that when NDVI is less than 0.05, it can be considered as a pure bare soil pixel with a vegetation coverage of 0; when NDVI is greater than 0.7, it can be considered as a pure vegetation pixel with a vegetation coverage of 1; when NDVI is between 0.05 and 0.7, it can be considered as a mixed pixel.

[0144] 2.4 Remote sensing inversion of atmospheric water vapor content in Pingtan Island

[0145] Traditionally, atmospheric water vapor content is mostly obtained based on empirical formulas or atmospheric radiation transfer model simulations, but the results are limited to point data, and it is difficult to reflect the spatial distribution differences of atmospheric water vapor content. Using remote sensing to invert atmospheric water vapor content can well make up for this shortcoming, and can accurately determine the atmospheric water vapor content to the pixel scale, which can greatly improve the accuracy of surface temperature inversion. The physical basis for inverting atmospheric water vapor content from MODIS data is to use reflected solar radiation to detect the amount of radiation absorption by water vapor. The total amount of water vapor can be obtained by the difference in solar radiation reflected between the absorption channel and the non-absorption channel of water vapor. MODIS data contains 36 bands, of which bands 17, 18 and 19 are atmospheric absorption bands, and bands 2 and 5 are atmospheric window bands. The maximum transmittance of the second band is close to 1, while the 19th band shows strong absorption of atmospheric water vapor. The MODIS data comes from the website (https: / / ladsweb.modaps.eosdis.nasa.gov / ), and the imaging time is selected to be closest to the above-mentioned Landsat image imaging time. The DEM data with a spatial resolution of 30 meters used in the embodiment of the present invention comes from the Geospatial Data Cloud (http: / / www.gscloud.cn / ). Therefore, the reflectivity of these two bands can be compared to estimate the atmospheric water vapor content. The calculation formula is:

[0146] w=(α-lnr w / β) 2

[0147] w =ρ(19) / ρ(2)

[0148] In the formula, τ w is the atmospheric transmittance; w is the atmospheric water vapor content, ρ(2) and ρ(19) represent the apparent reflectance of the 2nd and 19th bands of the MODIS data, respectively; α and β are constants with different values ​​for different surface types. Combined with the actual surface type in the study area, the parameters of the mixed surface are α=0.02 and β=0.651.

[0149] 2.5 Pingtan Island surface temperature inversion

[0150] This paper uses the Jiménez-Munoz single-channel algorithm to invert the surface temperature. Since this method was proposed, it has achieved good results in different research areas. Since the surface temperature is inverted using remote sensing technology, the inversion results cannot be verified pixel by pixel and in real time. The following are the research results of several people engaged in temperature inversion algorithm research to verify the availability of the inversion results of this paper. According to the research results of the first author Jiménez-Munoz et al., when the atmospheric water vapor content in the air is between 0.5g·cm-2 and 2g·cm-2, the error of the surface temperature inverted using the single-channel method should be between 1K and 2K. Figure 3 Table 5 is the statistical results of the inversion of land surface temperature for the six periods of data from 2000 to 2013.

[0151] Table 5 Statistics of temperature inversion results from 2000 to 2013

[0152]

[0153] 2.6 Surface temperature classification

[0154] The surface temperature has a strong continuity. Therefore, the temperature can be more intuitively expressed after the temperature is divided into grades. At the same time, the temperature classification can also be connected with the concept of "island" in the urban heat island effect. There are many methods for dividing the surface temperature grades, and the commonly used division methods include equal interval division method, standard deviation division method and natural break point method. Among them, the equal interval division method divides the data into several equal-sized ranges, and then automatically determines its value range by determining the number of division categories. It is more subjective, the division effect is not obvious, the intuitiveness is poor, and it is not conducive to subsequent analysis; the standard deviation division method considers the mean and standard deviation of the data at the same time, and divides the data level according to their different combinations. The natural break point classification method can identify the classification interval, can group similar values ​​most appropriately, and can maximize the differences between each class. The present invention uses the natural break point classification method to classify the surface temperature, and reclassifies the surface temperature into 7 categories, namely, extremely low temperature zone, low temperature zone, relatively low temperature zone, medium temperature zone, relatively high temperature zone, high temperature zone and extremely high temperature zone. like Figure 4 As shown, it can be seen that the densely populated areas of construction land in summer are high temperature areas or extremely high temperature areas, but not obvious in winter. Table 6 statistics the area of ​​each temperature zone in Pingtan Island from 2000 to 2013. In the winter of 2000-2010, the areas of extremely low temperature areas, low temperature areas and relatively low temperature areas all showed a downward trend; the areas of medium temperature areas, relatively high temperature areas, high temperature areas and extremely high temperature areas all showed an upward trend. Among them, the extremely low temperature area decreased by 15.03km 2 , a decrease of 46.26%; the high temperature area increased by 9.97km2, an increase of 27.25%. In the winter of 2000-2013, the area of ​​extremely low temperature area and medium temperature area showed a downward trend; the low temperature area, relatively low temperature area, relatively high temperature area, high temperature area and extremely high temperature area all showed an upward trend.

[0155] Table 6 Statistics of the area of ​​different temperature zones from 2000 to 2013

[0156]

[0157] 3. Analysis of the impact of land cover change on surface temperature in Pingtan Island:

[0158] The surface temperature corresponding to different land cover types is different. Using the previous land cover classification results, combined with the surface temperature map of Pingtan Island at the same time, the ArcGIS spatial analysis module is used to overlay the surface temperature map with the land type map to obtain the maximum, minimum, average, and standard deviation of the surface temperature of different seasonal land cover types in Pingtan Island from 2000 to 2013, and a single sample T test is performed to test whether the average surface temperature between different land use types reaches a significant level, where the significance level is set to 0.05. The present invention selects the six images of the study area in summer and winter as the analysis object, and the analysis method for the remaining periods is similar.

[0159] 3.1 Relationship between surface parameters and surface temperature

[0160] Vegetation, water areas, and impervious surfaces are typical types of surface cover. In order to quantify their relationship with surface temperature and avoid too general and broad analysis, it is necessary to introduce some specific indices that can represent surface information to explore their relationship with surface temperature. Among them, the Normalized Difference Vegetation Index (NDVI) is selected as a representative for vegetation research, the Modified Normalized Difference Water Index (MNDWI) is selected as a representative for water areas, and the Normalized Difference Building Index (NDBI) is selected as a representative for buildings.

[0161] 3.1.1 Relationship between NDVI and land surface temperature

[0162] When studying the relationship between surface temperature and the Normalized Difference Vegetation Index (NDVI), it is first necessary to obtain a set of data corresponding to the surface temperature and NDVI to perform a linear correlation analysis. It is worth noting that: based on multiple experiments, the present invention decided to use only NDVI (values ​​greater than 0) for correlation analysis. Because pixels with NDVI < 0 in the image are mostly water bodies or a few outliers, and the temperature variation of the water body is small, taking it into account will reduce the correlation between the data. The present invention uses the ArcGIS random point generation tool to select the inverted temperature values ​​and NDVI values ​​of several points with the same coordinates. A scatter plot is then generated using Matlab software, and a linear correlation analysis is performed based on the collected data to obtain a correlation trend graph of the surface temperature and NDVI from 2000 to 2013 as shown below. Figure 5 -, the fitting equation is shown in Table 7. It can be seen from Table 7 that the linear law is basically consistent, and the surface temperature decreases with the increase of NDVI value, indicating that vegetation coverage can affect the surface temperature, but the relationship between the two is worse in winter than in summer. In the sparse vegetation area, the NDVI value is small, and its relative temperature is high; in the lush vegetation area, the NDVI value is large, but its relative temperature is low.

[0163] Table 7 Regression equations between NDVI and MNDWI and surface temperature

[0164]

[0165] 3.1.2 Relationship between improved water index and surface temperature

[0166] When studying the relationship between surface temperature and water area, we first need to determine the MNDWI threshold to extract water bodies, and then based on the extracted water bodies, we need to collect a set of data corresponding to surface temperature and MNDWI to perform linear correlation analysis. This paper uses the ArcGIS random point generation tool to extract surface temperature values ​​and MNDWI values ​​of several points with the same coordinates on the temperature inversion result map and the MNDWI map. Figure 6 The correlation trend of surface temperature and MNDWI from 2000 to 2013 is shown in Figure 2. Figure 6 It can be seen that when MNDWI is greater than 0, the surface temperature decreases with the increase of MNDWI value, and when MNDWI is less than 0, the surface temperature increases with the increase of MNDWI value, but the seasonal performance is not obvious. This shows that the seasonal effect of water temperature change is not significant, which may be related to the nature of the water body. Even if the value of MNDWI changes greatly in the water-dense area, its relative temperature changes are small.

[0167] 3.1.3 Relationship between normalized building index and ground surface temperature

[0168] When studying the relationship between surface temperature and NDBI, it is necessary to collect a set of data corresponding to the surface temperature and NDBI and draw a line.

[0169] The present invention uses ArcGIS random point generation tool to extract the surface temperature and NDBI values ​​of several points with the same coordinates on the temperature inversion result map and the NDBI map. Then, based on this set of data, Matlab is used to generate the corresponding scatter plot and obtain the fitting equation. Figure 7 This is a graph showing the correlation between surface temperature and NDBI from 2000 to 2013. Figure 6 It can be seen that the surface temperature increases with the increase of NDBI value, which indicates that there are large differences in surface temperature between different land cover types due to NDBI values. Among them, the NDBI value of water area is smaller, and the corresponding surface temperature is lower.

[0170] 3.2 Analysis of the impact of seasonal differences in surface temperature

[0171] The inversion result of the surface temperature represents the instantaneous surface temperature when the satellite passes over. Since the study uses data from different years and seasons, the impact of year and seasonal differences on the surface temperature needs to be taken into account when exploring the impact of land cover changes on the surface temperature. In order to eliminate this impact, the present invention mainly uses surface temperature normalization, and only uses summer to summer and winter to winter changes in the research period, and uses the average temperature to participate in the calculation, which greatly reduces the impact of data year and seasonal differences on the surface temperature. The present invention simulates and analyzes the impact of year and seasonal differences on the surface temperature by introducing impermeable surface information, and partially eliminates the impact of seasonal differences for the analysis of the impact of land cover changes on the surface temperature.

[0172] 3.2.1 Extraction of impervious surface information of Pingtan Island

[0173] Mixed pixels are widely present in remote sensing images, which seriously affects the application of remote sensing technology in the acquisition of surface information. In order to give full play to the potential of images, mixed pixel decomposition technology has been generated and continuously developed. This invention uses ENVI software plug-in to extract impervious surface coverage information. The first author Xu Hanqiu found through statistics of research papers on impervious surface information extraction that the inversion accuracy of many impervious surface information can reach more than 85%. Figure 8 It shows the inversion results of impervious surface coverage of Pingtan Island from 2000 to 2013.

[0174] 3.2.2 Classification of impervious surfaces

[0175] In order to further quantify the urban land cover types and their relationship with the thermal environment, the present invention divides the impervious surface coverage into categories representing different urban development levels, groups the pixels using the equidistant classification method, and divides each pixel into 10 equal levels (0-10%, 10-20%, 30%, ..., 90-100%) according to its abundance value. This method is used to group the land surface temperature within each impervious surface grouping range, which is helpful for detailed analysis of the land surface temperature within different groups, and can analyze and compare the spatial distribution of urban thermal environment in areas with different urban densities.

[0176] 3.2.3 Impervious surface and surface temperature scenario simulation

[0177] The surface temperature changes with changes in urban landscape patterns, urban expansion, and seasonal changes. In order to quantify the impact of these changes on the urban thermal environment, it is necessary to distinguish the impact of urban expansion on the surface temperature of each region from the impact of seasonal changes. This paper simulates the surface temperature scenario and attempts to quantitatively determine how seasonal changes and urban expansion affect urban surface temperature. Although the relationship may be different in different seasons, there is a strong linear relationship between the average impervious surface proportion of different levels and their average surface temperature. Fig. 9 The relationship between different average surface temperatures and the average impervious surface ratio in Pingtan Island from 2000 to 2013 is shown. Overall, the relationship between the average surface temperature and the average impervious surface ratio in winter is more stable than that in summer. However, the correlation between the average surface temperature and the average impervious surface ratio in summer is stronger than that in winter.

[0178] 3.3 Impact of land type conversion on surface temperature

[0179] In order to explore the impact of the mutual conversion between land types on the surface temperature, this paper selects data from the summer and winter seasons of three periods (2000-2010, 2010-2013, and 2000-2013), focusing on analyzing the impact of the conversion of cultivated land and water areas and the conversion of construction land on the surface temperature. The analysis methods of other land types are similar. In order to quantitatively describe the effect of the change in surface temperature caused by the conversion between land types, this paper introduces the concept of temperature contribution percentage, and defines the average temperature after the conversion of land type i as T1. If the average temperature when land type i does not change is T2, then T1-T2 represents the temperature difference of the conversion, so the temperature contribution percentage can be calculated by the following public expression:

[0180]

[0181] In the formula, P represents the contribution ratio, i represents the land type; where i = [construction land, forest land, water area, sandy land, bare land, cultivated land, tidal flat]; T i Indicates the average temperature of the land type being transferred.

[0182] 3.4 Impact of land type conversion on surface temperature

[0183] First, the concept of distribution index (DI) is introduced to quantify the contribution of each land cover type to the thermal effect and determine its nature as a heat source or sink. The calculation formula is as follows:

[0184]

[0185] Where: S i represents the area of ​​the i-th land cover type; S ih represents the area of ​​high temperature (including high temperature and extremely high temperature) in the i-th land cover type; S represents the total area of ​​the study area; S h It indicates the area of ​​high temperature (including high temperature and extremely high temperature) surface temperature in the monitored area. If DI>1, it means that this land cover type has a great contribution to the thermal environment of the monitored area and should be defined as a heat source; if DI<1, it means that this land cover type has a small contribution to the thermal environment of the monitored area.

[0186] Fig.10 It shows the DI values ​​and high temperature area statistics of different land cover types in six periods of data from 2000 to 2013. Fig.11The spatial distribution of heat source areas in six data periods from 2000 to 2013 was statistically analyzed. Overall, the DI value of construction land is relatively high, among which the DI values ​​of construction land in December 2000, July 2010, July 2013 and December 2013 are all greater than 1, and the value in August 2001 is also in a dominant position. It can be concluded that construction land acts as a heat source for the thermal environment of the study area in any season. On the contrary, water areas and forests do not contribute significantly to the thermal environment of Pingtan Island. In December 2000, tidal flats and bare land also acted as heat sources, but their high temperature area was very small, accounting for about 1%, which shows that their contribution to the thermal environment is very small and can be ignored.

Claims

1. A method for studying land cover change and thermal environment impact based on multi-temporal data, characterized in that Here are the steps: 1) Obtain multi-temporal remote sensing images of the area to be monitored, conduct land cover classification mapping research on remote sensing data of different phases in the monitored area, and count the area of ​​each land cover type; study the change speed of land cover types in different phases in the monitored area, and set up research nodes to compare the changes in two different dynamics of single land cover dynamics and comprehensive land cover dynamics of different land cover types before and after the establishment of the monitored area, as well as the trend and periodicity of land cover changes; 2) Invert the surface temperature of remote sensing data of different phases in the monitored area, classify the temperature, and analyze the spatiotemporal changes of different temperature levels; 3) In exploring the impact of land cover and its changes on surface temperature; first, the temperature response of surface temperature of different land cover types in different seasons is analyzed; secondly, the relationship between surface characteristic parameters and surface temperature is semi-quantitatively studied by introducing them; then, after partially eliminating the seasonal impact of surface temperature by using impervious surface information, the contribution of land cover type changes to surface temperature is quantitatively analyzed; finally, by introducing distribution functions, the heat source area information of different land cover types in the study area is quantitatively analyzed, and the impact of land cover type changes on heat source areas is further explored through the trajectory of center of gravity movement.

2. The method for studying land cover change and thermal environment impact based on multi-temporal data according to claim 1 is characterized by: The detailed steps of conducting land cover classification mapping research on remote sensing data of different phases of the monitored area in step 1) and counting the area of ​​each land cover type include: A1) Based on the two classification systems, the USGS Land Use Land Cover Classification System and the Chinese Academy of Sciences Land Resources Classification System, and combined with the image features used for interpretation and the actual land cover situation of the monitored area, the land use / cover classification standard of the monitored area is constructed, and the spatial resolution of the monitored area data and seasonal differences are comprehensively considered to develop a land cover classification system for the study area; A2) Select training samples. The selection of training samples plays a decisive role in the image classification results. To ensure that the selected training samples are representative, the following requirements should be noted: the number of training samples for each type of land object must be appropriate. Too few training samples are not representative and it is difficult to ensure accuracy, while too many samples will increase the amount of calculation; the area of ​​the training samples should be large enough, and try not to choose the pixels near the boundaries between different types of land objects to ensure that the selected pixels are as pure as possible; the selected training samples should be evenly distributed in the study area; A3) Establish interpretation marks, based on the image features of the target and its surroundings on the remote sensing image - texture, size, tone, shape, pattern and the spatial arrangement and combination rules of the target on the image, and through the relationship between the objects, through comprehensive reasoning and in-depth analysis to identify the target and obtain the interpretation marks of the area to be monitored; A4) Add characteristic indexes, select Normalized Difference Vegetation Index (NDVI), Modified Normalized Difference Water Index (MNDWI), Residential Ratio Index (RRI), Enhanced Vegetation Index (EVI) and Digital Elevation Model (DEM) as characteristic indexes and add them to the classified images. These indexes can highlight the gaps between objects, increase the spectral evidence during classification, and effectively improve the classification accuracy. A5) Select the classification method. After comparing the minimum distance method, maximum likelihood method, neural network classification method, Mahalanobis distance method, support vector machine classification method, and random forest classification method, the random forest method with better classification performance in the study area was selected for classification. This method is a machine learning method that selects training samples, randomly constructs a decision tree model, conducts a comprehensive analysis of the decision tree classification results, and votes to obtain the final classification results. A6) Evaluate the classification accuracy. Only by analyzing the accuracy of the image classification results can the classifier determine whether the classification method is effective; the user can correctly and effectively obtain the available information in the classification results based on the accuracy of the classification results. There are generally two forms of accuracy verification: one is to make a qualitative evaluation of the classification results based on visual inspection, which is highly subjective; the other is to use the confusion matrix (also known as the error matrix) or ROC curve to quantitatively evaluate the classification results. The detailed steps of using the confusion matrix (also known as the error matrix) or the ROC curve in step A6) to quantitatively evaluate the classification results include: A61) Calculate the overall accuracy. Overall accuracy refers to the sum of correctly classified pixels divided by the total number of pixels. The sum of correctly classified pixels is the sum of the main diagonal values ​​of the confusion matrix. The total number of pixels is equal to the sum of the pixels of all real reference sources. Get the overall accuracy; A62) Calculate the producer precision. Producer precision refers to the ratio of the total number of pixels that are accurately classified to the total number of pixels of this type of features as independent training sample areas in each category. Get producer precision; A63) Calculate user precision. User precision means the probability that a pixel is classified into the correct object category. Get user precision; A64) Calculate the Kappa coefficient. The Kappa coefficient is an important indicator for evaluating accuracy. It comprehensively uses all the parameters of the confusion matrix. Get the Kappa coefficient, where X i+ and X +i represents the sum of the i column and the i row respectively, X ii represents the value on row i and column i, n represents the number of rows in the error matrix after classification, and M represents the total number of all ground object classification samples that need to be detected during classification; The detailed steps of studying the change speed of land cover types in different phases of the monitored area in step 1), setting up research nodes, comparing the changes of two different dynamic degrees of single land cover dynamic degree and comprehensive land cover dynamic degree of different land cover types before and after the establishment of the monitored area, and analyzing the change trend and periodicity of the land cover types include: B1) Calculate the dynamics of single land cover type change. The significance of single land cover dynamics is that it can clearly reflect the change range and speed of each land cover type, and can also reflect the type difference of change through comparison between types, providing a basis for subsequent driving force analysis. The dynamic degree of a certain land cover type during the monitoring period is obtained, where K is the dynamic degree of a certain land cover type during the monitoring period; U a To monitor the amount of a certain land cover type in the initial stage; U b is the number of a certain land cover type at the end of the monitoring period; T is the length of the research period. When T is set to years, the value of K is the annual change rate of a certain land cover type in the research area; B2) Calculate the comprehensive land cover dynamics of the monitored area. The comprehensive land cover dynamics is used to describe the change process of land cover types in the monitored area within a certain time range. The comprehensive land cover dynamics is obtained, where ΔLU i is the area of ​​the i-th land cover type at the start time of monitoring; ΔLU i-j It is the sum of the absolute values ​​of the area of ​​land cover type i converted to non-land cover type i and the area of ​​non-land cover type i converted to land cover type i during the monitoring period; T is the length of the monitoring period. When T is set to years, the value of LC is the annual change rate of comprehensive land cover in the study area; B3) Time series analysis - Dynamic Time Warping (DTW) and Autoregressive Integrated Moving Average (ARIMA) are introduced to study the changing trend and periodicity of land cover types; DTW is used to compare the spatial distribution and change patterns of land cover types at different time points, and the main patterns and trends of land cover type changes are identified by analyzing the similarity matrix; the time series data of each land cover type are tested for stationarity according to the ARIMA model, and if they are not stationary, they are subjected to differential processing until the series is stationary, and then the determined ARIMA model is used to estimate the parameters and test the model of the time series data of each land cover type to analyze the periodicity of land cover type changes; B4) Analyze the land conversion in the monitored area, including the conversion of construction land in and its spatial change analysis, the conversion of cultivated land out and its spatial change analysis, and the conversion of water area out and its spatial change analysis.

3. The method for studying land cover change and thermal environment impact based on multi-temporal data according to claim 1 is characterized by: The detailed steps of performing surface temperature inversion on remote sensing data of different phases in the monitored area in step 2), grading the temperature, and analyzing the spatiotemporal changes of different temperature levels include: C1) Calculate the average temperature of the pixel. The principle and algorithm of remote sensing inversion of the average temperature of the pixel are mainly based on Planck's law and the thermal infrared radiation transmission equation. Under ideal conditions, the spectral radiance of a black body is a function of its thermodynamic temperature and wavelength. The radiation information received by the satellite sensor not only includes the radiation energy transmitted from the surface through the atmosphere, but also includes the energy radiated outward by the atmosphere itself to reach the sensor. C2) to calculate the brightness temperature, and use the calibration coefficient to calibrate the DN value of the thermal infrared band. The result is that the DN value of the pixel is converted into the spectral radiation brightness of the upper boundary of the atmosphere. Get the brightness temperature of the area to be monitored, where K1 and K2 are constants, and L i is the spectral radiance of the upper boundary of the atmosphere; C3) Calculate the surface emissivity and select the NDVI threshold improved algorithm proposed by the first author Qin Zhihao to invert the surface emissivity; it not only considers the natural surface, but also the two types of surface: water bodies and towns; extract water bodies through the improved Normalized Difference Water Index (MNDWI) mask and assign a value of 0.995; among them, the town pixels are mixed pixels composed of buildings and green vegetation; the natural surface pixels are mixed pixels composed of different proportions of vegetation canopies and bare soil; C4) Remote sensing inversion of atmospheric water vapor content. The physical basis for inverting atmospheric water vapor content using MODIS data is to use reflected solar radiation to detect the magnitude of radiation absorption by water vapor. The total amount of water vapor can be obtained by the difference in solar radiation reflected between the absorption channel and the non-absorption channel of water vapor. According to w = (α - lnτ w / β) 2 With τ w =ρ(19) / β(2) to obtain the atmospheric water vapor content, where: w is the atmospheric water vapor content, ρ(2) and ρ(19) represent the apparent reflectance of the 2nd and 19th bands of MODIS data, respectively; α and β are constants; C5) to invert the surface temperature. The Jiménez-Munoz single-channel algorithm was used to invert the surface temperature. Since this method was proposed, it has achieved good results in different research areas. When the atmospheric water vapor content in the air is 0.5 g cm -2 ~2g·cm -2 When the temperature is between 0.04 and 0.17°, the error of inverting the surface temperature using the single-channel method should be between 1K and 2K. C6) Surface temperature classification was carried out by using remote sensing images at different times. To increase comparability, the surface temperature inversion data of different time series were used more accurately to analyze the spatiotemporal changes of surface temperature in the study area. The surface temperature was classified by using the mean value and standard deviation. The temperature was classified by calculating the mean value and standard deviation of the entire study area, and then the temperature was classified according to the mean value ± standard deviation. By comparing the relative high and low temperatures within the region, the influence of regional factors on temperature can be eliminated, eliminating interference factors for subsequent research on the relationship between land cover change and temperature.

4. The method for studying land cover change and thermal environment impact based on multi-temporal data according to claim 1 is characterized by: Analyzing the temperature response of the surface temperature of different land cover types in different seasons in the step 3) specifically refers to using the previous land cover classification results, combining the surface temperature map of the study area at the same time, and using the ArcGIS spatial analysis module to overlay the surface temperature map with the land type map to obtain the maximum, minimum, average, and standard deviation of the surface temperature of the land cover types in different seasons in the monitored area and perform a statistical test on a single sample, the purpose of which is to test whether the average surface temperature between different land use types reaches a significant level, wherein the significance level is set to 0.05; In step 3), introducing the semi-quantitative study of the relationship between the surface characteristic parameters and the surface temperature specifically refers to introducing some specific indexes that can represent the surface information to explore the relationship between the index and the surface temperature; wherein, the Normalized Difference Vegetation Index (NDVI) is selected as a representative for studying vegetation, the Modified Normalized Difference Water Index (MNDWI) is selected as a representative for studying water areas, and the Normalized Building Index (NDBI) is selected as a representative for studying buildings; After partially eliminating the seasonal influence of the surface temperature by using the impervious surface information in the step 3), the contribution ratio of the land cover type change to the surface temperature is quantitatively analyzed. Specifically, in order to further quantify the urban land cover type and its relationship with the thermal environment, the impervious surface coverage is divided into categories representing different urban development levels; the pixels are grouped by using the equidistant classification method, and each pixel is divided into 10 equal levels (0-10%, 10-20%, 30%, ..., 90-100%) according to its abundance value; the land surface temperature within each impervious surface grouping range is grouped by this method, which is helpful for detailed analysis of the land surface temperature in different groups, and the spatial distribution of urban thermal environment in areas with different urban densities can be analyzed and compared; In step 3), the distribution function is introduced to quantitatively analyze the heat source area information of different land cover types in the study area, and the impact of land cover type changes on the heat source area is further explored through the trajectory of center of gravity movement. Specifically, the concept of distribution index (DI) is introduced to quantify the contribution of each land cover type to the thermal effect and determine its nature as a heat source or sink. Get, where S i represents the area of ​​the i-th land cover type; S ih represents the area of ​​high temperature (including high temperature and extremely high temperature) in the i-th land cover type; S represents the total area of ​​the study area; S h It indicates the area of ​​high temperature (including high temperature and extremely high temperature) surface temperature in the study area. If DI>1, it means that this land cover type has a great contribution to the thermal environment of the study area and should be defined as a heat source. If DI<1, it means that this land type contributes very little to the thermal environment of the study area.

Images