High-temporal-spatial-resolution surface temperature reconstruction method used in cloud and mist environment

Abstract

The invention belongs to the field of land surface temperature reconstruction, and relates to a high-spatial-temporal-resolution land surface temperature reconstruction method used in a cloud environment, which comprises the following steps: calculating a first reference land surface temperature and a first land surface temperature residual error of a clear sky pixel under a low spatial resolution according to a first land surface temperature annual change model; seamless surface reflectance data under high spatial resolution are obtained; defining a relative solar radiation index based on the digital elevation model data with high spatial resolution; screening clear sky pixels participating in modeling of the XGBoost model; an XGBoost model is adopted, and a plurality of model relations between a low-spatial-resolution land surface temperature annual change model coefficient and a land surface attribute and between a day-by-day land surface temperature residual error and the land surface attribute are constructed; obtaining a high-spatial-resolution surface temperature annual change model coefficient and a day-by-day surface temperature residual error based on the constructed multiple model relationships, and obtaining a high-spatial-resolution daily surface temperature; the accuracy and the stability of a surface temperature reconstruction result are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of land surface temperature reconstruction, and specifically discloses a high spatiotemporal resolution land surface temperature reconstruction method in a cloud and fog environment. BACKGROUND

[0002] Land surface temperature (LST) is a key parameter in the process of dynamic balance of land surface energy, and has important scientific research and application value in meteorology, hydrology, ecology and geology. At present, satellite observation has become the only way to obtain land surface temperature at the pixel scale due to its wide range of observation and regular revisit. However, it is difficult to meet the needs of high temporal resolution and high spatial resolution at the same time in the current satellite remote sensing of land surface temperature. Satellites with high temporal resolution usually have low spatial resolution, and satellites with high spatial resolution have long revisit periods, which makes it difficult to meet the needs of continuous and high-precision monitoring. In addition, the current remote sensing inversion of land surface temperature relies on the thermal infrared band, and thermal infrared remote sensing is sensitive to atmospheric conditions and is easily disturbed by clouds and fog. Especially in mountainous areas with complex terrain, the available land surface temperature observation data is extremely limited due to the influence of clouds and fog all year round, which further exacerbates the difficulty of land surface temperature application in mountainous areas.

[0003] At present, some studies have proposed cloud-free reconstruction and spatial downscaling methods for land surface temperature to alleviate the problem of data missing caused by cloud and fog obstruction or insufficient spatial resolution in remote sensing observation. However, most of these methods focus on solving a single problem and are difficult to fully cope with the complex situation of temporal and spatial information missing in areas severely disturbed by clouds and fog. In addition, the step-by-step processing strategy of "cloud-free reconstruction first, then spatial downscaling" or "downscaling first, then temporal interpolation" is prone to error accumulation during processing, which ultimately reduces the accuracy and stability of the land surface temperature reconstruction results.

[0004] Therefore, the application provides a high spatio-temporal resolution land surface temperature reconstruction method in a cloud and fog environment. SUMMARY

[0005] The application aims to provide a high spatio-temporal resolution land surface temperature reconstruction method in a cloud and fog environment, solve the problem, improve the accuracy and stability of the land surface temperature reconstruction result, and specifically provide the following technical scheme.

[0006] Further, the step 2 comprises: constructing a first ground surface temperature annual variation model, calculating a first reference ground surface temperature through the first ground surface temperature annual variation model; calculating a daily first ground surface temperature residual based on a difference between the daily first ground surface temperature and the first reference ground surface temperature; the first ground surface temperature refers to an actual ground surface temperature obtained according to low spatial resolution data.

[0007] Further, the step 3 comprises: filling the missing part in the low spatial resolution reflectivity data to be fused to obtain seamless low spatial resolution reflectivity filling data; fusing the seamless low spatial resolution filling data and the high spatial resolution data through a space-time fusion technology, filling the missing part of the high spatial resolution reflectivity data through the reflectivity obtained by the fusion, and obtaining high spatial resolution reflectivity filling data; taking the high spatial resolution reflectivity filling data as one of the ground surface attribute data under the high spatial resolution.

[0008] Further, the calculation formula for filling the missing data in the low spatial resolution data to be fused is: ; In the formula, represents the ground surface reflectivity of the missing pixel, represents the ground surface reflectivity of the missing pixel at the previous clear sky pixel, represents the ground surface reflectivity of the missing pixel at the next clear sky pixel, represents the annual accumulated day of the missing pixel in a year, and respectively represent the annual accumulated days of the previous and next clear sky pixels corresponding to the missing pixel.

