CN114152350B - Earth surface temperature inversion method considering urban three-dimensional geometrical structure influence - Google Patents
Earth surface temperature inversion method considering urban three-dimensional geometrical structure influence Download PDFInfo
- Publication number
- CN114152350B CN114152350B CN202111494964.2A CN202111494964A CN114152350B CN 114152350 B CN114152350 B CN 114152350B CN 202111494964 A CN202111494964 A CN 202111494964A CN 114152350 B CN114152350 B CN 114152350B
- Authority
- CN
- China
- Prior art keywords
- pixel
- atmospheric
- surface temperature
- urban
- emissivity
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 40
- 230000005855 radiation Effects 0.000 claims abstract description 75
- 238000002834 transmittance Methods 0.000 claims abstract description 15
- 239000004576 sand Substances 0.000 claims abstract description 4
- 239000000463 material Substances 0.000 claims description 21
- 230000005540 biological transmission Effects 0.000 claims description 17
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 16
- 238000004364 calculation method Methods 0.000 claims description 14
- 230000000694 effects Effects 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 9
- 230000035699 permeability Effects 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims 1
- 238000000179 transient infrared spectroscopy Methods 0.000 claims 1
- 239000004567 concrete Substances 0.000 description 6
- 239000011521 glass Substances 0.000 description 4
- 230000003595 spectral effect Effects 0.000 description 4
- 239000004568 cement Substances 0.000 description 3
- 241000132092 Aster Species 0.000 description 2
- 238000000295 emission spectrum Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012950 reanalysis Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000005316 response function Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 244000025254 Cannabis sativa Species 0.000 description 1
- 241000196324 Embryophyta Species 0.000 description 1
- 239000010426 asphalt Substances 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
- G01J5/007—Radiation pyrometry, e.g. infrared or optical thermometry for earth observation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
- G01J2005/0077—Imaging
Abstract
The invention discloses a surface temperature inversion method considering the influence of an urban three-dimensional geometrical structure, which comprises the following steps of: step 1: obtaining flat earth surface temperature T by pixel-by-pixel inversion of formula (10)s_initialTaking the surface temperature at the moment as an initial value of the surface temperature of the adjacent pixel element; and 2, step: calculating the obtained sky visible factor, the obtained pixel effective emissivity of the urban ground surface, the obtained atmospheric transmittance, the atmospheric uplink radiation and the atmospheric downlink radiation, and the obtained radiance of the Landsat8 channel 10, simultaneously inputting a formula (9), and performing pixel-by-pixel inversion to obtain the urban ground surface temperature TsAnd taking the earth surface temperature at the moment as the earth surface temperature of the adjacent pixel for iterative computation of the radiance of the adjacent pixel; and step 3: and repeating the step 2 for iteratively calculating the surface temperature of the adjacent pixel element. The method has obvious advantages on the urban ground surface with a complex three-dimensional geometric structure, and can make up for larger errors caused by inverting the ground surface temperature in the prior art.
Description
Technical Field
The invention relates to a surface temperature inversion method considering the influence of an urban three-dimensional geometric structure, and belongs to the technical field of quantitative remote sensing.
Background
The surface temperature is a key factor of the interaction between the surface and the atmosphere and is also an important parameter for studying the surface energy balance. The knowledge of the space-time distribution of the earth surface temperature not only plays an important role in the earth surface energy process of regional and global scales, but also has important significance on urban climate research and urban thermal environment monitoring. The satellite remote sensing technology provides possibility for acquiring earth surface temperature with high space-time resolution in a large range. At present, scholars at home and abroad develop various satellite data surface temperature remote sensing inversion methods. However, most of these methods are directed to flat surfaces, and very little research is done on surfaces with complex geometries. Urban surface temperature is closely related to human activities, and urban land cover types have strong heterogeneity. The complex morphology of different types of buildings, artificial terrain and vegetation causes urban surface temperatures to exhibit strong spatial heterogeneity.
Under the influence of urban complex three-dimensional geometry, the inversion of urban surface temperature faces three challenges: (1) influence of multiple scattering and reflection inside the pixel on the emissivity of the urban pixel; (2) influence of urban three-dimensional geometrical structures on downward radiation of the atmosphere; (3) the neighboring picture elements radiate the effect on the target picture element. With the spatial resolution of the satellite sensor being higher and higher, the influence of the complex three-dimensional geometrical structure of the city on the inversion of the high-spatial-resolution thermal infrared data on the earth surface temperature is more and more prominent.
Most of the existing earth surface temperature inversion methods are directed at flat earth surfaces, are mainly suitable for natural earth surfaces such as bare soil, vegetation and the like, and do not consider the influence of complex underlying surfaces and earth surface geometrical structures on earth surface temperature inversion. The complex three-dimensional geometrical structure of the city causes the heat radiation transmission mode between the earth surface and the atmosphere to change, and influences the emissivity of an image element, the downlink radiation of the atmosphere, the radiation of an adjacent image element and the like, so that the earth surface temperature inversion method based on the conventional heat radiation transmission equation cannot be applied to the city earth surface with the complex three-dimensional geometrical structure, and the inversion accuracy of the temperature of the city earth surface is low. Accordingly, the prior art is deficient and needs improvement.
Disclosure of Invention
The invention aims to solve the technical problem of providing a surface temperature inversion method considering the influence of an urban three-dimensional geometric structure aiming at the defects of the prior art.
The technical scheme of the invention is as follows:
a surface temperature inversion method considering the influence of an urban three-dimensional geometrical structure comprises the following steps: the urban surface temperature inversion process comprises the following steps:
step 1: and (3) performing pixel-by-pixel inversion by using the obtained radiance of the Landsat8 channel 10, the calculated surface emissivity, the calculated atmospheric transmittance, the atmospheric uplink radiation and the atmospheric downlink radiation to obtain the flat surface temperature Ts_initialTaking the earth surface temperature at the moment as an initial value of the earth surface temperature of the adjacent pixel;
equation (10):
in the formula, Ts_initialObtaining flat surface temperature for inversion;
step 2: inputting the initial value of the surface temperature of the adjacent pixel into a formula (8) to obtain the radiance of the adjacent pixel, calculating the obtained sky visible factor, the obtained effective emissivity of the pixel of the urban surface, the obtained atmospheric transmittance, atmospheric uplink radiation and atmospheric downlink radiation and the obtained radiance of the Landsat8 channel 10, simultaneously inputting a formula (9), and performing pixel-by-pixel inversion to obtain the urban surface temperature TsAnd taking the earth surface temperature at the moment as the earth surface temperature of the adjacent pixel for iterative computation of the radiance of the adjacent pixel;
formula (8):
in the formula, subscripts A and B represent a target pixel and an adjacent pixel, respectively, LBIs the radiance of the adjacent pixel, theta is the angle between the connecting lines of the center points of the target pixel and the adjacent pixel, dSBIs the area of the adjacent picture element, rABThe distance between the target pixel and the adjacent pixel is defined, and N is the total number of the adjacent pixels;
formula (9):
in the formula, TsSurface temperature, B, for inversion taking into account the effects of urban three-dimensional geometry-1Is the inverse operation of the Planck function; l isTOAThe radiance received at the top of the atmosphere, LdFor atmospheric downlink radiation, LuFor upward radiation of the atmosphere, eeuIs the effective emissivity of the pixel on the surface of the city, tau is the atmospheric permeability, upsilonoutIs a sky visible factor L between Landsat8 pixelsadjRadiating for adjacent picture elements;
and step 3: and (3) repeating the step (2) for iteratively calculating the surface temperature of the adjacent pixel element until the difference between the surface temperatures of the adjacent pixel elements in the previous and next times is smaller than a certain value or the iteration times is larger than a certain value, and outputting a result to obtain the final urban surface temperature.
According to the surface temperature inversion method, the radiance of the Landsat8 channel 10 is obtained by the following steps: downloading on-satellite radiance data of a TIRS sensor carried on a Landsat8 satellite; data download website: https:// earth xplor. usgs. gov/; converting the on-satellite radiance data into an on-satellite radiance value by using a formula (1):
L＝gain×DN+offset (1)
in the formula, L is the on-satellite radiance, DN is the count value, gain and offset are gain and offset, respectively. For Landsat8 channel 10, the gain and offset were 0.00033420 and 0.1, respectively.
The earth surface temperature inversion method, the calculation method of the sky visible factor is as follows: the sky visible factor represents the visible degree of the earth's surface to the sky in the hemisphere space, that is, the ratio of the atmospheric downlink radiation received by the point to the atmospheric downlink radiation from the hemisphere space received when the point is flat and not shielded, the value range is 0 to 1, and the sky visible factor is calculated by the formula (2):
in the formula, gammaiIs the vertical elevation angle in a given search direction and n is the number of search directions.
The earth surface temperature inversion method is characterized in that the azimuth angle of the whole hemispherical space is dispersed into 16 directions which are uniformly distributed to calculate sky visible factors, wherein the sky visible factors are respectively 0 degree, 30 degrees, 45 degrees, 60 degrees, 90 degrees, 120 degrees, 135 degrees, 150 degrees, 180 degrees, 210 degrees, 225 degrees, 240 degrees, 270 degrees, 300 degrees, 315 degrees and 330 degrees.
The surface temperature inversion method comprises the following steps of:
(1) when the urban three-dimensional geometrical structure is not considered, calculating the pixel emissivity of the flat earth surface according to the weighted average value of the material emissivity of different ground components:
in the formula, epsilonefEmissivity of pixels for flat ground surface, fkIs the area ratio of the kth ground feature component in the pixel, εkIs the material emissivity of the kth terrestrial object component in the pixel;
(2) for urban earth surface with complex three-dimensional geometric structure, pixel effective emissivity is higher than pixel emissivity epsilon of flat earth surfaceefAnd calculating emissivity increment caused by multiple scattering and reflection among different surface object components in the pixel:
in the formula, epsiloneuEffective emissivity of pixel, delta epsilon, for urban surfacemThe emissivity increment caused by multiple scattering and reflection among different surface feature components in the pixel is shown, and N is the total number of adjacent pixels;
and (3) considering multiple scattering and reflection among different surface material components in the pixel, and finally obtaining the effective emissivity of the pixel on the urban earth surface:
in the formula, uinIs a sky visible factor 1-upsilon in the Landsat8 pixelinIs the part of visible sky shaded by buildings in the picture element; when the sky visible factor of the pixel is close to 0, the sky is almost invisible to the whole pixel, radiation is reflected for countless times in the pixel, the emission radiation of the pixel tends to 0, the pixel is approximate to a black body at the moment, and the effective emissivity of the pixel is 1; when the sky visibility factor of a picture element is equal to 1, meaning that the picture element is a flat ground surface without any building occlusion, the effective emissivity of the picture element is equal to the emissivity of the material.
The surface temperature inversion method also comprises the steps of generating an urban surface coverage type graph before calculating the surface emissivity:
generating a land cover type graph with the spatial resolution of 3m by using the high-spatial-resolution satellite image data with the spatial resolution of 3m and the building and water body vector data; in consideration of urban surface characteristics, the urban land coverage types are divided into five types of buildings, roads, grasslands, trees and water bodies; the satellite images with high spatial resolution and the building vector data are downloaded from an Arctiler platform; the high spatial resolution satellite imagery is from a QuickBird commercial satellite; the water vector data is obtained from an open street map; and superposing the building and water body vector data on the high-spatial-resolution satellite image, further distinguishing roads, grasslands and trees by using a supervision and classification method based on the high-spatial-resolution satellite image, and using the generated urban land cover type map for calculating the effective emissivity of the urban pixels.
The method for inverting the earth surface temperature comprises the following steps of: inputting the atmospheric profile data into atmospheric radiation transmission software MODTRAN to calculate atmospheric parameters such as atmospheric uplink radiation, atmospheric downlink radiation, atmospheric transmittance and the like;
the pixel-by-pixel calculation process of the atmospheric parameters comprises the following steps:
(1) selecting atmosphere profile data at the latest moment according to the transit time of the Landsat8 satellite;
(2) extracting height, temperature and relative humidity data of 37 atmospheric pressure layers of the ERA5 atmospheric profile, inputting the atmospheric profile into an atmospheric radiation transmission model MODTRAN, and calculating three parameters of atmospheric transmittance, atmospheric uplink radiation and atmospheric downlink radiation corresponding to each atmospheric profile;
(3) according to the surface elevation and the central longitude and latitude value of the pixel, performing linear interpolation on the atmospheric parameters of four grid points of ERA5 closest to the pixel to the same elevation, and then performing spatial linear interpolation on the atmospheric parameters of the same elevation plane;
(4) and according to the transit time of the Landsat8 satellite, performing linear interpolation on the atmospheric parameters subjected to elevation and spatial interpolation in time to finally obtain the atmospheric parameters corresponding to each Landsat8 pixel.
The invention has the following effective effects: by introducing city morphological parameters such as sky visible factors and the like into a conventional thermal radiation transmission equation and considering the influence of a city three-dimensional geometric structure on surface thermal radiation transmission, a thermal radiation transmission equation suitable for the city surface is established, and high-precision city surface temperature is inverted. Compared with the surface temperature inversion result of the conventional thermal radiation transmission equation, the method has the obvious advantages in the urban surface with the complex three-dimensional geometrical structure, and can make up for larger errors caused by the inversion of the surface temperature in the prior art. Fig. 4 shows that the surface temperature of the prior art inversion can cause an error of 1K in an urban area, which further confirms the effective effect of the present invention.
Drawings
FIG. 1 is a flow chart of a surface temperature inversion method considering the effects of urban three-dimensional geometry;
FIG. 2 is a digital surface model map of a study area;
fig. 3 is a map of the inversion result of the earth surface temperature considering the influence of the urban three-dimensional geometry, (a) 8 and 17 days in 2019, and (b) 12 and 7 days in 2019;
fig. 4 is a temperature difference diagram between the inversion result of the earth surface temperature considering the influence of the urban three-dimensional geometrical structure and the inversion result of the earth surface temperature based on the conventional thermal radiation transmission equation, (a) 8 and 17 days in 2019, and (b) 12 and 7 days in 2019;
Detailed Description
The present invention will be described in detail with reference to specific examples.
Step 1: on-board radiance product acquisition and data preprocessing
Download the on-board radiance data of TIRS (thermal information sensor) sensor carried on Landsat8 satellite. Data download website: https:// earth xplor. Converting the on-satellite radiance data into an on-satellite radiance value by using a formula (1):
L＝gain×DN+offset (1)
in the formula, L is the on-satellite radiance, DN is the count value, gain and offset are gain and offset, respectively. For Landsat8 channel 10, the gain and offset are 0.00033420 and 0.1, respectively.
And 2, step: sky visible factor calculation
The Sky Visibility Factor (SVF) is calculated using Digital Surface Model (DSM) data. DSM data with the spatial resolution of 3m is obtained from a natural resource satellite remote sensing cloud service platform (http:// www.sasclouds.com/EN). The data is generated by a three-dimensional satellite image of a domestic resource III, and height information of buildings and plants can be provided.
The sky visible factor represents the visible degree of the earth's surface to the sky in the hemisphere space, that is, the ratio of the atmospheric downlink radiation received by the point to the atmospheric downlink radiation from the hemisphere space received when the point is flat and not shielded, the value range is 0 to 1, and the sky visible factor is calculated by the formula (2):
in the formula, gammaiIs the vertical elevation angle in a given search direction and n is the number of search directions.
The sky visibility factors were calculated by discretizing the azimuth of the entire hemispherical space into 16 directions that were uniformly distributed, 0 °,30 °,45 °,60 °,90 °,120 °,135 °,150 °,180 °,210 °,225 °,240 °,270 °,300 °,315 ° and 330 °, respectively.
And 3, step 3: urban surface coverage type map generation
And generating a land cover type map with the spatial resolution of 3m by using the high-spatial-resolution satellite image data with the spatial resolution of 3m and the building and water body vector data. Urban land cover types are classified into five types of buildings, roads, lawns, trees, and water bodies in consideration of urban surface characteristics. High spatial resolution satellite imagery and building vector data were downloaded from the Arctiler platform (http:// www.arctiler.com /). The high spatial resolution satellite imagery is from QuickBird commercial satellites. Water body vector data is obtained from open street maps (https:// www.openstreetmap.org /). And superposing the building and water body vector data on the high-spatial-resolution satellite image, further distinguishing roads, grasslands and trees by using a supervision and classification method based on the high-spatial-resolution satellite image, and using the generated urban land cover type map for calculating the effective emissivity of the urban pixels.
And 4, step 4: urban pixel effective emissivity calculation
And calculating the material emissivity of typical ground objects of the city by using the ASTER spectral library and the city waterproof material spectral library. In addition to natural surfaces (grass, trees and water), urban surfaces are mainly covered by impervious surfaces, such as roads made of asphalt or concrete, buildings made of concrete and glass, and the like. In order to obtain the material emissivity of typical urban features, the material emissivity of buildings is set to the average emissivity of masonry, concrete and glass, and the material emissivity of roads is set to the average emissivity of cement, concrete and masonry. The material emissivity of the masonry, cement, concrete and glass is obtained by convolution calculation of a sensor response function of Landsat8 and emission spectrum curves of the masonry, cement, concrete and glass in a spectrum library of the urban impervious material. The emissivity of the materials of the grassland, the trees and the water body in the natural earth surface is obtained by convolution calculation of a sensor response function of Landsat8 and emission spectrum curves of the grassland, the trees and the water body in an ASTER spectral library. For a flat ground surface, the emissivity of the image element is calculated by the weighted average of the emissivity of materials of different surface feature components in the image element. For the urban earth surface with a complex three-dimensional geometric structure, the pixel effective emissivity of the urban earth surface needs to be calculated by considering multiple reflections and scattering among different ground object components inside the pixel. The process for calculating the effective emissivity of the pixel on the urban surface comprises the following steps:
(1) when the urban three-dimensional geometrical structure is not considered, calculating the pixel emissivity of the flat earth surface according to the weighted average value of the material emissivity of different ground components:
in the formula, epsilonefEmissivity of pixels for flat ground surface, fkIs the area ratio of the kth ground feature component in the pixel, εkIs the material emissivity of the kth terrestrial object component in the picture element.
(2) For urban terrains with complex three-dimensional geometries, the pixel effective emissivity passes through the pixel emissivity of a flat terrane (i.e.,. epsilon. in equation (3))ef) The emissivity increment caused by multiple scattering and reflection (namely cavity effect) between different ground object components in the pixel is calculated as follows:
in the formula, epsiloneuEffective emissivity of pixel, delta epsilon, for urban surfacemIs the emissivity increase caused by multiple scattering and reflection between different terrestrial components inside the picture element.
And (3) considering multiple scattering and reflection among different surface material components in the pixel, and finally obtaining the effective emissivity of the pixel on the urban earth surface:
in the formula, uinIs a sky visible factor of 1-upsilon inside Landsat8 picture elementinIs the part of the visible sky within the picture element that is obscured by the building. When the sky visible factor of the pixel is close to 0, it meansThe sky is almost invisible to the whole pixel, radiation is reflected for countless times in the pixel, the emission radiation of the pixel tends to be 0, the pixel is similar to a black body at the moment, and the effective emissivity of the pixel is 1; when the sky visibility factor of a picture element is equal to 1, meaning that the picture element is a flat ground surface without any building occlusion, the effective emissivity of the picture element is equal to the emissivity of the material.
And 5: atmospheric parameter pixel-by-pixel calculation
ERA5 reanalysis profile data is the latest generation of reanalysis data released by the European Central Medium-Range Weather projections (ECMWF). ERA5 realizes the improvement of time and space resolution on the basis of the ERA-Interim of the previous generation. The spatial resolution of the data was 0.25 ° x 0.25 °, and the temporal resolution was 1 hour. ERA5 provides atmospheric temperature and humidity profile data for 37 pressure layers of 1hPa to 1000 hPa. And inputting the atmospheric profile data into atmospheric radiation transmission software MODTRAN to calculate atmospheric uplink radiation, atmospheric downlink radiation and atmospheric transmittance atmospheric parameters.
The pixel-by-pixel calculation process of the atmospheric parameters comprises the following steps:
(1) selecting atmosphere profile data at the latest moment according to the transit time of the Landsat8 satellite;
(2) extracting height, temperature and relative humidity data of 37 atmospheric pressure layers of the ERA5 atmospheric profile, inputting the atmospheric profile into an atmospheric radiation transmission model MODTRAN, and calculating three parameters of atmospheric transmittance, atmospheric uplink radiation and atmospheric downlink radiation corresponding to each atmospheric profile;
(3) according to the surface elevation and the central longitude and latitude value of the pixel, performing linear interpolation on the atmospheric parameters of four grid points of ERA5 closest to the pixel to the same elevation, and then performing spatial linear interpolation on the atmospheric parameters of the same elevation plane;
(4) and according to the transit time of the Landsat8 satellite, performing linear interpolation on the atmospheric parameters subjected to elevation and spatial interpolation in time to finally obtain the atmospheric parameters corresponding to each Landsat8 pixel.
Step 6: urban surface temperature inversion
In the thermal infrared spectral range, the conventional thermal radiation transport equation is expressed as:
LTOA＝[εefB(Ts)+(1-εef)Ld]τ+Lu (6)
in the formula, LTOAFor the radiance received at the top of the atmosphere, epsilonefEmissivity of pixels of flat ground surface, B is Planck function, TsIs the surface temperature, tau is the atmospheric permeability, LdFor atmospheric downlink radiation, LuRadiating up to the atmosphere.
The complex three-dimensional geometry of urban terrains changes the conventional thermal radiation transmission mode. On one hand, the adjacent pixels shield part of atmospheric radiation from the hemispherical direction, so that the atmospheric downlink radiation received by the earth surface is reduced; on the other hand, radiation from adjacent picture elements increases the total radiation received by the sensor. Therefore, the thermal radiation transfer equation considering the influence of urban three-dimensional geometry is expressed as:
LTOA＝[εeuB(Ts)+(1-εeu)(υoutLd+Ladj)]τ+Lu (7)
in the formula, LTOAFor the radiance received at the top of the atmosphere, epsiloneuEffective emissivity of pixel on urban ground surfaceoutIs a sky visible factor L between Landsat8 pixelsadjFor adjacent pixel radiation, expressed as:
in which subscripts A and B represent a target pixel and an adjacent pixel, respectively, LBIs the radiance of the adjacent pixel, theta is the angle between the connecting lines of the center points of the target pixel and the adjacent pixel, dSBIs the area of the adjacent picture element, rABIs the distance between the target picture element and the neighboring picture element, and N is the total number of neighboring picture elements.
The urban surface temperature based on the inversion of the thermal radiation transmission equation considering the influence of the urban three-dimensional geometrical structure is expressed as follows:
in the formula, TsSurface temperature, B, for inversion taking into account the effects of urban three-dimensional geometry-1Is the inverse operation of the planck function.
As can be seen from equation (9), in the case where other surface and atmospheric parameters are known, the radiance of the neighboring pixels is obtained first. The precondition for calculating the radiance of the adjacent pixel is to acquire the surface temperature of the adjacent pixel. However, the surface temperature of the adjacent pixels is not known in advance. In order to solve the problem, the surface temperature obtained by inverting the conventional thermal radiation transmission equation is used as an initial value of the surface temperature of the adjacent pixel, and the final urban surface temperature is obtained by an iterative calculation method. The surface temperature based on the inversion of the conventional thermal radiation transport equation is expressed as:
in the formula, Ts_initialThe flat surface temperature obtained for inversion.
The urban surface temperature inversion process comprises the following steps:
(1) using the radiance of the Landsat8 channel 10 in the step 1, the surface emissivity calculated in the step 4, and the atmospheric transmittance, the atmospheric uplink radiation and the atmospheric downlink radiation calculated in the step 5, and performing pixel-by-pixel inversion by using a formula (10) to obtain the flat surface temperature Ts_initialAnd taking the surface temperature at the moment as an initial value of the surface temperature of the adjacent pixel element.
(2) Inputting the initial value of the surface temperature of the adjacent pixel into a formula (8) to obtain the radiance of the adjacent pixel, inputting the sky visible factor obtained in the step 2, the effective emissivity of the pixel of the urban surface calculated in the step 4, the atmospheric transmittance, the atmospheric uplink radiation and the atmospheric downlink radiation calculated in the step 5 and the radiance of the Landsat8 channel 10 in the step 1, and simultaneously inputting the radiance of the common pixelFormula (9), pixel-by-pixel inversion is performed to obtain the urban surface temperature TsAnd taking the surface temperature at the moment as the surface temperature of the adjacent pixel for iterative computation of the radiance of the adjacent pixel.
(3) And repeating the process (2) for iteratively calculating the surface temperature of the adjacent pixel element until the difference between the surface temperatures of the adjacent pixel element in the previous and next two times is less than 0.1K or the iteration times is more than 3, and outputting a result to obtain the final urban surface temperature.
FIG. 3 is a graph of urban surface temperature results inverted by the method of the present invention. The graph shows that the spatial distribution of the surface temperature is related to the building height and the type of surface coverage, wherein the surface temperature is higher in areas with lower building height and lower in water and vegetation covered areas. Fig. 4 shows that the earth surface temperature inversion method based on the conventional thermal radiation transmission equation can cause the highest error of 1K to the inversion result of the earth surface temperature, and simultaneously effectively illustrates the necessity of considering the influence of the urban three-dimensional geometrical structure when accurately inverting the urban earth surface temperature, thereby highlighting the important role of the invention.
It will be understood that modifications and variations can be made by persons skilled in the art in light of the above teachings and all such modifications and variations are intended to be included within the scope of the invention as defined in the appended claims.
Claims (6)
1. A surface temperature inversion method considering the influence of an urban three-dimensional geometrical structure is characterized by comprising the following steps: the urban surface temperature inversion process comprises the following steps:
step 1: calculating the surface emissivity:
(1) when the urban three-dimensional geometrical structure is not considered, calculating the pixel emissivity of the flat earth surface according to the weighted average value of the material emissivity of different ground components:
in the formula, epsilonefEmissivity of pixels for flat ground surface, fkIs the first in the picture elementArea ratio of k ground substance components,. epsilonkIs the material emissivity of the kth terrestrial object component in the pixel;
(2) for the urban ground surface with a complex three-dimensional geometrical structure, the pixel effective emissivity is equal to the pixel emissivity epsilon of the flat ground surfaceefAnd calculating emissivity increment caused by multiple scattering and reflection among different surface object components in the pixel:
in the formula, epsiloneuEffective emissivity of pixel, delta epsilon, for urban surfacemThe emissivity increment caused by multiple scattering and reflection among different surface feature components in the pixel is shown, and N is the total number of adjacent pixels;
and (3) considering multiple scattering and reflection among different surface material components in the pixel, and finally obtaining the effective emissivity of the pixel on the urban earth surface:
in the formula, uinIs a sky visible factor 1-upsilon in the Landsat8 pixelinIs the part of visible sky shaded by buildings in the picture element; when the sky visible factor of the pixel is close to 0, the sky is almost invisible to the whole pixel, radiation is reflected for countless times in the pixel, the emission radiation of the pixel tends to 0, the pixel is approximate to a black body at the moment, and the effective emissivity of the pixel is 1; when the sky visibility factor of the pixel is equal to 1, the pixel is a flat ground surface without any building shelter, and the effective emissivity of the pixel is equal to the emissivity of a material;
and (3) performing pixel-by-pixel inversion by using the obtained radiance of the Landsat8 channel 10, the calculated surface emissivity, the calculated atmospheric transmittance, the atmospheric uplink radiation and the atmospheric downlink radiation to obtain the flat surface temperature Ts_initialAnd the surface temperature at that time is taken asThe initial value of the surface temperature of the adjacent pixel element is obtained;
equation (10):
in the formula, Ts_initialFor flat surface temperature obtained by inversion, B-1As an inverse of the Planck function, LTOAThe radiance received at the top of the atmosphere, τ being the atmospheric transmittance, LdFor atmospheric downlink radiation, LuFor upward radiation of the atmosphere, eeuEffective emissivity of pixels on the surface of the city;
step 2: inputting the initial value of the surface temperature of the adjacent pixel element into a formula (8) to obtain the radiance of the adjacent pixel element, calculating the sky visible factor obtained by calculation, the effective emissivity of the pixel element of the urban surface, the atmospheric transmittance, the atmospheric uplink radiation and the atmospheric downlink radiation obtained by calculation and the radiance of the Landsat8 10 th channel obtained by calculation, simultaneously inputting a formula (9), and inverting pixel by pixel to obtain the radiance of the urban surface temperature TsAnd the city surface temperature at the moment is used as the surface temperature of the adjacent pixel for iterative computation of the radiance of the adjacent pixel;
formula (8):
in the formula, subscripts A and B represent a target pixel and an adjacent pixel, respectively, LBTheta is the radiance of the adjacent pixel, theta is the angle between the connecting lines of the center points of the target pixel and the adjacent pixel, dSBIs the area of the adjacent picture element, rABThe distance between the target pixel and the adjacent pixel is defined, and N is the total number of the adjacent pixels;
formula (9):
in the formula, TsSurface temperature, B, for inversion taking into account the effects of urban three-dimensional geometry-1As an inverse of the Planck function, LTOAThe radiance received at the top of the atmosphere, LdFor atmospheric downlink radiation, LuFor upward radiation of the atmosphere, eeuEffective emissivity of pixel on city surface, tau is atmospheric permeability, upsilonoutIs a sky visible factor L between Landsat8 pixelsadjRadiating for adjacent picture elements;
and 3, step 3: and (3) repeating the step (2) for iteratively calculating the surface temperature of the adjacent pixel element until the difference between the surface temperatures of the adjacent pixel elements in the previous and next times is smaller than a certain value or the iteration times is larger than a certain value, and outputting a result to obtain the final urban surface temperature.
2. The surface temperature inversion method of claim 1, wherein the method for obtaining the radiance of the Landsat8 channel 10 is as follows: downloading the on-satellite radiance data of a TIRS sensor carried on a Landsat8 satellite; data download website: https:// earth xplor. usgs. gov/; converting the starboard radiance data into starboard radiance values using equation (1):
L＝gain×DN+offset (1)
where L is the on-satellite radiance, DN is the count value, gain and offset are the gain and offset, respectively, and for Landsat8 channel 10, the gain and offset are 0.00033420 and 0.1, respectively.
3. The surface temperature inversion method of claim 1, wherein the sky visibility factor is calculated by: the sky visible factor represents the visible degree of the earth's surface to the sky in the hemisphere space, that is, the ratio of the atmospheric downlink radiation received by the point to the atmospheric downlink radiation from the hemisphere space received when the point is flat and not shielded, the value range is 0 to 1, and the sky visible factor is calculated by the formula (2):
in the formula, gammaiIs the vertical elevation angle in a given search direction and n is the number of search directions.
4. The surface temperature inversion method of claim 3, wherein the sky visibility factors are calculated by discretizing the azimuth of the entire hemispherical space into 16 directions that are uniformly distributed, 0 °,30 °,45 °,60 °,90 °,120 °,135 °,150 °,180 °,210 °,225 °,240 °,270 °,300 °,315 ° and 330 °, respectively.
5. The surface temperature inversion method of claim 3, further comprising, before calculating the surface emissivity, a step of generating an urban surface coverage type map:
generating a land cover type graph with the spatial resolution of 3m by using the high-spatial-resolution satellite image data with the spatial resolution of 3m and the building and water body vector data; in consideration of urban surface characteristics, the urban land coverage types are divided into five types of buildings, roads, grasslands, trees and water bodies; the satellite images with high spatial resolution and the building vector data are downloaded from an Arctiler platform; the high spatial resolution satellite imagery is from a QuickBird commercial satellite; the water vector data is obtained from an open street map; and superposing the building and water body vector data on the high-spatial-resolution satellite image, further distinguishing roads, grasslands and trees by using a supervision and classification method based on the high-spatial-resolution satellite image, and using the generated urban land cover type map for calculating the effective emissivity of the urban pixels.
6. The surface temperature inversion method according to claim 1, wherein the calculation method of the atmospheric transmittance, the atmospheric upward radiation and the atmospheric downward radiation comprises the following steps: inputting the atmospheric profile data into atmospheric radiation transmission software MODTRAN to calculate atmospheric uplink radiation, atmospheric downlink radiation and atmospheric transmittance atmospheric parameters;
the pixel-by-pixel calculation process of the atmospheric parameters comprises the following steps:
(1) selecting atmosphere profile data at the latest moment according to the transit time of the Landsat8 satellite;
(2) extracting height, temperature and relative humidity data of 37 atmospheric pressure layers of the ERA5 atmospheric profile, inputting the atmospheric profile into an atmospheric radiation transmission model MODTRAN, and calculating three parameters of atmospheric transmittance, atmospheric uplink radiation and atmospheric downlink radiation corresponding to each atmospheric profile;
(3) according to the surface elevation and the central longitude and latitude value of the pixel, performing linear interpolation on the atmospheric parameters of four grid points of ERA5 closest to the pixel to the same elevation, and then performing spatial linear interpolation on the atmospheric parameters of the same elevation plane;
(4) and according to the transit time of the Landsat8 satellite, performing linear interpolation on the atmospheric parameters subjected to elevation and spatial interpolation in time to finally obtain the atmospheric parameters corresponding to each Landsat8 pixel.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111494964.2A CN114152350B (en) | 2021-12-09 | 2021-12-09 | Earth surface temperature inversion method considering urban three-dimensional geometrical structure influence |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111494964.2A CN114152350B (en) | 2021-12-09 | 2021-12-09 | Earth surface temperature inversion method considering urban three-dimensional geometrical structure influence |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114152350A CN114152350A (en) | 2022-03-08 |
| CN114152350B true CN114152350B (en) | 2022-07-01 |
Family
ID=80454324
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111494964.2A Active CN114152350B (en) | 2021-12-09 | 2021-12-09 | Earth surface temperature inversion method considering urban three-dimensional geometrical structure influence |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114152350B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115186339B (en) * | 2022-07-05 | 2023-08-08 | 中国农业科学院农业资源与农业区划研究所 | Surface temperature and emissivity simultaneous inversion method considering influence of urban three-dimensional structure |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104748857A (en) * | 2015-03-05 | 2015-07-01 | 北京师范大学 | Method and system for inverting urban surface temperatures |
| CN109635477A (en) * | 2018-12-20 | 2019-04-16 | 中国农业科学院农业资源与农业区划研究所 | A kind of thermal infrared radiation mode for considering to close on pixel effect |
| CN111198162A (en) * | 2020-01-09 | 2020-05-26 | 胡德勇 | Remote sensing inversion method for urban surface reflectivity |
| CN112487346A (en) * | 2020-10-26 | 2021-03-12 | 中国农业科学院农业资源与农业区划研究所 | Mountain land surface temperature remote sensing retrieval method |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106909722B (en) * | 2017-02-10 | 2019-07-26 | 广西壮族自治区气象减灾研究所 | A kind of accurate inversion method of large area of temperature near the ground |
-
2021
- 2021-12-09 CN CN202111494964.2A patent/CN114152350B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104748857A (en) * | 2015-03-05 | 2015-07-01 | 北京师范大学 | Method and system for inverting urban surface temperatures |
| CN109635477A (en) * | 2018-12-20 | 2019-04-16 | 中国农业科学院农业资源与农业区划研究所 | A kind of thermal infrared radiation mode for considering to close on pixel effect |
| CN111198162A (en) * | 2020-01-09 | 2020-05-26 | 胡德勇 | Remote sensing inversion method for urban surface reflectivity |
| CN112487346A (en) * | 2020-10-26 | 2021-03-12 | 中国农业科学院农业资源与农业区划研究所 | Mountain land surface temperature remote sensing retrieval method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114152350A (en) | 2022-03-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Jones et al. | Storm-scale data assimilation and ensemble forecasting with the NSSL experimental Warn-on-Forecast system. Part II: Combined radar and satellite data experiments | |
| Islam et al. | A physics-based algorithm for the simultaneous retrieval of land surface temperature and emissivity from VIIRS thermal infrared data | |
| Ren et al. | Atmospheric water vapor retrieval from Landsat 8 thermal infrared images | |
| Ren et al. | Improving land surface temperature and emissivity retrieval from the Chinese Gaofen-5 satellite using a hybrid algorithm | |
| Hofierka et al. | Physically-based land surface temperature modeling in urban areas using a 3-D city model and multispectral satellite data | |
| CN112487346B (en) | Mountain land surface temperature remote sensing retrieval method | |
| CN109297605B (en) | Surface temperature inversion method based on mid-infrared and thermal infrared data | |
| Zhong et al. | Cross-Calibration of HJ-1/CCD over a desert site using landsat ETM $+ $ imagery and ASTER GDEM product | |
| Masiello et al. | Kalman filter physical retrieval of surface emissivity and temperature from SEVIRI infrared channels: A validation and intercomparison study | |
| Adderley et al. | The effect of radiometer placement and view on inferred directional and hemispheric radiometric temperatures of an urban canopy | |
| CN107389617A (en) | The inversion method and equipment of aerosol optical depth based on No. four satellites of high score | |
| CN114152350B (en) | Earth surface temperature inversion method considering urban three-dimensional geometrical structure influence | |
| Haurant et al. | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica | |
| Wang et al. | Estimating top-of-atmosphere daily reflected shortwave radiation flux over land from MODIS data | |
| CN114581349A (en) | Visible light image and infrared image fusion method based on radiation characteristic inversion | |
| Yang et al. | A physics-based algorithm to couple CYGNSS surface reflectivity and SMAP brightness temperature estimates for accurate soil moisture retrieval | |
| Liu et al. | Simultaneous retrieval of land surface temperature and emissivity from the FengYun-4A advanced geosynchronous radiation imager | |
| Zeng et al. | Land surface temperature and emissivity retrieval from nighttime middle and thermal infrared images of Chinese Fengyun-3D MERSI-II | |
| He et al. | Direct estimation of land surface albedo from simultaneous MISR data | |
| Oyoshi | Hourly LST monitoring with Japanese geostationary satellite MTSAT-1R over the Asia-Pacific region | |
| Majumdar et al. | Generation of emissivity and land surface temperature maps using MODIS TIR data for lithological mapping over the Singhbhum-Orissa Craton | |
| Zhou et al. | Land surface albedo estimation with Chinese GF-1 WFV data in Northwest China | |
| CN114296061A (en) | Cross calibration method based on multivariate variable detection and different radiation transmission models | |
| Song et al. | Estimation of land surface temperature using FengYun-2E (FY-2E) data: A case study of the source area of the yellow river | |
| Zhao et al. | An operational land surface temperature retrieval methodology for Chinese second-generation Huanjing disaster monitoring satellite data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |