A method for estimating surface temperature under clouds based on sun-cloud-satellite observation geometry

Abstract

The present invention belongs to the field of remote sensing application technology, and discloses a method for estimating the surface temperature under the cloud considering the geometry of sun-cloud-satellite observations, and specifically relates to an FY-3D MERSI-II cloud surface temperature estimation algorithm based on an integrated learning model. The estimation of the surface temperature under the cloud is mainly for the consideration of the geometric observation structure of the sun-surface-satellite of the FY-3D MERSI-II image, and the pixels under the cloud are divided into cloud-blocked pixels and cloud-covered pixels, and the surface temperature is estimated by using an integrated learning model for the different temperature formation mechanisms of the two parts of pixels. In the estimation process, the MODTRAN radiation transfer model is combined, and multiple surface temperature driving factors related to solar radiation, atmosphere, surface and precipitation meteorological conditions are integrated, and a set of process methods for realizing automatic calculation of surface temperature under the cloud is established, realizing the automatic all-weather surface temperature image estimation of the FY-3D MERSI-II image, and improving the applicability of the FY-3D MERSI-II high-resolution surface temperature data.

Description

Technical Field

[0001] The invention belongs to the technical field of remote sensing application, and in particular relates to a method for estimating surface temperature under clouds based on sun-cloud-satellite observation geometry. Background Art

[0002] Land surface temperature (LST) is a comprehensive representation of the energy exchange results of surface-atmosphere interaction at regional and global scales. It is widely used in the fields of surface energy balance, drought monitoring, farmland soil moisture management, crop yield estimation, numerical weather forecasting and climate change. Remote sensing observations are characterized by spatial continuity, short cycle and low cost, so they are an effective technical means to obtain spatial and temporal distribution information of LST at regional and even global scales. In actual inventions and applications, due to the convenience of measurement and mature clear sky LST inversion algorithms, surface temperature data are still mainly obtained through thermal infrared (TIR) ​​remote sensing observation data. At present, surface temperature products such as MODIS, GLASS, FY3C / VIRR (visible infrared scanning radiometer) based on thermal infrared remote sensing have been widely used in change detection, drought monitoring and climate change invention. However, TIR LST is not available in cloudy conditions. TIR signals cannot penetrate the clouds to reach the surface. The thermal radiation of cloud pixels observed by thermal infrared sensors mainly comes from the top of the clouds, rather than the thermal radiation of the surface under the clouds. The temperature of cloud pixels in thermal infrared remote sensing images is actually the cloud top temperature, not the surface temperature under the cloud. The lack of surface temperature information of cloud-covered pixels has largely destroyed the integrity and continuity of the spatial and temporal distribution of regional surface temperature, seriously hindering and limiting the application and development of surface temperature products. In many remote sensing applications, such as climate change investigations, evapotranspiration estimation, urban heat islands, etc., the problem of being unable to comprehensively and effectively monitor the invention area due to the inability to estimate the surface temperature of cloud pixels is often encountered. The lack of data caused by cloud cover is an important obstacle to the current application of surface temperature data. Comprehensively grasping the spatial pattern distribution and temporal evolution of all-weather LST has very important practical significance for the invention of related scientific fields such as surface resource and environmental monitoring, energy balance of the earth-atmosphere system and ecosystems.

[0003] In recent years, ensemble learning has become a new trend in the development of machine learning. By combining a series of weak learners into a new model with higher accuracy and stronger generalization ability, the defects of the original algorithm, such as single function and easy overfitting, can be improved. Random forest (RF), gradient boosted regression tree (GBRT) and autoregressive extreme gradient ascent algorithm (XGBoost) are emerging ensemble learning models developed in recent years. They are not only applicable to different scenarios, but also have outstanding accuracy. Among them, GBRT and XGBoost use boosting ensemble decision tree algorithms, while random forest (RF) uses bagging ensemble decision tree algorithms. These three ensemble algorithms have been widely used in various fields such as biology, medicine, geoscience, transportation, energy, agricultural engineering and industry, and have achieved relatively ideal classification or regression effects. The problem of estimating the surface temperature under the cloud belongs to a multivariate regression problem, so these three ensemble regression algorithms are used to simulate the surface temperature under the cloud of FY-3D MERSI-II.

[0004] The spatial distribution of clouds on remote sensing images is related to the solar altitude angle and the satellite observation angle (Wang Fei, 2017; WANG et al., 2017). When the solar incidence angle is inconsistent with the satellite observation angle, there is a position deviation between the clouds observed by the satellite and the clouds at the actual position (SHI et al., 2017). The larger the relative angle between the solar incidence direction and the sensor observation direction, the greater the deviation. Specifically, in practical applications, remote sensing data obtains cloud positions through cloud detection to determine the distribution of clear sky and clouds and determine whether the pixel is suitable for clear sky or cloudy weather algorithms. However, the cloud pixels detected by remote sensing images are only the projection of the real cloud cluster in the sensor observation direction, and do not correspond to the orthographic projection position of the cloud cluster on the surface, nor do they completely correspond to the area of ​​influence of the cloud on the surface radiation. The actual surface area affected by the cloud corresponds to the ground projection area formed by the entire cloud cluster in the direction of the solar incidence angle, of which the visible part is the shadow area on the remote sensing image. This effect of the solar-cloud-satellite geometry on the surface parameters is called the SCSG (Solar cloud satellite geometry) effect (WANG et al., 2017; WANG et al., 2019). At present, these effects are often ignored in the invention of surface temperature estimation under clouds, and the solar-cloud-satellite geometry effect (SCSG) of surface temperature and radiation flux has not received enough attention.

[0005] Based on the above invention background, this invention will estimate the surface temperature under the cloud for FY-3D MERSI-II data. First, the clear sky surface temperature is inverted based on the Two Factor Split Window Algorithm (TFSWA), and the sun-cloud-satellite geometric effect of surface temperature and radiation flux is considered. Based on the integrated learning model, an algorithm for estimating the surface temperature under the cloud suitable for FY-3D MERSI-II data is developed. Finally, the FY-3D MERSI-II all-weather surface temperature data is constructed to promote the further development of thermal infrared remote sensing and enhance the comprehensive benefits and service capabilities of domestic satellites.