[0009] Further, the step 4 comprises: calculating the solar radiation received by the top of the atmosphere based on the solar constant; determining the sunset and sunrise times based on the hour angle; determining the solar radiation received by the rugged mountain surface based on the solar radiation received by the top of the atmosphere, the sunset and sunrise times, and the solar incident angle under the mountain terrain; determining the solar radiation received by the flat surface based on the solar radiation received by the top of the atmosphere, the sunset and sunrise times, and the solar zenith angle; determining the relative solar radiation index based on the ratio of the solar radiation received by the rugged mountain surface to the solar radiation received by the flat surface.

[0010] Further, the calculation formula of the solar radiation received by the top of the atmosphere is: ; The calculation formula of the sunset and sunrise times is: ; ; ; The formula for calculating the solar radiation received by the rugged mountainous surface is: ; The formula for calculating the solar radiation received by the flat surface is: ; The formula for calculating the relative solar radiation index is: ; In the formula, represents the relative solar radiation index on the dth day of the year, and respectively represent the solar radiation received by the rugged mountainous surface and the flat surface; is the solar radiation received by the top of the atmosphere on the dth day of the year, cos represents the cosine function, dt represents the derivative with respect to time t, and t represents the time variable, is the solar zenith angle, is the solar incident angle under the mountainous terrain; and respectively represent the sunset and sunrise times of the sun; is the hour angle corresponding to the sunrise and sunset of the sun; is the latitude; is the solar declination; d represents the date variable; represents the circular constant, is the solar constant, and cos represents the cosine function.

[0011] Further, step 5 comprises: determining the number of pixel thresholds of the clear-sky pixels participating in the XGBoost model training; taking the geometric center of the study area as the center, searching for the pixel threshold number of clear-sky pixels from short to long in distance, obtaining a clear-sky pixel set for participating in the XGBoost model training; and taking the distance between the center and the pixel farthest from the clear-sky pixel set as the radius of the dynamic window.

[0012] Further, the step 6 comprises: establishing a first model relationship between the first land surface attribute and the first land surface temperature residual at a low spatial resolution by an XGBoost model; the first land surface attribute comprises latitude, elevation, slope, sky view factor, small-scale self-heating factor, solar incidence angle, normalized vegetation index, land surface reflectivity, land use type and cross-polarization ratio; establishing a second model relationship between the second land surface attribute and the annual average temperature parameter in the first land surface temperature annual variation model at a low spatial resolution by an XGBoost model; the second land surface attribute comprises latitude, elevation, slope, sky view factor, small-scale self-heating factor, relative solar index, annual average normalized vegetation index, annual average reflectivity and annual average cross-polarization ratio; establishing a third model relationship between the third land surface attribute and the annual amplitude parameter in the first land surface temperature annual variation model at a low spatial resolution by an XGBoost model; the third land surface attribute comprises latitude, elevation, slope, relative solar index, annual average normalized vegetation index, normalized vegetation index variation, annual average reflectivity, annual average cross-polarization ratio and cross-polarization ratio variation.

[0013] Further, the calculation formula of the sky view factor is: ; In the formula, SVF represents the sky view factor, Di represents the total number of divided directions, k represents the direction variable, represents the ground blocking angle in the kth direction; sin represents the sine function; The calculation formula of the small-scale self-heating factor is: ; In the formula, SSP represents the small-scale self-heating factor, i represents the pixel variable, m and n are the pixel sizes of the coarse resolution and fine resolution DEM data respectively, and the size of the coarse resolution m is greater than the size of the fine resolution n; represents the number of fine resolution pixels contained in the coarse resolution pixel, and respectively represent the slope of the coarse resolution pixel and the slope of the fine resolution pixel, and j represents the fine resolution pixel variable in the coarse resolution pixel; The calculation formula of the cross-polarization ratio is: ; In the formula, q is the cross-polarization ratio, and respectively represent the backscattering coefficients of vv polarization and vh polarization in the SAR data.

[0014] Further, step 7 comprises: determining a plurality of high spatial resolution land surface property data based on the cloud-free land surface reflectance data under high spatial resolution and the high resolution land surface property resampling data; processing the plurality of high spatial resolution land surface property data through a plurality of model relationships to obtain a high spatial resolution land surface temperature annual variation model coefficient and a second land surface temperature residual; constructing a second land surface temperature annual variation model based on the high spatial resolution land surface temperature annual variation model coefficient, and calculating a second reference land surface temperature; and determining a high spatial resolution daily land surface temperature based on the sum of the second reference land surface temperature and the second land surface temperature residual.

[0015] The present application has the following advantages and beneficial effects: The present application has the following advantages and beneficial effects:

[0016] The present application has the following advantages and beneficial effects:

[0017] The present application has the following advantages and beneficial effects: BRIEF DESCRIPTION OF DRAWINGS

[0018] Figure 1 A flowchart of a high spatiotemporal resolution land surface temperature reconstruction method suitable for a cloud and fog environment is provided for the embodiments of the present application; Figure 2 A high elevation map of a study area and a ground station site diagram for verification are provided for the present application; Figure 3a An image of parameter a in the ATC model in the study area under high spatial resolution is provided for the present application; Figure 3b An image of parameter b in the ATC model in the study area under high spatial resolution is provided for the present application; Figure 4a A comparison diagram of the reconstructed land surface temperature and the ground station temperature under all conditions is provided for the present application; Figure 4b A comparison chart of reconstructed surface temperature and ground station temperature for the case of cloud cover greater than 0.99; Figure 5 A graph of the average RMSE of XGBoost training under different dynamic window radii; Figure 6 A low-resolution surface temperature map of the study area before reconstruction provided by the present application is disturbed by clouds and fog; Figure 7 A high-resolution all-weather surface temperature map of the study area after reconstruction provided by the present application; Figure 8 A RSRI graph of the relative solar index of the study area provided by the present application. DETAILED DESCRIPTION

[0019] In order to make the purpose, technical scheme and advantages of the embodiments of the present application clearer, the technical scheme in the embodiments of the present application will be described clearly and completely below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. The components of the embodiments of the present application described and shown in the drawings herein can be arranged and designed in various different configurations.

[0020] Figure 1 A flowchart of the high-temporal and high-spatial resolution surface temperature reconstruction method suitable for cloud and fog environment provided by the embodiments of the present application is shown. As shown in Figure 1 The present application discloses a high-temporal and high-spatial resolution surface temperature reconstruction method for cloud and fog environment, which includes the following contents: Step 1: Collect remote sensing image data; the remote sensing image data includes low spatial resolution data and high spatial resolution data.

[0021] Step 2: Calculate the first reference surface temperature and the first surface temperature residual of the clear sky pixel under low spatial resolution according to the first surface temperature annual variation model; the first surface temperature annual variation model refers to the surface temperature annual variation model (Annual Temperature Cycle, ATC) constructed by the low spatial resolution data. The first reference surface temperature refers to the surface temperature output by the first surface temperature annual variation model under low spatial resolution. The first surface temperature residual refers to the surface temperature residual under low spatial resolution. For example, the difference between the surface temperature fitted by the first surface temperature annual variation model and the MOD11A1 surface temperature.

[0022] In some embodiments, step 2 comprises: constructing a first land surface temperature annual variation model, and calculating a first reference land surface temperature by the first land surface temperature annual variation model; and calculating a first land surface temperature residual based on a difference between the daily first land surface temperature and the first reference land surface temperature; the first land surface temperature refers to an actual land surface temperature obtained according to low spatial resolution data.

[0023] Step 3: based on an ESTARFM (Enhanced STARFM, Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model) spatio-temporal fusion model, fusing low spatial resolution reflectivity data and high spatial resolution reflectivity data to obtain seamless land surface reflectivity data at high spatial resolution.

[0024] In some embodiments, step 3 comprises: filling in missing data in the low spatial resolution data to be fused to obtain low spatial resolution filled-in data; and fusing the low spatial resolution filled-in data and the high spatial resolution data by a spatio-temporal fusion technique, filling in missing parts of the high spatial resolution data by the reflectivity data obtained by the fusion to obtain high spatial resolution filled-in data; the low spatial resolution filled-in data refers to complete low spatial resolution data obtained by filling in missing data in the low spatial resolution data; the high spatial resolution filled-in data refers to seamless high spatial resolution data obtained by filling in original high spatial resolution data using the synthesized high resolution data obtained by the ESTARFM; and the high spatial resolution reflectivity filled-in data is used as one of the land surface attribute data at high spatial resolution.

[0025] Step 4: defining a relative solar radiation index (RSRI) based on high spatial resolution digital elevation model (DEM) data; the relative solar radiation index is used to evaluate the difference between the solar radiation received by a mountainous land surface and the solar radiation received by a flat land surface.

[0026] In some embodiments, step 4 comprises: calculating solar radiation received by the top of the atmosphere based on a solar constant; determining sunset and sunrise times based on a time angle; determining solar radiation received by a rugged mountainous land surface based on the solar radiation received by the top of the atmosphere, the sunset and sunrise times, and a solar incident angle under the mountainous terrain; determining solar radiation received by a flat land surface based on the solar radiation received by the top of the atmosphere, the sunset and sunrise times, and a solar zenith angle; and determining the relative solar radiation index based on a ratio of the solar radiation received by the rugged mountainous land surface to the solar radiation received by the flat land surface.

[0027] Step 5: Select clear-sky pixels for XGBoost model modeling using a dynamic window centered on the study area. The dynamic window can be sorted by the distance of each pixel from the center of the study area. The distance can refer to the distance between the geometric center of the study area and the geometric center of the pixel. For example, reconstruct the land surface temperature based on the MOD11A1 1 km daily land surface temperature product, and select clear-sky pixels according to the quality control flag of MOD11A1. XGBoost models the two parameters related to annual average temperature and annual amplitude in the ATC model and the residual in the ATC fitting process, respectively. Because the spatial distribution of clear-sky pixels selected is different due to different cloud cover each day, the radius of the dynamic window is set dynamically according to the cloud cover each day.

[0028] In some embodiments, step 5 includes: determining a pixel threshold number of clear-sky pixels participating in XGBoost model training; centering on the geometric center of the study area, searching for the pixel threshold number of training clear-sky pixels from short to long distance to obtain a clear-sky pixel set; and taking the distance between the center and the farthest pixel in the clear-sky pixel set as the radius of the dynamic window. The pixel threshold number refers to the total number of clear-sky pixels participating in modeling. The training clear-sky pixel refers to the clear-sky pixel used to participate in XGBoost model training. For example, the available clear-sky pixels are searched step by step outward according to the distance of the pixel from the geometric center of the study area, and the distance between the current farthest pixel from the center and the center is taken as the radius of the dynamic window.

[0029] Step 6: Use the XGBoost model (extreme gradient boosting model) to build a low spatial resolution land surface temperature annual variation model coefficient and a plurality of model relationships between daily land surface temperature residuals and land surface attributes. Use XGBoost as a regression model to establish the nonlinear relationship between the two parameters of the ATC model and the land surface attributes and the LST residuals and the land surface attributes. The two parameters of the ATC model include a and b, which represent the parameters related to the annual average LST and the LST amplitude, respectively.

[0030] In some embodiments, step 6 includes: establishing a first model relationship between a first land surface attribute and a first land surface temperature residual at a low spatial resolution by an XGBoost model; the first land surface attribute includes latitude (Lat), elevation, slope (Slope), sky view factor (SVF), small-scale self-heating factor (SSP), solar incident angle (SEI), normalized vegetation index (NDVI), land surface reflectance (SR, including red, green, blue, and near-infrared bands), land use type, and cross-polarization ratio (q). determining a second model relationship between the second land surface properties and the annual mean temperature parameter in the first land surface temperature annual variation model at a low spatial resolution by the XGBoost model; the second land surface properties include latitude (Lat), elevation (Elv), slope (Slope), sky view factor (SVF), small-scale self-heating factor (SSP), relative sun index (RSRI), annual mean normalized difference vegetation index (aNDVI), annual mean reflectance (aSR, including red band, green band, blue band and near-infrared band) and annual mean cross-polarization ratio (aq).

[0031] determining a third model relationship between the third land surface properties and the annual amplitude parameter in the first land surface temperature annual variation model at a low spatial resolution by the XGBoost model; the third land surface properties include latitude (Lat), elevation (Elv), slope (Slope), relative sun index (RSRI), annual mean normalized difference vegetation index (aNDVI), normalized difference vegetation index variation (vNDVI), annual mean reflectance (aSR, including red band, green band, blue band and near-infrared band), annual mean cross-polarization ratio (aq) and cross-polarization ratio variation (vq).

[0032] Step 7: applying the plurality of model relationships to the high spatial resolution land surface property data respectively to obtain the high spatial resolution land surface temperature annual variation model coefficients and the daily land surface temperature residual, and then obtaining the daily land surface temperature at a high spatial resolution in the study area. For example, by obtaining the high spatial resolution ATC model coefficients, a second land surface temperature annual variation model is constructed, and the daily high spatial resolution LST is obtained by adding the daily high spatial resolution LST residual to the reference value of the daily high spatial resolution LST calculated by the second land surface temperature annual variation model. The second land surface temperature annual variation model is used to obtain the reference land surface temperature at a high spatial resolution (i.e., the second reference land surface temperature).

[0033] In some embodiments, step 7 includes: determining a plurality of high spatial resolution land surface property data based on the cloud-free land surface reflectance data at a high spatial resolution and the high resolution land surface property resampling data; processing the plurality of high spatial resolution land surface property data by the plurality of model relationships to obtain the high spatial resolution land surface temperature annual variation model coefficients and the second land surface temperature residual; constructing a second land surface temperature annual variation model based on the high spatial resolution land surface temperature annual variation model coefficients and calculating the second reference land surface temperature; and determining the daily land surface temperature at a high spatial resolution based on the sum of the second reference land surface temperature and the second land surface temperature residual. The second land surface temperature residual refers to the high spatial resolution land surface temperature residual. The second land surface temperature annual variation model is used to calculate the reference land surface temperature at a high spatial resolution. The second reference land surface temperature refers to the reference land surface temperature at a high spatial resolution.

[0034] Embodiment 1 As an example, as shown in FIG. 1, taking the Wanglang Nature Reserve in Sichuan as an example, the embodiment provides a high-spatial-temporal resolution land surface temperature reconstruction method in cloud and fog environment, including construction of a land surface temperature annual variation model, acquisition of high-spatial-temporal resolution cloud-free land surface reflectance data, definition of relative solar radiation index, dynamic window screening of XGBoost model input pixels, model training, and acquisition of high-spatial resolution cloud-free land surface temperature. The research area and ground stations for verification are as shown in FIG. 2. Figure 1 Figure 2

[0035] Step 1: Collect remote sensing image data, including high-resolution remote sensing images and low-resolution remote sensing images. Low-spatial resolution data includes thermal infrared data (MOD11A1), land surface reflectance data (MOD09A1), and normalized vegetation index NDVI data (MOD13A2), which are resampled to 1 km. High-spatial resolution data includes land surface reflectance data (Landsat-8 / 9, Sentinel-2), DEM data (Alos), SAR data (Sentinel-1), and land classification data (CLCD), which are resampled to 30 m. Then resample DEM, SAR, and land classification data to 1 km to meet the needs of building a model at low spatial resolution.

[0036] Step 2: According to the Annual Temperature Cycle (ATC) model, calculate the reference LST and residual of clear-sky pixels at low spatial resolution per day, where the daily land surface temperature can be represented as: ; In the formula, T(d) is the land surface temperature on the dth day of the year, is the LST reference value calculated by the ATC model, is the land surface temperature residual. The ATC model can be represented as: ; In the formula, represents the annual mean temperature parameter, represents the annual amplitude parameter, is the phase shift relative to the equinox, N is the length of the annual cycle, usually set to 365 days, ignoring leap day, e represents the exponential function; NDVI(d) represents the normalized vegetation index on the dth day of the year; sin represents the sine function; d represents the date variable; represents the circular constant. For a given pixel, the parameters a and b are fitted by the least squares algorithm using existing clear-sky LST observation data and NDVI data to construct the annual cycle model.​​

[0037] For this embodiment, the MOD11A1 land surface temperature product is used, the clear-sky pixel observation data is screened through the quality control layer, and the MOD13A2 NDVI product is used to obtain the NDVI corresponding to the clear-sky observation LST.

[0038] Step 3: Considering that the study area is affected by clouds and fog, there are serious missing data in the high-resolution land surface reflectance data, the ESTARFM spatio-temporal fusion technology is used to fill the missing part of the high-resolution land surface reflectance data. The high-resolution land surface reflectance data uses Landsat-8 / 9 and Sentinel-2 data, and the low-resolution land surface reflectance data used for fusion uses the MOD09A1 8-day composite reflectance product. For the areas seriously affected by cloud and fog, the missing data of the 8-day composite reflectance product is still missing, and the missing data is filled before spatio-temporal fusion: ; In the formula, represents the land surface reflectance of the missing pixel, represents the land surface reflectance of the previous clear-sky pixel of the missing pixel in the previous period, represents the land surface reflectance of the next clear-sky pixel of the missing pixel in the next period; represents the annual accumulated day of the missing pixel in a year, and respectively represent the annual accumulated day of the clear-sky pixel of the previous image and the next image of the missing pixel.

[0039] For the ESTARFM spatio-temporal fusion method, the low spatial resolution reflectance data corresponding to the target time and each group of corresponding high spatial resolution and low spatial resolution land surface reflectance images before and after the target time are needed. For the missing part of the Landsat-8 / 9 reflectance data, the required land surface reflectance data is screened out for spatio-temporal fusion, and the fused reflectance data is used to fill the missing part of the original data to obtain high spatial resolution cloud-free reflectance data every 8 days.

[0040] Step 4: In order to measure the influence of mountain terrain on land surface temperature, the relative solar radiation index (RSRI) is defined to evaluate the difference between the solar radiation received by the mountain surface and the solar radiation received by the flat surface: ; In the formula, represents the relative solar radiation index on the dth day in a year, and respectively represent the solar radiation received by the mountainous rugged surface and the flat surface, ignoring the influence of atmospheric transmittance, wherein: ; ; wherein, is the solar radiation received by the top of the atmosphere on the dth day of the year, cos represents the cosine function, dt represents the derivative with respect to time t, t represents the time variable, is the solar zenith angle, is the solar incident angle under the complex terrain of the mountainous area; and respectively represent the sunset and sunrise times of the sun.

[0041] ; ; ; ; ; ; ; ; ; wherein, is the hour angle, tan represents the tangent function, represents the latitude, represents the solar declination, sin represents the sine function, and d represents the date variable; represents the circular constant, is the solar constant, cos represents the cosine function, represents the solar azimuth angle, slope and Aspect respectively represent the slope and aspect of the ground. Thus, the relative solar radiation index RSRI of each pixel of the study area in a year can be represented as: ; In this embodiment, the relative solar index is calculated by ALOS digital elevation data, Figure 8 is the RSRI spatial distribution image of the study area calculated in this embodiment.

[0042] Step 5: Taking the geometric center of the study area as the center, the available clear-sky pixels are searched step by step outward according to the distance of the pixel from the study area, and the clear-sky pixels available for participating in the XGBoost model training reach a pre-set threshold number. The distance of the pixel farthest from the center of the study area is the radius of the dynamic window.

[0043] Step 6: Using the XGBoost model, the model relationship between the low spatial resolution ATC model coefficients (annual mean temperature and annual amplitude), daily LST residual and land surface attributes in step 1 was constructed respectively.

[0044] For daily LST residual, first, the relationship between latitude (Lat), elevation (Elv), slope (Slope), sky view factor (SVF), small-scale self-heating factor (SSP), solar incident angle (SEI), normalized difference vegetation index (NDVI), land surface reflectance (SR, including red, green, blue and near-infrared bands), land use type, cross-polarization ratio (q) and residual was established at 1 km resolution.

[0045] ; In the formula, SVF represents the sky view factor, Di indicates that Di directions are divided in total, k represents the direction variable, represents the ground blocking angle in the kth direction; sin represents the sine function.

[0046] ; In the formula, SSP represents the small-scale self-heating factor, i represents the pixel variable, and m and n are the pixel sizes of the coarse resolution and fine resolution DEM data, respectively, and the size of the coarse resolution m is larger than that of the fine resolution n; represents the number of fine resolution pixels contained in the coarse resolution pixel, which is a positive integer, and the value is represents the number of fine resolution pixels contained in a coarse resolution pixel, and respectively represent the slope of the coarse resolution pixel and the slope of the fine resolution pixel, and j represents the fine resolution pixel variable in the coarse resolution pixel. For the ALOS DEM data used in this embodiment, the 12.5 m resolution DEM data is resampled to 500 m and 1 km respectively as the fine resolution pixel and coarse resolution pixel in the SSP calculation process.

[0047] ; In the formula, q is the cross-polarization ratio, which is one of the parameters containing soil moisture information, and respectively represent the backscattering coefficients of vv and vh polarization in SAR data.

[0048] The DEM data in this example comes from ALOS 12.5 m DEM data, the surface reflectance data comes from MOD09A1 surface reflectance product, the NDVI comes from MOD13A2 NDVI product, the land use data comes from 30 m CLCD data, and the SAR data comes from Sentinel-1. The spatial resolutions of the above data are unified to 1 km through resampling, and the model relationship between these surface attributes and LST residual is established at 1 km spatial resolution through XGBoost.

[0049] For the ATC model parameter a, the relationship between latitude (Lat), elevation (Elv), slope (Slope), sky view factor (SVF), small-scale self-heating factor (SSP), relative solar index (RSRI), annual average normalized vegetation index (aNDVI), annual average reflectivity (aSR, including red, green, blue, and near-infrared bands), and a is established at 1 km resolution. The annual average normalized vegetation index, the annual average reflectivity, and the average cross-polarization ratio are directly obtained by averaging the cross-polarization ratio data calculated from MOD13A2, MOD09A1, and Sentinel-1.

[0050] For the ATC model parameter β, the relationship between latitude (Lat), elevation (Elv), slope (Slope), relative solar index (RSRI), annual average normalized vegetation index (aNDVI), normalized vegetation index variation (vNDVI), annual average reflectivity (aSR, including red, green, blue, and near-infrared bands), annual average cross-polarization ratio (aq), and β is established at 1 km resolution. The normalized vegetation index and the average cross-polarization ratio are directly obtained by calculating the difference between the maximum and minimum values of the cross-polarization ratio data calculated from MOD13A2 and Sentinel-1.

[0051] The model relationship between the ATC model coefficients, the daily LST residual, and the surface attributes is established respectively, and the established model is used for the corresponding high-resolution surface attribute data to obtain the high-resolution ATC model coefficients and the daily LST residual. In the surface attribute data, the reflectivity data and the NDVI data come from the results of ESTARFM data fusion, and the remaining surface attribute data comes from the resampling of high-resolution surface attribute data. Figure 3a and Figure 3b The spatial distribution of high-resolution ATC model parameters a and β in the study area in 2022 is shown respectively.

[0052] Step 6: Based on the high-resolution ATC model coefficients obtained in step 5, calculate the daily LST reference value in the study area, and obtain the high spatial resolution daily LST image in the study area by combining the daily LST residual.

[0053] By Figure 4a As can be seen, by comparing with the measurement data of the three ground stations of Wanglang Station, the RMSE of the reconstructed high-resolution land surface temperature is 3.78 K. By Figure 4b As can be seen, on the days with serious cloud cover (cloud cover greater than 0.99), by comparing with the measurement data of the three ground stations of Wanglang Station, the RMSE of the reconstructed high-resolution land surface temperature is 3.92 K, and the reconstruction accuracy has a small amount of decline compared with the annual reconstruction result. Figure 5 The average RMSE of using the XGBoost model to establish the relationship between the daily LST residual and the surface attribute under different dynamic window radii is shown. When the dynamic window radius is less than 100 km, the average RMSE of the XGBoost model and the RMSE between the measured temperature of the ground station and the reconstructed land surface temperature are the minimum values. As the cloud and fog interference gradually increases, the radius of the dynamic window also increases, and the average RMSE of the XGBoost model and the RMSE between the measured temperature of the ground station and the reconstructed land surface temperature also increase accordingly. The Figure 6 The MOD11A1 land surface temperature product image on the 15th of each month in 2022 is shown. Before reconstruction, the study area is seriously interfered by cloud and fog, and there is a large amount of data missing. The Figure 7 The land surface temperature after reconstruction on the 15th of each month in 2022 is shown. It can be found that the spatial resolution of the reconstructed land surface temperature has been significantly improved, and the part affected by cloud and fog interference and resulting in data missing has also been effectively filled.

[0054] The above is only a preferred embodiment of the present application and is not used to limit the present application. For those skilled in the art, the present application can have various modifications and changes. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A method for reconstructing land surface temperature with high temporal and spatial resolution in a foggy environment, characterized in that: include: Step 1: Collect remote sensing image data; remote sensing image data includes low spatial resolution data and high spatial resolution data; Step 2: Calculate the first reference surface temperature and the first surface temperature residual for clear-sky pixels at low spatial resolution every day based on the first surface temperature annual variation model; the first surface temperature annual variation model refers to the surface temperature annual variation model constructed using low spatial resolution data; Step 3: Based on the spatiotemporal fusion model, the low spatial resolution surface reflectance data and the high spatial resolution surface reflectance data are fused to obtain seamless surface reflectance data at high spatial resolution; Step 4: Define the relative solar radiation index based on the high spatial resolution digital elevation model data; The relative solar radiation index is used to evaluate the difference in solar radiation received by mountainous surfaces and flat surfaces; Step 5: Centering the study area, use a dynamic window to filter clear sky pixels for XGBoost modeling. Step 6: Use the XGBoost model to construct low spatial resolution surface temperature annual variation model coefficients and multiple model relationships between daily surface temperature residuals and surface attributes; Step 7: Apply multiple model relationships to high-spatial-resolution surface attribute data to obtain high-spatial-resolution surface temperature annual variation model coefficients and daily surface temperature residuals, and obtain high-spatial-resolution daily surface temperature in the study area.

2. The high temporal and spatial resolution surface temperature reconstruction method for foggy environment according to claim 1 is characterized in that: Step 2 includes: Constructing a first surface temperature annual variation model, and calculating a first reference surface temperature through the first surface temperature annual variation model; The daily first surface temperature residual is calculated based on the difference between the daily first surface temperature and the first reference surface temperature; the first surface temperature refers to the actual surface temperature obtained based on low spatial resolution data.

3. The method for reconstructing land surface temperature with high temporal and spatial resolution in a foggy environment according to claim 1, characterized in that: Step 3 includes: Filling the missing parts in the low spatial resolution reflectance data to be fused to obtain seamless low spatial resolution reflectance filling data; The seamless low spatial resolution filling data is fused with the high spatial resolution data through the spatiotemporal fusion technology, and the missing parts of the high spatial resolution reflectivity data are filled with the reflectivity obtained by fusion to obtain the high spatial resolution reflectivity filling data; The high spatial resolution reflectance filling data is used as one of the surface attribute data at high spatial resolution.

4. The high temporal and spatial resolution surface temperature reconstruction method for foggy environment according to claim 3 is characterized in that: The calculation formula for filling the missing data in the low spatial resolution data to be fused is: ； Where, represents the surface reflectance of the missing pixel, Represents the surface reflectance of the missing pixel in the previous clear sky pixel, Represents the surface reflectance of the missing pixel in the next clear sky pixel, represents the cumulative number of missing pixels in a year, and They represent the annual days corresponding to the clear sky pixel before and after the missing pixel respectively.

5. The method for reconstructing land surface temperature with high temporal and spatial resolution in a foggy environment according to claim 1, characterized in that: Step 4 includes: Calculate the solar radiation received by the top of the atmosphere based on the solar constant; Determine the time of sunset and sunrise based on the hour angle; Determine the solar radiation received by the mountainous rugged surface based on the solar radiation received at the top of the atmosphere, the time of sunset and sunrise, and the angle of incidence of the sun in mountainous terrain; Determine the solar radiation received by a flat surface based on the solar radiation received at the top of the atmosphere, the times of sunset and sunrise, and the solar zenith angle; The relative solar radiation index is determined based on the ratio of solar radiation received by the rugged surface of the mountain to that received by the flat surface.

6. The method for reconstructing land surface temperature with high temporal and spatial resolution in a foggy environment according to claim 5, characterized in that: The calculation formula for the solar radiation received by the top of the atmosphere is: ； The formula for calculating the time of sunset and sunrise is: ； ； ； The calculation formula for the solar radiation received by the rugged surface of the mountain is: ； The calculation formula for the solar radiation received by a flat surface is: ； The calculation formula for the relative solar radiation index is: ; Where, represents the relative solar radiation index on the dth day of the year, and They represent the solar radiation received by the rugged surface of the mountain and the flat surface respectively; is the solar radiation received by the top of the atmosphere on the dth day of the year, cos represents the cosine function, dt represents the derivative with respect to time t, and t represents the time variable. is the solar zenith angle, is the solar incidence angle in mountainous terrain; and Represents the time of sunset and sunrise respectively; The hour angle corresponding to the sunrise and sunset of the sun; is latitude; is the solar declination; d represents the date variable; represents pi, is the solar constant, and cos represents the cosine function.

7. The method for reconstructing land surface temperature with high temporal and spatial resolution in a foggy environment according to claim 1, characterized in that: Step 5 includes: Determine the pixel threshold number of clear sky pixels participating in XGBoost model training; Taking the geometric center of the study area as the center, the clear sky pixels with a pixel threshold number are gradually searched from short to long distances, and the clear sky pixel set is obtained for participating in the XGBoost model training; The distance between the center and the farthest pixel in the clear sky pixel set is used as the radius of the dynamic window.

8. The method for reconstructing land surface temperature with high temporal and spatial resolution in a foggy environment according to claim 1, characterized in that: Step 6 includes: establishing a first model relationship between a first surface attribute and a first surface temperature residual at low spatial resolution using an XGBoost model; the first surface attribute includes latitude, elevation, slope, sky visibility factor, small-scale self-heating factor, solar incidence angle, normalized difference vegetation index, surface reflectance, land use type, and cross-polarization ratio; A second model relationship between a second surface attribute and the annual average temperature parameter in the first surface temperature annual variation model is established at low spatial resolution using the XGBoost model; the second surface attribute includes latitude, elevation, slope, sky visibility factor, small-scale self-heating factor, relative solar index, annual average normalized difference vegetation index, annual average reflectance, and annual average cross-polarization ratio; The third model relationship between the third surface attribute and the annual amplitude parameter in the first surface temperature annual change model is established at low spatial resolution through the XGBoost model; the third surface attribute includes latitude, elevation, slope, relative solar index, annual average normalized vegetation index, normalized vegetation index change, annual average reflectance, annual average cross-polarization ratio and cross-polarization ratio change.

9. The method for reconstructing land surface temperature with high temporal and spatial resolution in a foggy environment according to claim 8, characterized in that: The calculation formula of sky visibility factor is: ； Where SVF represents the sky visibility factor, Di represents the total number of divided directions, and k represents the direction variable. represents the horizon occlusion angle in the 𝑘th direction; sin represents the sine function; The calculation formula of the small-scale self-heating factor is: ； Where SSP represents the small-scale self-heating factor, i represents the pixel variable, m and n are the pixel sizes of the coarse-resolution and fine-resolution DEM data, respectively. The size of the coarse-resolution m is larger than that of the fine-resolution n. Indicates the number of fine-resolution pixels contained in a coarse-resolution pixel, and They represent the slope of the coarse-resolution pixel and the slope of the fine-resolution pixel respectively, and j represents the fine-resolution pixel variable in the coarse-resolution pixel; The calculation formula of the cross-polarization ratio is: ； Where q is the cross-polarization ratio, and represent the backscattering coefficients of vv polarization and vh polarization in SAR data respectively.

10. The method for reconstructing land surface temperature with high temporal and spatial resolution in a foggy environment according to claim 1, characterized in that: Step 7 includes: Determine a variety of high spatial resolution surface attribute data based on cloud-free surface reflectance data at high spatial resolution and high-resolution surface attribute resampling data; Through multiple model relationships, a variety of high spatial resolution surface attribute data are processed to obtain high spatial resolution surface temperature annual variation model coefficients and second surface temperature residuals; A second surface temperature annual variation model is constructed based on the high spatial resolution surface temperature annual variation model coefficients, and a second reference surface temperature is calculated; A daily land surface temperature with high spatial resolution is determined based on the sum of the second reference land surface temperature and the second land surface temperature residual.

Images