A method for surface temperature inversion in African desert regions

Abstract

The invention relates to a method for inverting land surface temperature in desert areas in Africa. The method uses target satellite remote sensing image data and water vapor content data to respectively simulate a land surface emissivity matrix and a water vapor content matrix to improve the precision of the land surface emissivity and the water vapor content. The method combines a universal single-channel algorithm and a near-infrared band image in a target satellite remote sensing image to invert the land surface temperature in a target area. Different combinations of homogenized vegetation index thresholds are used to target the land surface coverage in the area to improve the land surface temperature inversion precision. The proposed design method will provide reference and reference for the inversion of large-scale land surface temperature, especially the inversion of land surface temperature in desert arid areas.

Description

Technical Field

[0001] The invention relates to a surface temperature inversion method for African desert areas, belonging to the technical field of surface temperature inversion. Background Art

[0002] At present, there is a lack of high-precision surface temperature data in desert areas, mainly due to the insufficient accuracy of surface emissivity and atmospheric water vapor content data. In the existing technology, the inversion of surface temperature using Landsat thermal infrared data mainly includes three parts: surface emissivity calculation, atmospheric water vapor content calculation and inversion algorithm. Among them, the important parameter of surface emissivity is not only related to the physical properties of the material itself, but also to a series of factors such as land cover type, vegetation structure, soil moisture, etc. in the region; the parameter of atmospheric water vapor content varies in different climate zones and even microclimate zones, so there is a problem of applicability of the surface temperature inversion model in geographically heterogeneous areas.

[0003] In addition, existing large-scale surface temperature inversion methods all use a fixed NDVI threshold combination, and use NCEP AWV data with lower temporal and spatial resolution when simulating atmospheric water vapor content. However, the fixed NDVI threshold combination cannot adapt to different geographical regions and different seasonal changes, and the NCEP AWV data with lower temporal and spatial resolution may lead to inaccurate atmospheric water vapor content in the simulation, affecting the accuracy of surface temperature inversion. Summary of the invention

[0004] The technical problem to be solved by the present invention is to provide a surface temperature inversion method for African desert areas, which can improve the accuracy of surface emissivity and the temporal and spatial resolution and accuracy of atmospheric water vapor content.

[0005] In order to solve the above technical problems, the present invention adopts the following technical solution: The present invention designs a surface temperature inversion method for African desert areas, which is used to obtain the surface temperature distribution of the target area corresponding to the target historical moment, comprising the following steps:

[0006] Step A. Obtain a target satellite remote sensing image containing a red band and a near-infrared band at a target historical moment corresponding to the target area, and obtain the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image; then obtain the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset based on the target surface emissivity dataset corresponding to the target historical moment of the target area, and then proceed to step B;

[0007] Step B. Taking the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image as the main analysis, and the normalized vegetation index of each pixel with respect to the target surface emissivity data set as the auxiliary analysis, the surface emissivity of each pixel in the target satellite remote sensing image is obtained, and the surface emissivity matrix of all pixels in the target satellite remote sensing image is formed, and then entering step C;

[0008] Step C. Obtain water vapor content data at each sampling time point at a preset sampling frequency within a historical time range corresponding to a time span of the target area including the target historical moment, and obtain the water vapor content value of each pixel in the target satellite remote sensing image corresponding to the target historical moment, to form a water vapor content matrix of all pixels in the target satellite remote sensing image, and then proceed to step D;

[0009] Step D. increasing and updating the on-board radiance of the near-infrared band image in the target satellite remote sensing image, and then calculating and obtaining the on-board brightness temperature of the image capturing sensor position corresponding to the target satellite remote sensing image, and then proceeding to step E;

[0010] Step E. Based on the surface emissivity matrix and water vapor content matrix of all pixels in the target satellite remote sensing image, the on-board radiance of the near-infrared band image in the target satellite remote sensing image, and the on-board brightness temperature of the image shooting sensor position corresponding to the target satellite remote sensing image, a universal single-channel algorithm is applied to calculate the surface temperature matrix of all pixels in the target satellite remote sensing image, which is the surface temperature distribution of the target area corresponding to the target historical moment.

[0011] As a preferred technical solution of the present invention: in step A, for the target satellite remote sensing image, according to the following formula:

[0012] (1)

[0013] Obtain the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image ,in, , Represents the number of pixels in the target satellite remote sensing image. Indicates the first The normalized difference vegetation index of the target satellite remote sensing image is Indicates the first The reflectance value of the pixel with respect to the red band, Indicates the first The reflectance value of each pixel in the near-infrared band;

[0014] According to the target surface emissivity dataset of the target area corresponding to the target historical moment with respect to the global surface emissivity dataset, the average value of the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to each band in the target surface emissivity dataset is obtained as the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset. , Indicates the first The normalized difference vegetation index of the target surface emissivity dataset for each pixel.

[0015] As a preferred technical solution of the present invention: in the step A, a Landsat8 satellite remote sensing image of the target area corresponding to the target historical moment is obtained to form a target satellite remote sensing image of the target area corresponding to the target historical moment including a red band and a near-infrared band.

[0016] As a preferred technical solution of the present invention: the step B includes the following steps B1 to B4;

[0017] Step B1. For each pixel in the target satellite remote sensing image, the normalized vegetation index of the pixel with respect to the target satellite remote sensing image is calculated. , according to the following formula:

[0018] (2-1)

[0019] Obtain the vegetation coverage of the pixel with respect to the target satellite remote sensing image ,in, , Represents the number of pixels in the target satellite remote sensing image. Indicates the first The vegetation coverage of pixels in the target satellite remote sensing image is Target satellite remote sensing image The normalized vegetation index of pixels on the target satellite remote sensing image, Represents the preset soil NDVI threshold for the target satellite remote sensing image. Indicates the NDVI threshold of the preset vegetation for the target satellite remote sensing image;

[0020] At the same time, the normalized vegetation index of the pixel for the target surface emissivity dataset , according to the following formula:

[0021] (2-2)

[0022] Obtain the vegetation coverage of the pixel with respect to the target surface emissivity dataset ,in, , Represents the number of pixels in the target satellite remote sensing image. Indicates the first The vegetation coverage of pixels with respect to the target surface emissivity dataset, Target satellite remote sensing image Normalized Difference Vegetation Index of pixels with respect to the target surface emissivity dataset;

[0023] Step B2. Definition =0, the surface emissivity of the pixel in the target satellite remote sensing image is equal to the preset pure water body emissivity, and the definition =1, the surface emissivity of the pixel in the target satellite remote sensing image is equal to the preset pure vegetation emissivity, and the definition and The pixels in the corresponding target satellite remote sensing image are mixed pixels;

[0024] Step B3. For each mixed pixel in the target satellite remote sensing image, execute the following steps B3-1 to B3-3 to obtain the surface emissivity of each mixed pixel;

[0025] Step B3-1. Follow the formula below:

[0026] (3-1)

[0027] (3-2)

[0028] Obtain the background emissivity of the mixed pixel corresponding to the 13th band of the target surface emissivity dataset , background emissivity in band 14 ,in, , Represents the number of mixed pixels in the target satellite remote sensing image. Indicates the first The mixed pixel corresponds to the background emissivity of the 13th band of the target surface emissivity dataset. Indicates the first The mixed pixel corresponds to the background emissivity of the 14th band of the target surface emissivity dataset. represents the emissivity of the 13th band of the global surface emissivity dataset, represents the emissivity of the 14th band of the global surface emissivity dataset, represents the pure vegetation emissivity with respect to the global land surface emissivity dataset, Indicates the first The vegetation coverage of the mixed pixels with respect to the target surface emissivity dataset;

[0029] Step B3-2. According to the following formula:

[0030] (4)

[0031] Obtain the background emissivity of the mixed pixel in the near infrared band of the target satellite remote sensing image , Indicates the first The background emissivity of the mixed pixel in the near-infrared band of the target satellite remote sensing image is , , They represent the regression coefficients corresponding to the target satellite remote sensing images;

[0032] Step B3-3. According to the following formula:

[0033] (5)

[0034] Obtain the emissivity of the mixed pixel in the near infrared band of the target satellite remote sensing image , as the surface emissivity of the mixed pixel, where Indicates the first The emissivity of the mixed pixels in the near-infrared band of the target satellite remote sensing image is Indicates the first The vegetation coverage of the mixed pixels on the target satellite remote sensing image, It represents the pure vegetation emissivity in the near-infrared band of the target satellite remote sensing image;

[0035] Step B4. The surface emissivity matrix of all pixels in the target satellite remote sensing image is constructed from the surface emissivity of all pixels in the target satellite remote sensing image.

[0036] As a preferred technical solution of the present invention: based on the Landsat8 satellite remote sensing image of the target area corresponding to the target historical moment obtained in step A, a target satellite remote sensing image of the target area corresponding to the target historical moment including a red band and a near infrared band is formed, in step B3-2, =0.682, =0.2578, =0.0584.

[0037] As a preferred technical solution of the present invention: the step C includes the following steps C1 to C3;

[0038] Step C1. Obtain the water vapor content data of each sampling time point at a preset sampling frequency within the historical time range corresponding to the time span of the target area and including the target historical moment, and determine whether the target historical moment is a sampling time point. If yes, obtain the water vapor content value of each pixel in the target satellite remote sensing image corresponding to the sampling time point, and form a water vapor content matrix of all pixels in the target satellite remote sensing image; otherwise, proceed to step C2;

[0039] Step C2. For each pixel in the target satellite remote sensing image, use the following formula:

[0040] (6)

[0041] Get the water vapor content value of the pixel corresponding to the target historical moment ,in, , Represents the number of pixels in the target satellite remote sensing image. Indicates the first Pixels correspond to the target historical moment The water vapor content value, Representation and target historical moments The sampling time points adjacent to the historical time direction, Representation and target historical moments The sampling time points in the adjacent future time direction, Indicates the first Pixels correspond to sampling time points The water vapor content value, Indicates the first Pixels correspond to sampling time points Then proceed to step C3;

[0042] Step C3. The water vapor content matrix of all pixels in the target satellite remote sensing image is constructed by corresponding to the water vapor content value of the target historical moment of all pixels in the target satellite remote sensing image.

[0043] As a preferred technical solution of the present invention: in the step C, the MERRA-2 product is used to obtain the water vapor content data of each sampling time point at a preset sampling frequency within the historical time range corresponding to the target area and including the target historical moment.

[0044] As a preferred technical solution of the present invention: the step D includes the following steps D1 to D2;

[0045] Step D1. For the near-infrared band image in the target satellite remote sensing image, according to the ee.Algorithms. Landsat.calibrated Radiance function of the GEE platform, the following formula is used:

[0046] (7)

[0047] Improve and update the on-board radiance of near-infrared band images in target satellite remote sensing images ,in, Represents the pixel value matrix of all pixels in the near-infrared band image of the updated target satellite remote sensing image. Represents the pixel value matrix of all pixels in the near-infrared band image of the target satellite remote sensing image before updating, Represents the preset gain constant under the ee.Algorithms. Landsat.calibrated Radiance function of the GEE platform. represents the preset bias constant under the ee. Algorithms. Landsat.calibrated Radiance function of the GEE platform, and then proceeds to step D2;

[0048] Step D2. Follow the formula below:

[0049] (8)

[0050] Obtain the on-board brightness temperature of the image capturing sensor position corresponding to the target satellite remote sensing image ,in, , They respectively represent the calibration constants of the image shooting sensor corresponding to the target satellite remote sensing image with respect to the near-infrared band.

[0051] As a preferred technical solution of the present invention: the step E includes the following steps E1 to E2;

[0052] Step E1. Based on the water vapor content matrix of all pixels in the target satellite remote sensing image , according to the following formula:

[0053] (9)

[0054] Get with Related parameters , , ,in, , , , , , , , , Represents the coefficient of the satellite remote sensing image type to which the target satellite remote sensing image belongs, and then proceeds to step E2;

[0055] Step E2. Follow the formula below:

[0056] (10)

[0057] Calculate the surface temperature matrix of all pixels in the target satellite remote sensing image , which is the surface temperature distribution of the target area at the target historical moment, where It represents the on-board radiance of the near-infrared band image in the target satellite remote sensing image. Indicates the on-board brightness temperature of the target satellite remote sensing image at the image capturing sensor location. Represents the surface emissivity matrix of all pixels in the target satellite remote sensing image, It represents the quantity equivalent to the wavelength of the image shooting sensor corresponding to the target satellite remote sensing image.

[0058] As a preferred technical solution of the present invention: based on the Landsat8 satellite remote sensing image of the target area corresponding to the target historical moment obtained in step A, in step E1, =0.04019, =0.02916, =1.01523, =-0.38333, =-1.50294, =0.20324, =0.00918, =1.36072, =-0.27514.

[0059] The surface temperature inversion method for African desert areas described in the present invention has the following technical effects compared with the prior art by using the above technical solution:

[0060] The present invention designs a surface temperature inversion method for African desert areas. The target satellite remote sensing image data and water vapor content data are used to simulate the surface emissivity matrix and the water vapor content matrix respectively, so as to improve the accuracy of the surface emissivity and the water vapor content. The surface temperature of the target area is inverted by combining a universal single-channel algorithm and the near-infrared band image of the target satellite remote sensing image. Different homogenized vegetation index threshold combinations are used to improve the surface temperature inversion accuracy according to the surface coverage of the area. The proposed design method will provide a reference and reference for the inversion of large-scale surface temperature, especially the inversion of surface temperature in desert arid areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 It is a schematic diagram of an embodiment area of ​​the method for inverting land surface temperature in desert areas of Africa designed by the present invention;

[0062] Figure 2 is a schematic diagram comparing the root mean square error of the design method of the present invention and three prior art methods;

[0063] Figure 3 This is a schematic diagram of the accuracy comparison corresponding to 24 threshold combinations in the total study area;

[0064] Figure 4 It is the combined zoning map of the highest threshold accuracy of the total study area in each group of NDVIveg. DETAILED DESCRIPTION

[0065] The specific implementation modes of the present invention will be further described in detail below in conjunction with the accompanying drawings.

[0066] The present invention designs a surface temperature inversion method for African desert areas, which is used to obtain the surface temperature distribution of the target area corresponding to the target historical moment. In actual application, the specific design executes the following steps A to E.

[0067] Step A. Obtain a target satellite remote sensing image containing a red band and a near-infrared band at a target historical moment corresponding to the target area, and obtain the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image; then obtain the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset based on the target surface emissivity dataset corresponding to the target historical moment of the target area, and then proceed to step B.

[0068] In practical applications, in the above step A, the Landsat8 satellite remote sensing image of the target area corresponding to the target historical moment is obtained to form the target satellite remote sensing image of the target area corresponding to the target historical moment including the red band and the near infrared band. For the target satellite remote sensing image, the following formula is used:

[0069] (1)

[0070] Obtain the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image ,in, , Represents the number of pixels in the target satellite remote sensing image. Indicates the first The normalized vegetation index of pixels on the target satellite remote sensing image, Indicates the first The reflectance value of the pixel with respect to the red band, Indicates the first The reflectance value of each pixel in the near-infrared band.

[0071] Then, according to the target surface emissivity dataset of the target area corresponding to the target historical moment with respect to the global surface emissivity dataset, the average value of the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to each band in the target surface emissivity dataset is obtained as the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset. , Indicates the first The normalized difference vegetation index of the target surface emissivity dataset for each pixel.

[0072] Step B. Take the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image as the main analysis, and use the normalized vegetation index of each pixel with respect to the target surface emissivity data set as the auxiliary analysis to obtain the surface emissivity of each pixel in the target satellite remote sensing image, and form the surface emissivity matrix of all pixels in the target satellite remote sensing image, and then enter step C.

[0073] In actual application, the above step B is specifically designed to execute the following steps B1 to B4.

[0074] Step B1. For each pixel in the target satellite remote sensing image, the normalized vegetation index of the pixel with respect to the target satellite remote sensing image is calculated. , according to the following formula:

[0075] (2-1)

[0076] Obtain the vegetation coverage of the pixel with respect to the target satellite remote sensing image ,in, , Represents the number of pixels in the target satellite remote sensing image. Indicates the first The vegetation coverage of pixels in the target satellite remote sensing image is Target satellite remote sensing image The normalized vegetation index of pixels on the target satellite remote sensing image, Indicates the preset soil NDVI threshold value of the target satellite remote sensing image = 0.2, Indicates that the preset vegetation NDVI threshold value for the target satellite remote sensing image is 0.5.

[0077] At the same time, the normalized vegetation index of the pixel for the target surface emissivity dataset , according to the following formula:

[0078] (2-2)

[0079] Obtain the vegetation coverage of the pixel with respect to the target surface emissivity dataset ,in, , Represents the number of pixels in the target satellite remote sensing image. Indicates the first The vegetation coverage of pixels with respect to the target surface emissivity dataset, Target satellite remote sensing image The normalized difference vegetation index of the target surface emissivity dataset for each pixel.

[0080] Step B2. Definition =0, the surface emissivity of the pixel in the target satellite remote sensing image is equal to the preset pure water body emissivity, and the definition =1, the surface emissivity of the pixel in the target satellite remote sensing image is equal to the preset pure vegetation emissivity, and the definition and The pixels in the corresponding target satellite remote sensing image are mixed pixels.

[0081] Step B3. For each mixed pixel in the target satellite remote sensing image, execute the following steps B3-1 to B3-3 to obtain the surface emissivity of each mixed pixel.

[0082] Step B3-1. Follow the formula below:

[0083] (3-1)

[0084] (3-2)

[0085] Obtain the background emissivity of the mixed pixel corresponding to the 13th band of the target surface emissivity dataset , background emissivity in band 14 ,in, , Represents the number of mixed pixels in the target satellite remote sensing image. Indicates the first The mixed pixel corresponds to the background emissivity of the 13th band of the target surface emissivity dataset. Indicates the first The mixed pixel corresponds to the background emissivity of the 14th band of the target surface emissivity dataset. represents the emissivity of the 13th band of the global surface emissivity dataset, represents the emissivity of the 14th band of the global surface emissivity dataset, Indicates the pure vegetation emissivity of the global surface emissivity dataset = 0.986, Indicates the first The vegetation coverage of the mixed pixels with respect to the target surface emissivity dataset.

[0086] Step B3-2. According to the following formula:

[0087] (4)

[0088] Obtain the background emissivity of the mixed pixel in the near infrared band of the target satellite remote sensing image , Indicates the first The background emissivity of the mixed pixel in the near-infrared band of the target satellite remote sensing image is , , They represent the regression coefficients corresponding to the target satellite remote sensing images. In practical applications, in step A, the target area is composed of Landsat8 satellite remote sensing images corresponding to the target historical moment, which includes the red band and the near-infrared band. Here =0.682, =0.2578, =0.0584.

[0089] Step B3-3. According to the following formula:

[0090] (5)

[0091] Obtain the emissivity of the mixed pixel in the near infrared band of the target satellite remote sensing image , as the surface emissivity of the mixed pixel, where Indicates the first The emissivity of the mixed pixels in the near-infrared band of the target satellite remote sensing image is Indicates the first The vegetation coverage of the mixed pixels on the target satellite remote sensing image, Represents the pure vegetation emissivity in the near-infrared band of the target satellite remote sensing image.

[0092] Step B4. The surface emissivity matrix of all pixels in the target satellite remote sensing image is constructed from the surface emissivity of all pixels in the target satellite remote sensing image.

[0093] Step C. Obtain the water vapor content data of each sampling time point at the preset sampling frequency within the historical time range corresponding to the time span of the target area including the target historical moment, and obtain the water vapor content value of each pixel in the target satellite remote sensing image corresponding to the target historical moment, to form a water vapor content matrix of all pixels in the target satellite remote sensing image, and then enter step D.

[0094] In actual application, the above step C is specifically designed to execute the following steps C1 to C3.

[0095] Step C1. Apply the MERRA-2 product to obtain the water vapor content data of each sampling time point at the preset sampling frequency within the historical time range corresponding to the time span of the target area including the target historical moment, and determine whether the target historical moment is a sampling time point. If so, obtain the water vapor content value of each pixel in the target satellite remote sensing image corresponding to the sampling time point to form a water vapor content matrix of all pixels in the target satellite remote sensing image; otherwise, proceed to step C2.

[0096] In practical applications, since the spatial resolution of MERRA-2 products is 0.625°×0.5°, it is necessary to perform spatial interpolation on MERRA-2 products to obtain the atmospheric water vapor content corresponding to each pixel in the remote sensing image of the target area. Here, the ee. Image. resample function provided by the GEE platform is specifically designed, and the interpolation mode is selected as bilinear interpolation.

[0097] Step C2. For each pixel in the target satellite remote sensing image, use the following formula:

[0098] (6)

[0099] Get the water vapor content value of the pixel corresponding to the target historical moment ,in, , Represents the number of pixels in the target satellite remote sensing image. Indicates the first Pixels correspond to the target historical moment The water vapor content value, Representation and target historical moments The adjacent sampling time points in the historical time direction, Representation and target historical moments The sampling time points in the adjacent future time direction, Indicates the first Pixels correspond to sampling time points The water vapor content value, Indicates the first Pixels correspond to sampling time points Then enter step C3.

[0100] Step C3. The water vapor content matrix of all pixels in the target satellite remote sensing image is constructed by corresponding to the water vapor content value of the target historical moment of all pixels in the target satellite remote sensing image.

[0101] Step D: Increase and update the on-board radiance of the near-infrared band image in the target satellite remote sensing image, and then calculate and obtain the on-board brightness temperature of the image shooting sensor position corresponding to the target satellite remote sensing image, and then enter step E.

[0102] In practical application, the above step D is specifically designed to execute the following steps D1 to D2.

[0103] Step D1. For the near-infrared band image in the target satellite remote sensing image, according to the ee.Algorithms. Landsat.calibrated Radiance function of the GEE platform, the following formula is used:

[0104] (7)

[0105] Improve and update the on-board radiance of near-infrared band images in target satellite remote sensing images ,in, Represents the pixel value matrix of all pixels in the near-infrared band image of the updated target satellite remote sensing image. Represents the pixel value matrix of all pixels in the near-infrared band image of the target satellite remote sensing image before updating, Represents the preset gain constant under the ee.Algorithms. Landsat.calibrated Radiance function of the GEE platform. It represents the preset bias constant under the ee. Algorithms. Landsat.calibrated Radiance function of the GEE platform, and then goes to step D2.

[0106] Step D2. Follow the formula below:

[0107] (8)

[0108] Obtain the on-board brightness temperature of the image capturing sensor position corresponding to the target satellite remote sensing image ,in, , They represent the calibration constants of the image capturing sensor corresponding to the target satellite remote sensing image with respect to the near-infrared band. For the near-infrared band, =774.89 W / m2·sr·um, =1321.08K.

[0109] Step E. Based on the surface emissivity matrix and water vapor content matrix of all pixels in the target satellite remote sensing image, the on-board radiance of the near-infrared band image in the target satellite remote sensing image, and the on-board brightness temperature of the image shooting sensor position corresponding to the target satellite remote sensing image, a universal single-channel algorithm is applied to calculate the surface temperature matrix of all pixels in the target satellite remote sensing image, which is the surface temperature distribution of the target area corresponding to the target historical moment.

[0110] In practical application, the above step E is specifically designed to execute the following steps E1 to E2.

[0111] Step E1. Based on the water vapor content matrix of all pixels in the target satellite remote sensing image , according to the following formula:

[0112] (9)

[0113] Get with Related parameters , , ,in, , , , , , , , , represents the coefficient of the satellite remote sensing image type to which the target satellite remote sensing image belongs, and then proceeds to step E2; and in conjunction with step A, obtains the Landsat8 satellite remote sensing image of the target area corresponding to the target historical moment, =0.04019, =0.02916, =1.01523, =-0.38333, =-1.50294, =0.20324, =0.00918, =1.36072, =-0.27514.

[0114] Step E2. Follow the formula below:

[0115] (10)

[0116] Calculate the surface temperature matrix of all pixels in the target satellite remote sensing image , which is the surface temperature distribution of the target area at the target historical moment, where It represents the on-board radiance of the near-infrared band image in the target satellite remote sensing image. Indicates the on-board brightness temperature of the target satellite remote sensing image at the image capturing sensor location. Represents the surface emissivity matrix of all pixels in the target satellite remote sensing image, It represents the quantity equivalent to the wavelength of the image shooting sensor corresponding to the target satellite remote sensing image.

[0117] The surface temperature inversion method designed by the present invention for the African desert region is used. In order to carefully reflect the performance of different NDVI threshold combinations in the algorithm, a total of 16 verification areas with three main surface cover types in the total study area are specially selected. The detailed information is as follows Figure 1 shown.

[0118] The accuracy of the design method of the present invention, SMW Algorithm and PSC Algorithm was verified with USGS Landsat LST in the total study area and each surface coverage type verification area. The results are as follows: Figure 2 As shown in the figure, in the desert area, the RMSE value of the PSC algorithm is generally high, especially in the verification area 3, reaching 2.11K, which is the highest error in all verification areas; the RMSE of the SMW algorithm in most desert areas is smaller than that of PSC and GSC, and the value is more stable; the RMSE value of the GSC algorithm in the desert area is generally the lowest.

[0119] connect Figure 1 ,like Figure 2 As shown in the figure, the RMSE value of the built-up area is lower than that of the desert area, but there is still a certain difference; in the validation areas 1 and 6, the RMSE value of the SMW algorithm is the largest, while the PSC algorithm performs relatively well in all areas, especially in the validation areas 2 and 3, where the RMSE is small; the RMSE value of the GSC algorithm is at a medium level in most areas.

[0120] In addition, the RMSE values ​​in the farmland area vary greatly among the algorithms; the SMW algorithm has a higher RMSE in validation areas 2 and 7, but lower values ​​in other areas. The PSC algorithm performs relatively poorly in validation area 2, with a higher RMSE, but is relatively stable in other areas; the GSC algorithm performs well in validation areas 1, 3, and 5, with the lowest RMSE value.

[0121] Therefore, in the entire study area, the SMW algorithm has the highest RMSE value, followed by the PSC algorithm, and the GSC algorithm has the lowest overall RMSE value. In addition, different LST accuracy verification information was obtained by combining different NDVI thresholds. The results are as follows: Figure 3 shown.

[0122] connect Figure 1 In the validation areas with different land cover types, the threshold combination with the lowest RMSE value in the total study area within each NDVIVveg group was selected for validation of different land cover types. The results are as follows:

[0123] like Figure 4 As shown, the RMSE value in the desert area is generally low, and the changes in different NDVI thresholds in validation areas 1 and 2 have little effect on the results, with the RMSE value being around 0.2 K. In validation area 3, the RMSE value rises slightly, and the performance differences of different NDVI threshold combinations begin to emerge, with the lowest error when NDVIsoil = 0.75 and NDVIveg = 0.2 (green).

[0124] The overall RMSE value of the built-up area is not high, but the RMSE differences of different NDVI threshold combinations are relatively obvious in validation areas 2 and 3. In validation area 2, NDVIsoil = 0.75 and NDVIveg = 0.2.

[0125] The error under the (green) combination is the smallest. In general, as the NDVIsoil threshold increases (from 0.5 to 0.75), the RMSE value decreases in the urban construction area.

[0126] The RMSE value of the cropland area is the highest among all land surface types, especially in validation areas 1 and 2, where the error exceeds 1.5K. In the cropland area, the RMSE values ​​under the combination of NDVIsoil = 0.55, NDVIveg = 0.05 (light blue) and NDVIsoil = 0.6, NDVIveg = 0.2 (dark blue) are relatively low, which performs well. As the NDVIsoil threshold increases, the error decreases, but this trend is not obvious in validation areas 4, 5, and 6.

[0127] The present invention is designed to be used for a surface temperature inversion method in African desert areas. The target satellite remote sensing image data and water vapor content data are used to simulate the surface emissivity matrix and the water vapor content matrix respectively, so as to improve the accuracy of the surface emissivity and the water vapor content. The surface temperature of the target area is inverted by combining a universal single-channel algorithm and near-infrared band images in the target satellite remote sensing image. Different homogenized vegetation index threshold combinations are used to target the surface coverage of the area, so as to improve the surface temperature inversion accuracy. The proposed design method will provide a reference and reference for the inversion of large-scale surface temperature, especially the inversion of surface temperature in desert arid areas.

[0128] The embodiments of the present invention are described in detail above with reference to the accompanying drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A surface temperature inversion method for African desert areas, used to obtain the surface temperature distribution of a target area corresponding to a target historical moment, characterized in that: The steps include: Step A. Obtain a target satellite remote sensing image containing a red band and a near infrared band at a target historical moment corresponding to the target area, and obtain a normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image; then obtain a normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset based on a target surface emissivity dataset corresponding to the target historical moment of the target area, and then proceed to step B; Step B. Taking the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image as the main analysis, and coordinating the normalized vegetation index of each pixel with respect to the target surface emissivity data set as the auxiliary analysis, the surface emissivity of each pixel in the target satellite remote sensing image is obtained, and the surface emissivity matrix of all pixels in the target satellite remote sensing image is formed, and then entering step C; In the above step B, firstly, based on the definition of the surface emissivity of pixels with vegetation coverage of 0 and 1 as the preset pure water body emissivity and the preset pure vegetation emissivity, for mixed pixels with vegetation coverage between 0 and 1, the formula and Obtain the background emissivity ε of the mixed pixel corresponding to the 13th band for the target surface emissivity dataset A13,n,bare , background emissivity ε of band 14 A14,n,bare , where 1≤n≤N, N represents the number of mixed pixels in the target satellite remote sensing image, ε A13 , ε A14 , ε A,veg They represent the 13th band emissivity, 14th band emissivity, and pure vegetation emissivity of the global surface emissivity dataset, respectively. A,n It represents the vegetation coverage of the nth mixed pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset; Then according to the formula ε LNIR,n,bare =c0ε A13,n,bare +c1ε A14,n,bare +c2, obtain the background emissivity ε of the nth mixed pixel in the target satellite remote sensing image with respect to the near infrared band in the target satellite remote sensing image LNIR,n,bare , c0, c1, c2 represent the regression coefficients corresponding to the target satellite remote sensing image; Finally, according to the formula ε LNIR,n =PV L,n ·ε LNIR,veg +(1-PV L,n )·ε LNIR,n,bare , obtain the emissivity ε of the mixed pixel about the near infrared band in the target satellite remote sensing image LNIR,n , where PV L,n represents the vegetation coverage of the nth mixed pixel in the target satellite remote sensing image with respect to the target satellite remote sensing image, ε LNIR,veg It represents the pure vegetation emissivity in the near-infrared band of the target satellite remote sensing image; The surface emissivity matrix of all pixels in the target satellite remote sensing image is formed by the surface emissivity of all pixels in the target satellite remote sensing image; Step C. Obtain water vapor content data at each sampling time point at a preset sampling frequency within a historical time range corresponding to a time span of the target area including the target historical moment, and obtain the water vapor content value of each pixel in the target satellite remote sensing image corresponding to the target historical moment, to form a water vapor content matrix of all pixels in the target satellite remote sensing image, and then proceed to step D; Step D. improving and updating the on-board radiance of the near-infrared band image in the target satellite remote sensing image, and then calculating and obtaining the on-board brightness temperature of the image shooting sensor position corresponding to the target satellite remote sensing image, and then entering step E; Step E. Based on the surface emissivity matrix and water vapor content matrix of all pixels in the target satellite remote sensing image, as well as the on-board radiance of the near-infrared band image in the target satellite remote sensing image and the on-board brightness temperature of the image shooting sensor position corresponding to the target satellite remote sensing image, a universal single-channel algorithm is applied to calculate and obtain the surface temperature matrix of all pixels in the target satellite remote sensing image, that is, the surface temperature distribution of the target area corresponding to the target historical moment; The above step E includes the following steps E1 to E2; Step E1. According to the water vapor content matrix w of all pixels in the target satellite remote sensing image, the following formula is used: Obtain parameters ψ1, ψ2, ψ3 related to w, where c 11 、c 12 、c 13 、c 21 、c 22 、c 23 、c 31 、c 32 、c 33 Represents the coefficient of the satellite remote sensing image type to which the target satellite remote sensing image belongs, and then proceeds to step E2; Step E2. According to the following formula: Calculate the surface temperature matrix T of all pixels in the target satellite remote sensing image s , which is the surface temperature distribution of the target area at the target historical moment, where L λ represents the on-board radiance of the near-infrared band image in the target satellite remote sensing image, BT represents the on-board brightness temperature of the image shooting sensor position corresponding to the target satellite remote sensing image, ε represents the surface emissivity matrix of all pixels in the target satellite remote sensing image, and b γ It represents the quantity equivalent to the wavelength of the image shooting sensor corresponding to the target satellite remote sensing image.

2. A method for inverting land surface temperature in desert areas of Africa according to claim 1, characterized in that: In step A, for the target satellite remote sensing image, the following formula is used: Obtain the normalized vegetation index NDVI of each pixel in the target satellite remote sensing image L,k , where 1≤k≤K, K represents the number of pixels in the target satellite remote sensing image, ρ LR,k represents the reflectance value of the kth pixel in the target satellite remote sensing image with respect to the red band, ρ LNIR,k It represents the reflectance value of the kth pixel in the target satellite remote sensing image with respect to the near-infrared band; According to the target surface emissivity dataset of the target area corresponding to the target historical moment with respect to the global surface emissivity dataset, the average value of the normalized vegetation index of each pixel in the target satellite remote sensing image with respect to each band in the target surface emissivity dataset is obtained as the normalized vegetation index NDVI of each pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset. A,k .