[0006] In general, the method for estimating the surface temperature under clouds developed in the present invention has good applicability to different surface properties and climate zones, and can provide a method reference for estimating the surface temperature under clouds from polar-orbiting satellite data. Summary of the invention

[0007] The present invention belongs to the technical field of remote sensing application, and discloses a new method for estimating the surface temperature under the cloud based on the sun-cloud-satellite geometry mechanism. For all-weather remote sensing estimation of the surface temperature, the pixels are divided into clear sky pixels and cloud pixels, and the cloud pixels are further divided into cloud-blocked pixels and cloud-covered pixels. Considering the sun-cloud-satellite geometry effect of the surface temperature, the cloud pixels on the remote sensing image are divided into two types: cloud coverage (the incident solar radiation on the surface is affected by clouds, area B) and cloud blocking (the incident solar radiation on the surface is not affected by clouds, area C), and the LST under the cloud is estimated respectively. Among them, area C is the direct sunlight area, and the theoretical clear sky surface temperature represents the surface temperature of the cloud pixels in this area; the cloud coverage area (area B) is the real projection of the cloud in the direction of solar observation, and the solar radiation is affected by cloud blocking. Figure 2 The diagram shows the geometric relationship between sun, cloud and satellite observation. The surface temperature under the cloud has the complexity of multi-factor coupling, and the integrated regression model has superior performance in mining the characteristics of nonlinear and non-stationary processes. The present invention provides a new method for estimating the surface temperature under the cloud based on the integrated regression learning model, which mainly includes the following steps:

[0008] Step 1: Radiometric calibration and geometric correction of remote sensing images

[0009] Radiometric calibration is a necessary step in quantitative remote sensing applications, which converts the dimensionless DN value recorded by the sensor's original observations into spectral reflectance or radiance data with actual physical meaning. Collect L1 satellite image data, calculate the on-board reflectance of the reflection channel and the brightness temperature of the thermal emission channel, and perform geometric correction.

[0010] Step 2: Clear sky surface temperature inversion

[0011] The Planck function is linearly expanded, and a two-factor split window algorithm (TFSWA) is proposed after a series of derivations. The LST is estimated using a series of linear combinations of two brightness temperature channels, the corresponding band-averaged surface emissivity, and atmospheric transmittance, as shown below:

[0012] T s =A0+A1T 24 -A2T 25 (1-1)

[0013] A0=a 24 [D 25 (1-C 24 -D 24 ) / (C 24 D 25 -C 25 D 24 )]-a 25 [D 24 (1-C 25 -D 25 ) / (C 24 D 25 -C 25 D 24 )] (1-2)

[0014] A1=1+D 24 / (C 24 D 25 -C 25 D 24 )+b 24 [D 25 (1-C 24 -D 24 ) / (C 24 D 25 -C 25 D 24 )] (1-3)

[0015] A2=D 24 / (C 24 D 25 -C 25 D 24 )+b 25 [D 24 (1-C 25 -D 25 ) / (C 24 D 25 -C 25 D 24 )] (1-4)

[0016] C 24 =ε 24 τ 24n (1-5)

[0017] D 24 =[1-τ 24n ][1+(1-ε 24 )τ 24n ] (1-6)

[0018] C 25 =ε 25 τ 25n (1-7)

[0019] D 25 =[1-τ 25n ][1+(1-ε 25 )τ 25n ] (1-8)

[0020] Where, T s is the surface temperature (K); A0, A1, A2 are intermediate parameters, T i is the brightness temperature, τ in is the atmospheric transmittance, i is the thermal infrared band, ε i is the corresponding surface emissivity in the thermal infrared band, a 24 、a 25 , b 24 , b 25 is a constant.

[0021] Step 3: Clear sky, cloud shadow, cloud occlusion and cloud coverage area detection

[0022] For estimating the temperature of cloud pixels, it is necessary to distinguish between cloud and clear sky pixels based on cloud detection products. Different methods are used to estimate the surface temperature under the cloud according to whether the surface solar radiation corresponding to the pixel under the cloud is affected by the weakening effect of clouds. It is necessary to distinguish between cloud occlusion and cloud coverage. Therefore, it is necessary to detect the specific locations corresponding to cloud shadows, cloud occlusions and cloud coverage. In this process, orthorectification of clouds is a key step to correctly characterize the geometric relationship between the sun, cloud and satellite (SHI et al., 2017). According to the geometric relationship between the sun, cloud and satellite, the actual position of the cloud can be directly determined by using the cloud height, the zenith angle and azimuth of the sun and satellite, and the position of the cloud on the remote sensing image as input data. The spatial position of the surface shadow can be derived based on the orthorectified position (WANG et al., 2017). The angle and cloud height are derived from the geometric positioning file and the cloud height product. The orthographic position of the cloud is determined using formula (1-9), and then the orthographic projection position of the corresponding cloud shadow is estimated using formula (1-10) (WANG et al., 2017; WANG et al., 2019).

[0023]

[0024] In the formula, (Xcld , Y cld ) is the orthographic projection of the cloud on the surface; (X shw , Y shw ) is the projection position of the cloud shadow on the ground; (x map ,y map ) is the cloud position in the remote sensing image, i.e., the projected position of the cloud cluster along the satellite observation direction; and H is the cloud top height above the ground, which is calculated by subtracting the ground height (DEM product) from the cloud top height (cloud height product) of MERSI-II. v and φ v is the sensor observation zenith angle and azimuth angle at the moment of satellite transit, θ s and φ s They are the zenith angle and azimuth of the sun when the sensor is observing. When calculating the offset position, the azimuth needs to be converted to between -180° and 180°. The east, west, south, and north directions are: 90°, -90°, 180°, and 0° respectively.

[0025] Step 4: Shortwave net radiation inversion and correction of observation geometry effects

[0026] The earth's surface is often covered by clouds, which is one of the main factors of global or regional weather or climate change. Clouds transmit, absorb and reflect solar radiation, and emit long-wave radiation themselves. Clouds control the solar incident radiation and atmospheric long-wave radiation. Therefore, clouds have a decisive influence on the balance of ground radiation flux, surface temperature changes and surface energy balance. JIN et al. (2000) proposed a neighboring pixel interpolation algorithm (NP) based on the surface energy balance model (SEB) to estimate the surface temperature of pixels under clouds. The basic assumption of this algorithm is that the surface temperature difference between cloud-covered pixels and similar clear sky pixels adjacent in time and space is mainly caused by the difference in solar radiation. The surface temperature of pixels under clouds is estimated by establishing a linear relationship between surface long / short-wave radiation, latent / sensible heat flux and surface temperature. The present invention uses short-wave net radiation as the radiation characteristic variable for estimating the surface temperature under clouds.

[0027] (1) Shortwave net radiation inversion

[0028] The present invention adopts a parameterized algorithm to invert shortwave net radiation. During the calculation process, only the solar zenith angle, atmospheric water vapor, relative azimuth, atmospheric water vapor content and albedo of the top of the atmosphere need to be input as input data to invert the shortwave net radiation. This method has the advantages of simple calculation and low dependence on cloud inversion.

[0029] The shortwave net radiation inversion first converts the narrowband reflectivity ρ iConvert it to TOA short-band broadband albedo γ, and then use the parameterization method to obtain the empirical coefficient to calculate NSSR. Based on the MODTRAN model, a linear conversion formula between narrow-band reflectivity and broadband albedo can be established:

[0030] γ=b0+b1ρ1+b2ρ2+b3ρ3+b4ρ4+b5ρ5+b6ρ6+b7ρ7 (1-11)

[0031] Where γ is the satellite’s shortwave broadband albedo, ρ i is the on-board narrow-band reflectivity of MERSI-II band i, b i is the fitting parameter corresponding to MERSI-II band i. The variation of albedo with azimuth angle is smaller than that with zenith angle, so three rough intervals ranging from 0° to 60°, 60° to 120°, and 120° to 180° are selected to represent the approximate reflections in the forward, side, and reverse directions, respectively.

[0032] For a given pair of observation zenith angle (VZA) and solar zenith angle (SZA), use formula (5-20) to establish the relationship between the obtained coefficients b0~b7 and VZA

[0033]

[0034] In the formula, c 1i -c 4i are constants for given SZA = 0°, 10°, 20°, 30°, 40°, 50°, 60° and 70° respectively. Conversion parameters for different SZAs, VZAs and relative azimuth angles used here. Using MODTRAN5.2 with typical atmospheric conditions and surface boundary parameters as input, a simulation database is established through the radiation transfer model. Based on the simulation database, a parameter lookup table corresponding to different observation zenith angles, solar altitude angles and relative azimuth angles under clear sky and cloud conditions is established, and the wide-band albedo on the satellite is obtained by linear interpolation through searching the lookup table.

[0035] The shortwave net radiation is expressed as follows:

[0036] NSSR=a s E0μ / d 2 (1-13)

[0037] Where E0 is the solar irradiance on the satellite (1357W / m 2 ), d is the distance between the Earth and the Sun. Both parameters are expressed in astronomical units, a s is the surface absorption coefficient.

[0038] Surface absorption coefficient a s It can be expressed as the equation of the wide band albedo r of the top of the atmosphere:

[0039] a s (μ,ω,γ)=α′-β′γ (1-14)

[0040] The intercept α′ and slope β′ in equation (1-14) are respectively expressed as:

[0041] α′=1-α1μ -1 -a2μ -x -[1-exp(-μ)](a3+a4ω y )μ -1 (1-15)

[0042]

[0043] Where a1-a7 are constants for the various surfaces, μ is the cosine of SZA, and the values ​​of x, y, and z are 0.5.

[0044] (2) Correction of geometric effect of shortwave net radiation

[0045] The shortwave net radiation is a key input parameter for estimating the surface temperature under clouds based on the surface energy balance model. The geometric structure of the sun-cloud-satellite has a universal effect on surface radiation. The shortwave net radiation data is affected by the SCSG geometric effect and has obvious overestimation or underestimation. The distribution of clouds is offset to a certain extent on remote sensing images due to the SCSG effect. The steps for correcting the NSSR geometric effect are as follows: First, the cloud and its shadow range are determined pixel by pixel mainly using equations 1-9 and 1-10, and the position relationship between each cloud pixel and its shadow projection is calibrated. The NSSR of the shadow area (area A) pixels and the cloud covered area (area B) cloud pixels are assigned to their corresponding cloud projection area NSSR. First, for all cloud covered area pixels in the image, the above process is iteratively performed pixel by pixel. For the cloud-blocked area (area C), these pixels can be illuminated by the sun, so the NSSR of these pixels should be very similar to that of the clear sky pixels in the close distance. In this paper, the NSSR of these cloud-blocked pixels is filled by interpolation calculation of neighboring clear sky pixels or HANTS algorithm, and finally the shortwave radiation data NSSR after geometric effect correction is generated. corrected-cld .

[0046] Step 5: Feature variable selection and preprocessing

[0047] (1) Feature variable selection

[0048] The temperature characteristic factor is the necessary input data for training data and regression reconstruction of cloud pixel LST, which is used to realize the entire reconstruction process. The surface temperature is synergistically affected by geographical location, vegetation, topography, micro-topography, surface effective radiation, meteorological information, etc.

[0049] The present invention selects characteristic variables for the regression operation of the surface temperature under the cloud from the aspects of surface properties, solar radiation, atmospheric conditions and sensors. Among all the factors related to the surface temperature, elevation is the preferred important influencing factor. The correlation between altitude and temperature is a quantifiable parameter. The influence of elevation on the change of surface temperature is particularly significant in plateaus and mountainous areas. When the altitude decreases, the influence on the distribution of radiation energy is weakened. DEM, slope and aspect are selected as surface terrain characterization parameters. Vegetation selectively absorbs and reflects solar radiation, regulates latent heat and sensible heat exchange to affect the surface thermal radiation information. Vegetation parameters are indispensable, and NVI, EVI and LAI are selected as vegetation parameters. In addition, the daily variation curve of LST is related to the short-wave net radiation flux of surface absorption and heat generation, and the close relationship between the daily surface temperature cycle and the net solar radiation (short-wave net radiation) of the surface has also been determined. Therefore, in terms of solar radiation, the surface short-wave net radiation and the solar altitude angle are considered. According to the shortwave net radiation inversion process, the input data of the surface shortwave net radiation inversion include solar constant, solar altitude and azimuth, atmospheric water vapor (WVC) and sensor observation angle, which shows that this characteristic variable is a comprehensive representation of geographical location, solar radiation, sensor observation conditions and atmospheric conditions. Table 1 gives the time resolution, unit and other detailed information of each characteristic variable.

[0050] Table 1 Characteristic variables for estimating surface temperature under clouds

[0051]

[0052] (2) Harmonic optimization processing of NVI, EVI and LAI data

[0053] NVI, EVI and LAI are important variables that indicate the characteristics of surface vegetation cover. Due to the influence of atmospheric conditions, aerosols and sensor failures during remote sensing data observation, there are certain data missing and low-quality data in the MERSI-II LAI data product, which seriously affects the model training application in the estimation of surface temperature under clouds. In order to obtain higher-precision model simulation results, the HANTS algorithm is used to fill the missing pixels of NVI, EVI and LAI data, smooth and filter the low-quality data caused by cloud pollution or other factors, and generate continuous and complete long-term series of high-quality training data.

[0054] (3) Collinearity analysis of characteristic variables

[0055] In order to select high-quality feature variables, it is necessary to check whether there is collinearity between variables. Based on the Seaborn statistical analysis visualization module, the PEARSON feature correlation matrix heat map between feature variables is calculated, and variables with large correlation are removed.

[0056] Step 6: Model building and parameter optimization

[0057] (1) Determination of training data set

[0058] In order to accurately represent the complex relationship between LST and characteristic variables, a fitting regression model from each characteristic variable to the surface temperature was constructed. High-quality clear sky and shadow data were selected as training data sets, and the surface temperature of cloud pixels was used as the target variable to establish a regression model. Based on Python's third-party modules Scikit-learn and XGBoost, three regression models, RF, GBRT and XGBoost, were used to train the surface temperature data sets of clear sky and shadow pixels. The train_test_split function was used to randomly extract and divide the sample data into two parts, 80% training partitions and 20% test partitions, which were applied to the training model and the test model respectively.

[0059] (2) Parameter setting and optimization

[0060] To prevent overfitting, the number of model regression trees (max-depth) in the three models was adjusted, and the scores of the training data were checked through 10-fold cross validation. To check the performance of the model, the training time was checked for different training samples.

[0061] Step 7: Train the clear sky and shadow area based on the ensemble learning model to establish a regression prediction model between the feature variables and LST

[0062] Considering the multi-factor coupling characteristics of the surface temperature under the cloud, the integrated regression algorithm has excellent performance in feature mining of nonlinear transformation processes and non-stationary characteristics. On the basis of fully considering the sun-cloud-satellite geometric effect of surface temperature and radiation, three integrated regression algorithms, RF, GBRT and XGBoost, are used to establish the LST regression prediction model for the clear sky area and the shadow area with the target variable (FY-3D MERSI-II surface temperature) and the pre-processed feature variables through feature training simulation. At the same time, based on the Python third-party extension package scikit-learn feature variable importance calculation method (PermutationImportance), the importance score and ranking of each feature variable are calculated.

[0063] Step 8: Input the image to be calculated and use the regression model to estimate the LST under the cloud in the cloud-covered and cloud-covered areas

[0064] The sub-cloud LST regression model established in the previous step is used to predict the surface temperature of all pixels in the cloud-blocked and cloud-covered areas in the invention area, and the CLDAS assimilated LST data are used to evaluate the accuracy of the estimated sub-cloud surface temperature.

[0065] Compared with other methods for estimating surface temperature under clouds, its advantages are:

[0066] 1. The integrated regression algorithm is used to estimate the surface temperature under the cloud, which has a faster calculation speed and improves the efficiency of data processing.

[0067] 2. The algorithm is not affected by factors such as the underlying surface type and cloud coverage size, which improves the applicability of cloud pixel LST estimation.

[0068] 3. The integrated regression model of three decision tree types is used to estimate the surface temperature under the cloud, which can improve the accuracy of model fitting and prevent overfitting through parameter optimization. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 Flowchart of the integrated regression algorithm for subcloud surface temperature considering surface temperature and radiation SCSG effects;

[0070] Figure 2 Schematic diagram of the geometric relationship between sun-cloud-satellite observations;

[0071] Figure 3 FY-3D MERSI-II cloud shadow detection and pixel statistics of cloud occlusion and cloud coverage areas;

[0072] Figure 4 Pixel classification results of SCSG effect in FY-3D MERSI-II in four seasons;

[0073] Figure 5 Schematic diagram of NSSR geometric effect correction;

[0074] Figure 6 Shortwave net radiation data after correction for geometric effects;

[0075] Figure 7 Pearson feature correlation matrix heat map of machine learning regression model;

[0076] Figure 8 Distribution diagram of input characteristic variables and LST density in clear sky and shadow area on May 2, 2020;

[0077] Fig. 9 Maximum tree depth and model validation set score;

[0078] Fig.10 The time required for the three models to train different numbers of samples;

[0079] Fig.11 Training accuracy and time taken;

[0080] Fig.12 Importance of features of surface temperature regression operation under cloud for RF, GBRT and Xgboost algorithms;

[0081] Fig.13 FY-3D MERSI-II surface temperature image under clouds based on three integrated regression models;

[0082] Fig.14 Scatter plot of subcloud surface temperature and CLDAS temperature based on FY-3D MERSI-II data. DETAILED DESCRIPTION

[0083] In order to implement the proposed method, we take the FY-3D medium resolution imager MERSI-II as an example to conduct an experiment on estimating the surface temperature under clouds. We need to perform geometric registration between the target image and the reference image and determine whether the pixel is clear sky or cloud covered. The process has three parts: radiometric calibration and geometric correction of the reference image, clear sky surface temperature inversion, cloud detection to identify cloud occlusion, cloud covered pixels, shadow pixels and clear sky pixels, selection and preprocessing of feature variables, model establishment and training, and estimation of surface temperature under clouds.

[0084] FY-3D is the fourth satellite of the FY-3 series of my country's second-generation polar-orbiting meteorological satellites. The Medium Resolution Spectral Imager II (MERSI-II) is one of the core instruments on FY-3D. It integrates the functions of the original Fengyun-3 satellite imaging instruments MERSI and VIRR. It has a total of 25 channels, of which channels 1-19 are solar reflection channels (0.4-2.1μm), 20-25 are infrared emission channels (3.8-12.5μm), channels 1-4, 24-25 have a spatial resolution of 250 meters, and channels 5-23 have a spatial resolution of 1000m. FY-3D MERSI-II currently mainly provides L1, L2 and L3 data products, a total of 25 data sets, and the data storage format is the HDF5 scientific data set format.

[0085] The present invention takes the sun-cloud-satellite geometric effect as the starting point and develops the FY-3D MERSI-II under-cloud surface temperature integrated regression algorithm. Firstly, starting from the influence of solar radiation, surface and atmospheric conditions on the surface temperature, 7 variables are selected as the surface thermal radiation correlation factors, and the key characteristic variable shortwave net radiation is corrected by the sun-cloud-satellite geometric effect, and each characteristic variable is preprocessed by collinearity analysis. Three integrated regression models, Random Forest (RF), Gradient Boosting Tree (GBRT) and Autoregressive Extreme Gradient Boosting (XGBoost), are used. The FY-3D MERSI-II clear sky and cloud shadow area data are used as training data sets, and the surface temperature of cloud occlusion and cloud coverage area is simulated and predicted as the target variable. A set of cloud LST integrated regression algorithms that fully consider the surface temperature and radiation flux sun-cloud-satellite geometric effects are established. Based on the algorithm, four days of data from different seasons in North China in 2020 are used to predict the surface temperature under the cloud. Based on FY-3D MERSI-II, the surface temperature estimation of clear sky and under the cloud requires MERSI-II geometric positioning data and some data products (cloud height, cloud detection, NVI and LAI ten-day products) as auxiliary data. The present invention is further described in detail in conjunction with the accompanying drawings. Table 2 shows the FY-3D MERSI-II dataset used for all-weather surface temperature estimation.

[0086] Table 2 FY-3D MERSI-II dataset used for all-weather surface temperature estimation

[0087]

[0088] Step 1: FY-3D MERSI-II data preprocessing

[0089] According to the selected scope of the North China Invention Area, the L1-level satellite image data of FY-3D MERSI-II is collected. This part of the work mainly includes three parts: on-board reflectivity calculation of the reflection channel, brightness temperature calculation of the thermal emission channel, and geometric correction.

[0090] 1. On-board reflectivity calculation

[0091] MERSI-II L1 data uses unsigned integer data to store the original observed pixels to reduce the data size, and the grayscale value is output after multi-detector normalization correction. Radiometric calibration is a necessary step in quantitative remote sensing applications, which converts the dimensionless DN value recorded by the sensor's original observations into spectral reflectance or radiance data with actual physical meaning.

[0092] To calculate the on-board reflectivity of the solar reflection band CH1 to 19, the DN value must first be converted into the visible light calibration reflectivity Ref. The data sets EV_1KM_RefSB and EV_250_Aggr.1KM_RefSB are the normalized DN values ​​of the 1km reflection channel (CH5 to 19) and the 250m reflection channel (CH1 to 4), respectively. The calibration reflectivity conversion formula is as follows:

[0093] dn=DN*Slope+Intercept (1-17)

[0094] Ref=Cal2*dn 2 +Cal1*dn+Cal0 (1-18)

[0095] Where Cal0, Cal1 and Cal2 are the calibration coefficients of the corresponding channels in the dataset VIS_Cal_Coeff (corresponding to the 1st, 2nd and 3rd columns respectively), and Slope and Intercept are the attributes of the datasets EV_1KM_RefSB and EV_250_Aggr.1KM_RefSB.

[0096] Then the visible light calibration reflectance Ref is converted into the entrance pupil radiance L toa :

[0097]

[0098] Combined with the entrance pupil radiance L toa It can be converted into apparent reflectance ρ toa :

[0099]

[0100] Where E0 is the solar irradiance outside the atmosphere of the channel, which can be obtained from the attribute Solar_Irradiance in the L1 dataset, and DES is taken from the attribute EarthSun Distance Ratio in the L1 dataset. μ is cosθ, and θ is the solar zenith angle, which can be obtained from the SolarZenith attribute of the MERSI-II L1 scientific dataset. Since the storage method is scaled integer data, the actual solar zenith angle θ needs to be calculated using the corresponding slope and intercept coefficients according to the following formula.

[0101] θ=SolarZenith×Slope+Intercept (1-21)

[0102] 2. Brightness temperature calculation

[0103] FY-3D MERSI-II data EV_1KM_Emissive and EV_250_Aggr.1KM_Emissive are the amplified radiance values ​​of the 1km transmission channel (CH20-23) and the 250m transmission channel (CH24-25), respectively. The gain (scales) and offset (offsets) of each band i are obtained from the attributes, and the radiance I of band i is calculated according to the following formula i (W·m -2 ·sr -1 μm -1 )(Xu Na et al., 2019):

[0104] I i =(DN-offset) scales (1-22)

[0105] The radiance of the thermal infrared channel (CH24-25) is converted to the blackbody brightness temperature. The equivalent blackbody brightness temperature T is calculated by the inverse transformation of the Planck function through the equivalent center wavelength effect_Center_WaveLength (unit: μm) and the channel radiance. ei , the calculation formula is as follows:

[0106]

[0107] In the formula, λ i is the effective center wavelength of the i-th (i=24, 25) band, which can be obtained according to the attribute variable effect_Center_WaveLength (unit μm) in the L1 file; From Table 3-1, it can be seen that the center wavelength λ of the 24th and 25th channels of MERSI-II i Take λ respectively 24 =11.8μm and λ 25 =12.0μm; C1 and C2 are the first and second spectral constants, respectively, and C1=1.19104356×10 -16 W / m 2 and C2=1.4387685×10 4 μm·K. Finally, the channel brightness temperature correction coefficient (TbbCorr_Coeff) needs to be used to convert T ei Convert to channel blackbody brightness temperature T i :

[0108] T i =A*T ei +B (1-24)

[0109] Where, T i is the brightness temperature of the mersi-ii band i (i = 24, 25), i.e., the brightness temperature T corresponding to thermal infrared channels 24 and 2524 and T 25 ; A and B are the attribute variables TBB_Trans_Coefficient_A and TBB_Trans_Coefficient_B in the L1 file.

[0110] 3. Geometric correction

[0111] FY-3D MERSI-II data adopts multi-detector scanning technology, and the geometric characteristics of sensor imaging, terrain undulation, surface curvature and satellite jitter lead to the generation of image geometric distortion. The present invention adopts the Geographic Lookup Table (GLT) method to perform system positioning. Specifically, the geometric positioning information of the corresponding data input by the geometric positioning file is used to generate a geographical location lookup table. The actual geographical location of each pixel in the original data can be inferred through the lookup table. However, due to the influence of many factors such as terrain, surface and atmospheric refraction, this preliminary system correction result still has a geometric positioning error of 2-10 pixels, which cannot meet the application requirements. To this end, the present invention adopts the remote sensing image automatic matching technology based on the Moravec corner detection algorithm to perform geometric precision correction on the preliminary corrected data. The algorithm effectively reduces the local geometric distortion of the image through automatic control point selection, improves the geometric positioning accuracy of the data, and the error is basically controlled within 3 pixels. The geometric correction accuracy can basically meet the requirements. Accurate geometric positioning lays a good geometric foundation for the smooth development of the subsequent all-weather surface temperature estimation work of FY-3D MERSI-II data. Accurate geometric positioning has laid a good geometric foundation for the smooth implementation of the subsequent all-weather surface temperature estimation work of FY-3D MERSI-II data. Prepare NDVI, LAI, DEM, slope, aspect, solar altitude angle, satellite observation angle and other data corresponding to the area and time of the data to ensure the consistency of data projection and resolution.

[0112] Step 2: FY-3D MERSI-II clear sky surface temperature inversion

[0113] The Planck function is linearly expanded, and a two-factor split window algorithm (TFSWA) is proposed after a series of derivations. A series of linear combinations of two brightness temperature channels, the corresponding band-averaged surface emissivity and atmospheric transmittance are used to estimate LST. Between -50° and 50°, the constant term a of the MERSI-II TFSWA algorithm band 24 24 =-53.477, b 24 =0.3951; the constant term of band 25 is a 25 =-57.087, b 25 =0.4292.

[0114] Step 3: Cloud shadow, cloud occlusion and cloud coverage area detection

[0115] MERSI-II cloud detection products mainly provide MERSI-II satellite data cloud and clear sky distribution information, which is a basic auxiliary data for the study of surface temperature estimation under clouds. The present invention uses the cloud detection results obtained by comprehensive analysis of the FY-3D MERSI-II cloud detection product by the multi-band (3, 4, 6, 7, 19, 20, 21, 24, 25) feature threshold method. According to whether the surface solar radiation corresponding to the pixel under the cloud is affected by the weakening effect of the cloud, different methods are used to estimate the surface temperature under the cloud, and it is necessary to detect the specific positions corresponding to the cloud shadow, cloud occlusion and cloud coverage area. Therefore, combined with the angle and cloud height data of FY-3D MERSI-II, the method in step three is used to carry out cloud occlusion, cloud coverage and cloud shadow detection for a total of 4 days of data in different seasons in 2020 in North China. Figure 3 The statistical results of cloud shadow, cloud occlusion, cloud cover and clear sky pixels of FY-3D MERSI-II are shown in Figure 2. Figure 4 The pixel classification results of four images in different seasons in 2020 are shown.

[0116] Step 4: Shortwave net radiation inversion and geometric effect correction

[0117] Shortwave net radiation is an important characteristic parameter in the inversion of surface temperature under clouds. This time, a parameterized algorithm was used to invert the FY-3D MERSI-II shortwave net radiation, and the shortwave net radiation was corrected for SCSG geometric effects. MERSI-II has two sets of visible light bands, of which bands 8 to 14 have a narrow bandwidth (20nm) and are mainly used for inventions such as ocean color and plankton. Bands 1 to 7 are designed for land applications. The inversion of the albedo on the MERSI-II satellite mainly uses these seven reflected solar radiation (RSB) channels and one water vapor absorption channel (Table 3).

[0118] Table 3 MERSI-II reflected solar radiation bands used for NSSR inversion

[0119]

[0120] The shortwave net radiation inversion first converts the narrowband reflectivity ρ i Convert it to TOA short-band broadband albedo γ, and then use the parameterization method to obtain the empirical coefficient to calculate NSSR. Based on the MODTRAN model, a linear conversion formula between narrow-band reflectivity and broadband albedo can be established:

[0121] γ=b0+b1ρ1+b2ρ2+b3ρ3+b4ρ4+b5ρ5+b6ρ6+b7ρ7+b 19 ρ 19 (1-25)

[0122] Where γ is the satellite’s shortwave broadband albedo, ρ i is the on-board narrow-band reflectivity of MERSI-II band i, b i is the fitting parameter corresponding to MERSI-II band i. The variation of albedo with azimuth angle is smaller than that with zenith angle, so three rough intervals ranging from 0° to 60°, 60° to 120°, and 120° to 180° are selected to represent the approximate reflections in the forward, side, and reverse directions, respectively.

[0123] Table 4 Main parameter settings in MODTRAN5.2 atmospheric radiation transfer simulation

[0124]

[0125]

[0126] According to the parameterization process in step 4, the FY-3D MERSI-II shortwave net radiation can be estimated.

[0127] Step 5: Feature variable selection and preprocessing

[0128] (1) Feature variable selection

[0129] The MERSI-II LAI leaf area index 10-day product and the NVI vegetation index 10-day product data are mainly used for similar pixel search and as feature variables of the integrated regression algorithm. The two data sets are the GLL equal longitude and latitude projection and the Hammer projection, and are stored in a global 10°×10° frame with a spatial resolution of 0.01°.

[0130] MERSI-II vegetation index ten-day product data, including NDVI, EVI, sensor and satellite angle information, and quality control information. MERSI-II leaf area index ten-day product mainly includes leaf area index and corresponding data quality control information. Cloud top height refers to the height of the cloud top, in hPa or m. All data are uniformly converted to geometric height units for calculating cloud shadow positions and geometric effect correction of radiation data. All auxiliary data are unified into WGS84 geographic coordinate system, 0.01° resolution and other longitude and latitude projections.

[0131] (2) Variable preprocessing

[0132] a. Correction of geometric effect of shortwave net radiation

[0133] Considering the universality of the SCSG effect of radiation data, the shortwave net radiation, a key characteristic variable used to estimate the surface temperature under clouds, is corrected for the SCSG geometric effect. For specific methods, refer to step 4. Figure 5A schematic diagram of the geometric effect of shortwave radiation data is shown, which shows the MERSI-II NSSR observed over North China on May 2, 2020. Figure 5 (a) is an RGB true color image, and the two straight lines m and n are used to invent the interpolation profile curves before and after NSSR correction. Figure 5 (b) is the classification result of each type of pixel: cloud shadow area in area A, cloud covered area in area B, cloud blocked area in area C, and clear sky area in area D. In order to better characterize the influence of the sun-cloud-satellite geometry on NSSR, the NSSR corrected by SCSG geometry effect is subtracted from the original NSSR data. Figure 5 (c). From the difference map, we can see that the shadow area (Area A) in the image is negative, and the cloud-blocked area (Area C) is positive, which means that if the geometric effect of the sun-cloud-satellite is ignored, the NSSR will be underestimated in the cloud-blocked area (Area C) and overestimated in the shadow area (Area A). The specific size of this difference varies with the physical and optical characteristics of the cloud. In the example given, this difference can reach about 800W / m 2 , accounting for about 80% of the clear-sky NSSR. Figure 5 (d) and (e) are the NSSR differences between the two lines m and n (NSSR corrected-eld -NSSR clear ) profile. The dashed line on the left corresponds to the cloud-blocked area of ​​the original data, which is an underestimated area. The dashed line on the right corresponds to the shadow area of ​​the original data, which is an overestimated area of ​​the NSSR. Figure 6 The images of shortwave net radiation data after geometric effect correction at different times in four seasons (20200218, 20200502, 20200715 and 20201113) are given, as well as the corresponding cloud detection results, i.e., pixel classification results.

[0134] b. Harmonic optimization processing of NVI, EVI and LAI data

[0135] In order to obtain higher-precision model simulation results, the HANTS algorithm is used to fill in the missing pixels of the FY-3D MERSI-II NVI, EVI and LAI data, and the low-quality data caused by cloud pollution or other factors are smoothed and filtered to generate a continuous and complete long-term series of high-quality training data.

[0136] c. Collinearity analysis of characteristic variables

[0137] In order to select high-quality feature variables, it is necessary to check whether there is collinearity between the variables. Based on the Seaborn statistical analysis visualization module, the PEARSON feature correlation matrix heat map between the above 9 feature variables is calculated ( Figure 7), it can be seen that the correlation between SZA and NSSR is above 0.7, and there is a collinearity problem, so the SZA variable is removed. There is no obvious collinearity problem with other variables. In addition, the feature importance analysis was performed, and it was found that the aspect score was the lowest. Therefore, in order to reduce the impact of noise, the aspect feature variable was removed, and 7 valid feature variables were determined. On May 2, 2020, there were 802,468 clear and shadow pixels; these pixels were widely distributed covering different land cover types, vegetation coverage, terrain characteristics, and solar radiation conditions. Based on the observation data on May 2, 2020, the relationship between the characteristic variables and temperature in the clear sky and cloud shadow area LST, that is, the remote sensing visible area, was analyzed. Figure 8 The density distribution of the seven characteristic variables and LST determined in the clear sky and shadow area on May 2, 2020 is shown. It can be seen from the figure that the distribution of surface temperature and various parameters is relatively concentrated, and there is no obvious linear relationship between the surface temperature under the cloud and the parameters. This complex relationship confirms that the surface temperature is jointly affected by many surface variables, and the coupling effect between various factors leads to this complexity.

[0138] Step 6: Model building and parameter optimization

[0139] (1) Determination of training data set

[0140] In order to accurately represent the complex relationship between LST and characteristic variables, a fitting regression model from each characteristic variable to the surface temperature was constructed. High-quality clear sky and shadow data on May 2, 2020 were selected as the training data set, and the surface temperature of cloud pixels was used as the target variable to establish a regression model. Based on Python's third-party modules scikit-learn and XGBoost library, three regression models, RF, GBRT and XGBoost, were used to train the surface temperature data set of clear sky and shadow pixels. In this process, the train_test_split function was used to randomly extract and divide the sample data into two parts, 80% training partition and 20% test partition, which were applied to the training model and the test model respectively.

[0141] (2) Parameter setting and optimization

[0142] In order to prevent overfitting, the number of model regression trees (max-depth) parameters in the three models were adjusted, and the scores of the training data were checked through 10-fold cross validation. The results are as follows Fig. 9 As max-depth increases from 0 to 20, the validation score of the model increases. When max-depth = 8, 5, and 4 are used for RF, GBRT, and XGBoost, the model scores reach the extreme values. max_depth = 8, 5, and 4 are selected as the optimal parameters for RF, GBRT, and XGBoost. In order to check the performance of the model, the training time ( Fig.10 ), for the same number of training samples, GBRT takes the shortest time, while XGBoost takes the longest time. Fig.11 The figure shows the time taken by each model to achieve different accuracies. From the analysis of the figure, we can see that the three integrated tree methods, random forest and XGBoosting algorithms, consume relatively long time but achieve relatively high accuracy.

[0143] Step 7: Regression calculation of surface temperature under clouds

[0144] For the central and eastern regions of Asia, several images were selected in four seasons of 2020 to conduct surface temperature reconstruction experiments under clouds and evaluate the proposed method. Based on the feature variable importance calculation method (Permutation Importance) of the Python third-party extension package scikit-learn, the importance score and ranking of each feature variable were calculated. Fig.11 The calculation results of the feature importance of three regression models are given. As can be seen from the figure, among the parameters used by the three models, the top three feature variables are consistent. The parameter with the greatest impact on the surface temperature under the cloud is DEM, and the feature importance ranges from 0.87 of XGBoost to 0.95 of RF, followed by NWDI (0.76-0.82) and NSSR (0.35-0.39), and finally the vegetation index and leaf area index. Different DEMs will change the optical path length of solar radiation. In addition, the thickness of the atmosphere that the sun's rays pass through when reaching the ground, the higher the altitude, the thinner the atmosphere, and the stronger the radiation. The elevation condition determines the energy distribution, so it has a greater impact on the surface thermal radiation in the local area. Shortwave net radiation is a comprehensive representation of solar radiation, clouds and the surface, and is one of the dominant factors in the change of surface temperature. In general, in the invention area, the condition that has the greatest impact on the surface temperature under the cloud is the terrain condition, followed by NDWI and surface shortwave net radiation, and finally the vegetation coverage condition. Fig.12 FY-3DMERSI-II image of surface temperature distribution under clouds estimated based on three machine learning regression model algorithms.

[0145] The accuracy of the algorithm was verified using the surface temperature data from the China Meteorological Administration Land Data Assimilation System (CLDAS). Fig.14 The scatter plot of the subcloud surface temperature and CLDAS temperature based on FY-3D MERSI-II data is shown. Compared with the CLDAS surface temperature data, the results show that the RMSE accuracy of this method is about 2.84K-5.43K. The results show that the subcloud surface temperature integrated regression algorithm considering the geometric effect of sun-cloud-satellite observation has a good accuracy in estimating the surface temperature, which can better describe the impact of terrain and vegetation on the surface temperature.

Claims

1. A method for estimating surface temperature under clouds based on sun-cloud-satellite observation geometry, characterized in that: The following steps are involved: S1: Collect and preprocess the L1 satellite image data of FY-3D MERSI-II; S2: clear sky surface temperature inversion; S3: cloud shadow, cloud occlusion and cloud coverage area detection; S4: Use parameterized algorithm to invert the FY-3D MERSI-II shortwave net radiation and perform SCSG geometric effect correction on the shortwave net radiation; The step S4 specifically comprises the following steps: The shortwave net radiation inversion first converts the narrowband reflectivity ρ i Convert it to TOA short-band broadband albedo γ, and then use the parameterization method to obtain the empirical coefficient to calculate the NSSR; based on the MODTRAN model, a linear conversion formula between narrow-band reflectivity and broadband albedo can be established: γ=b0+b1ρ1+b2ρ2+b3ρ3+b4ρ4+b5ρ5+b6ρ6+b7ρ7+b 19 r 19 (1-11) Where γ is the satellite’s shortwave broadband albedo, ρ i is the on-board narrow-band reflectivity of MERSI-II band i, b i is the fitting parameter corresponding to MERSI-II band i; Then, the cloud and its shadow range are determined pixel by pixel, and the positional relationship between each cloud pixel and its shadow projection is calibrated. The NSSR of the cloud-blocked area pixel is generated by interpolation calculation using the Hants algorithm. The NSSR of the cloud-covered area is geometrically corrected, and finally the shortwave radiation data corrected for geometric effects are generated. S5: Select characteristic variables of the subcloud surface temperature regression operation and preprocess them; S6: Parameter setting and optimization for three integrated regression models; The step S6 specifically includes: Three regression models, RF, GBRT and XGBoost, are used to train the temperature data sets of clear sky and shadow pixels. The train_test_split function is used to randomly extract sample data and divide them into two parts, training partition and test partition, which are applied to the training model and the test model respectively. S7: For the LST in clear sky and shadow area, an integrated regression model of each characteristic variable to the surface temperature under the cloud is established; The step S7 specifically includes: Three integrated regression algorithms, RF, GBRT and XGBoost, were used to perform regression operations with FY-3D MERSI-II subcloud surface temperature as the target variable, and the accuracy of the estimated subcloud surface temperature was evaluated using CLDAS assimilated LST data. S8: Input the image to be calculated and use the established regression model to estimate the LST of cloud occlusion and cloud coverage areas.

2. The method for estimating surface temperature under clouds based on sun-cloud-satellite observation geometry according to claim 1, characterized in that: The step S1 specifically includes the following steps: S11: Calculate the on-star reflectivity; S12: Calculate brightness temperature; S13: The system positioning work is carried out based on the Geographic Lookup Table (GLT) method, and the Moravec corner detection algorithm is used to perform geometric precision correction on the preliminary corrected data.

3. The method for estimating surface temperature under clouds based on sun-cloud-satellite observation geometry according to claim 1, characterized in that: The step S2 specifically includes the following steps: LST is estimated using a series of equations from two brightness temperature channels, the corresponding band-averaged surface emissivity, and atmospheric transmittance.

4. The method for estimating surface temperature under clouds based on sun-cloud-satellite observation geometry according to claim 1, characterized in that: The step S3 specifically comprises the following steps: Determine the orthographic position of the cloud: Then estimate the orthographic projection position corresponding to the cloud shadow: In the formula, (X cld , Y cld ) is the orthographic projection of the cloud on the surface; (X shw , Y shw ) is the projection position of the cloud shadow on the ground; (x map ,y map ) is the cloud position in the remote sensing image, i.e., the projected position of the cloud cluster along the satellite observation direction; and H is the cloud top height above the ground, which is calculated by subtracting the ground height from the cloud top height of MERSI-II; θ v and is the sensor observation zenith angle and azimuth angle at the moment of satellite transit, θ s and are the zenith angle and azimuth angle of the sun when the sensor observes.

5. According to the method for estimating the surface temperature under the clouds based on the sun-cloud-satellite observation geometry according to claim 1, the characteristic variables of the regression operation of the surface temperature under the clouds are specifically a combination of one or more of the following variables: leaf area index, normalized vegetation index, enhanced vegetation index, digital ground elevation, slope, aspect, normalized water index, solar altitude angle, and shortwave net radiation.

6. The method for estimating surface temperature under clouds based on sun-cloud-satellite observation geometry according to claim 1, characterized in that: The step S8 specifically includes: The constructed regression model is used to predict the surface temperature under clouds in the invention area, and the predicted cloud pixel LST is merged with the clear sky LST to obtain a seamless all-weather LST image of the invention area.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 6 is implemented.

Images