3. A method for inverting land surface temperature in desert areas of Africa according to claim 1 or 2, characterized in that: In the step A, a Landsat8 satellite remote sensing image of the target area corresponding to the target historical moment is obtained to form a target satellite remote sensing image of the target area corresponding to the target historical moment including a red band and a near infrared band.

4. The method for inverting land surface temperature in desert areas of Africa according to claim 1, characterized in that: The step B also includes the following step B1; Step B1. For each pixel in the target satellite remote sensing image, the normalized vegetation index NDVI of the pixel with respect to the target satellite remote sensing image is calculated. L,k , according to the following formula: Obtain the vegetation coverage PV of the pixel on the target satellite remote sensing image L,k , where 1≤k≤K, K represents the number of pixels in the target satellite remote sensing image, PV L,k Indicates the vegetation coverage of the kth pixel in the target satellite remote sensing image, NDVI L,k Normalized Difference Vegetation Index (NDVI) of the kth pixel in the target satellite remote sensing image soil Represents the NDVI threshold of the preset soil for the target satellite remote sensing image, NDVI veg Indicates the NDVI threshold of the preset vegetation for the target satellite remote sensing image; At the same time, the Normalized Difference Vegetation Index (NDVI) of the target surface emissivity dataset for the pixel A,k , according to the following formula: Obtain the vegetation coverage PV of the pixel with respect to the target surface emissivity dataset A,k , where PV A,k Represents the vegetation coverage of the kth pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset, NDVI A,k The normalized vegetation index of the kth pixel in the target satellite remote sensing image with respect to the target surface emissivity dataset.

5. A method for inverting land surface temperature in desert areas of Africa according to claim 4, characterized in that: Based on the Landsat8 satellite remote sensing image of the target area corresponding to the target historical moment obtained in step A, a target satellite remote sensing image of the target area corresponding to the target historical moment including a red band and a near-infrared band is constructed, c0=0.682, c1=0.2578, c2=0.0584.

6. The method for inverting land surface temperature in desert areas of Africa according to claim 1, characterized in that: The step C includes the following steps C1 to C3; Step C1. Obtain water vapor content data at each sampling time point at a preset sampling frequency within a historical time range corresponding to a time span of the target area including the target historical moment, and determine whether the target historical moment is a sampling time point. If yes, obtain the water vapor content value of each pixel in the target satellite remote sensing image corresponding to the sampling time point, and form a water vapor content matrix of all pixels in the target satellite remote sensing image; otherwise, proceed to step C2; Step C2. For each pixel in the target satellite remote sensing image, use the following formula: Get the water vapor content value AWV of the pixel corresponding to the target historical moment k,t , where 1≤k≤K, K represents the number of pixels in the target satellite remote sensing image, t1 represents the sampling time point adjacent to the target historical moment t in the direction of historical time, and t2 represents the sampling time point adjacent to the target historical moment t in the direction of future time. It represents the water vapor content value of the kth pixel in the target satellite remote sensing image corresponding to the sampling time point t1. Indicates the water vapor content value of the kth pixel in the target satellite remote sensing image corresponding to the sampling time point t2; then enters step C3; Step C3: The water vapor content matrix of all pixels in the target satellite remote sensing image is constructed by corresponding to the water vapor content value of the target historical moment of all pixels in the target satellite remote sensing image.

7. A method for inverting land surface temperature in desert areas of Africa according to claim 1 or 6, characterized in that: In step C, the MERRA-2 product is used to obtain water vapor content data at each sampling time point at a preset sampling frequency within a historical time range corresponding to a time span of the target area including the target historical moment.

8. The method for inverting land surface temperature in desert areas of Africa according to claim 1, characterized in that: The step D includes the following steps D1 to D2; Step D1. For the near-infrared band image in the target satellite remote sensing image, according to the ee.Algorithms.Landsat.calibrated Radiance function of the GEE platform, the following formula is used: L λ =gian×DN+offset (7) Improve and update the on-board radiance L of the near-infrared band image in the target satellite remote sensing image λ , where DN represents the pixel value matrix of all pixels in the near infrared band image of the target satellite remote sensing image before updating, gian represents the preset gain constant under the ee.Algorithms.Landsat.calibrated Radiance function of the GEE platform, and offset represents the preset bias constant under the ee.Algorithms.Landsat.calibrated Radiance function of the GEE platform, and then enters step D2; Step D2. According to the following formula: The on-board brightness temperature BT of the image capturing sensor position corresponding to the target satellite remote sensing image is obtained, wherein K1 and K2 respectively represent the calibration constants of the image capturing sensor corresponding to the target satellite remote sensing image with respect to the near-infrared band.

9. The method for inverting land surface temperature in desert areas of Africa according to claim 1, characterized in that: Based on the Landsat8 satellite remote sensing image of the target area corresponding to the target historical moment obtained in step A, in step E1, c 11 =0.04019, c 12 =0.02916, c 13 =1.01523, c 21 =-0.38333, c 22 =-1.50294, c 23 =0.20324, c 31 =0.00918, c 32 =1.36072, c 33 =-0.27514.

Images