Evaluation of Three Parametric Models for Estimating Directional Thermal Radiation from Simulation, Airborne, and Satellite Data

Abstract

An appropriate model to correct thermal radiation anisotropy is important for the wide applications of land surface temperature (LST). This paper evaluated the performance of three published directional thermal radiation models—the Roujean–Lagouarde (RL) model, the Bidirectional Reflectance Distribution Function (BRDF) model, and the Vinnikov model—at canopy and pixel scale using simulation, airborne, and satellite data. The results at canopy scale showed that (1) the three models could describe directional anisotropy well and the Vinnikov model performed the best, especially for erectophile canopy or low leaf area index (LAI); (2) the three models reached the highest fitting accuracy when the LAI varied from 1 to 2; and (3) the capabilities of the three models were all restricted by the hotspot effect, plant height, plant spacing, and three-dimensional structure. The analysis at pixel scale indicated a consistent result that the three models presented a stable effect both on verification and validation, but the Vinnikov model had the best ability in the erectophile canopy (savannas and grassland) and low LAI (barren or sparsely vegetated) areas. Therefore, the Vinnikov model was calibrated for different land cover types to instruct the angular correction of LST. Validation with the Surface Radiation Budget Network (SURFRAD)-measured LST demonstrated that the root mean square (RMSE) of the Moderate Resolution Imaging Spectroradiometer (MODIS) LST product could be decreased by 0.89 K after angular correction. In addition, the corrected LST showed better spatial uniformity and higher angular correlation.

1. Introduction

Land surface temperature (LST) is an indispensable parameter of the land surface energy budget, widely used in climate change, evapotranspiration, urban heat island studies, and other fields [1,2]. However, there are still several difficulties in providing accurate and harmonized LST datasets for long-term Earth observation. Directionality of thermal radiation is such an obstacle, and is described as temperature differences among different viewing geometries for the same ground object and observation time [3,4]. Many experimental studies based on various platforms for various land cover types have demonstrated that the temperature differences are significant and even reach up to 15 K [5,6,7,8,9,10,11], dramatically reducing comparability of LSTs retrieved from multisensor data. Therefore, quantification and correction of thermal radiation directional effects is an urgent and meaningful undertaking.

Previous efforts have developed several models to estimate directional anisotropy. For example, Guillevic et al. extended the Discrete Anisotropy Radiative Transfer (DART) model to the thermal infrared (TIR) domain to simulate the TIR radiative budget and upward spectral radiance of heterogeneous canopies [12] and used the model to analyze directional viewing effects on satellite LST products over sparse vegetation cover [13]. Working on the basis of the four-stream radiative transfer formalism, Verhoef et al. developed the Scattering by Arbitrarily Inclined Leaves (4SAIL) model which relates top-of-canopy thermal radiance directionality with the temperatures of sunlit and shaded soil and sunlit and shaded leaves [14]. Duffour et al. demonstrated the capability of the Soil Canopy Observation of Photochemistry and Energy fluxes (SCOPE) model to simulate TIR directional anisotropy [15] and then explored the driving factors of directional variability using the model [16]. The above-mentioned models are physics-based and can be used to obtain accurate simulated results but have low computational efficiency. Moreover, they need too many spatial and temporal parameters describing the land surface properties, which are difficult to measure. All of these shortcomings greatly hinder their practical applications.

By contrast, parametric models are very attractive owing to their simplicity, robustness, and ability to be applied in any spatial scale [17]. Relative studies are still rare and learn more or less from the visible/near-infrared (VNIR) domain. Although physical meanings of radiation are somewhat different compared with in the VNIR domain, sunlit/shaded elements also have different values resulting in different temperatures in the TIR domain. Therefore, Lagouarde and Irvine adapted a reflectance hotspot model to describe the thermal infrared directional signature by replacing the reflectance with the surface temperature [8]. The model was validated using airborne measurements and the SCOPE model, and produced a hopeful result [17]. With the same methodology, Peng et al. [18] identified the most suitable kernel-driven Bidirectional Reflectance Distribution Function (BRDF) model for fitting and extrapolating directional anisotropy using the airborne multi-angle thermal infrared dataset; the model was also tested against the 4SAIL model [19]. Following the traditional structure of BRDF models, Vinnikov et al. developed an empirical model based on a linear combination of emissivity and solar kernels using one full year of simultaneous observations by two Geo-stationary Earth Orbit (GEO) satellites at five Surface Radiation Budget Network (SURFRAD) stations [4]. Currently, the model is the only one applied for angular correction of satellite-derived LSTs [20].

Field experiments can reveal thermal radiation directional features and provide valuable data for model validation. There are various relevant experimental studies including ground-based [6,7], tower-based [5,11], and airborne platforms [8,10,21]. However, each observation is only for a single canopy with a fixed structure. Anisotropy data with a variety of situations is scarce and a large simulated dataset is therefore required for extensive evaluation of parametric models [17]. At the pixel scale, current studies mainly use GEO satellites to develop multi-angle datasets because of their high temporal resolution [4,20,22]. In turn, the datasets have a very limited range of viewing angles, small observation area, and complicated spatial heterogeneity. Low Earth Orbit (LEO) satellites observe across the globe and their broader coverage may make up for the deficiency in their temporal resolution. More importantly, LEO satellites have more rich viewing geometries and land cover types. Therefore, multi-angle datasets developed using LEO satellites can provide more comprehensive validation for models.

The main objective of this paper is to evaluate the performance of three directional thermal radiation parametric models. At the canopy scale, a large anisotropy dataset generated by the 4SAIL model and an airborne multi-angle dataset acquired by the Wideangle infrared Dual-mode line/area Array Scanner (WiDAS) system are used as the reference. At the pixel scale, the first LEO satellite multi-angle LST dataset established using Moderate Resolution Imaging Spectroradiometer (MODIS) and Advanced Along-Track Scanning Radiometer (AATSR) products is applied to assess the models. Section 2 describes the data including the synthetic, airborne, satellite, and SURFRAD datasets. Section 3 shows the details of the three parametric models, while Section 4 presents the evaluations of the models. Section 5 displays the application of LST angular correction. Section 6 presents the discussion, and conclusions are drawn in the last section.

2. Data

2.1. Simulation Dataset

As a widely validated and used model, the 4SAIL model is deemed well suited for turbid canopies [23]. Because the DART model is suitable for discrete canopies and the SCOPE model—a combination of several models—uses precisely the 4SAIL model to simulate radiative transfer, we chose the 4SAIL model as the data generator. The model divides the canopy into sunlit and shaded soil and sunlit and shaded leaves to simulate top-of-canopy (TOC) thermal radiation for different viewing angles. It can also analyze the influences of leaf area index (LAI), leaf inclination distribution function (LIDF), and hotspot effect on TOC thermal radiation. For more information, readers can refer to the paper that introduced the 4SAIL model [14].

We use the same input parameters as did Verhoef et al. [14], owing to the fact that the situation was well validated with ground measurements. Different canopy structures are created by prescribing four LIDFs (erectophile, spherical, plagiophile, planophile) and four LAIs (0.5, 1, 2, 4) with a hotspot parameter of 0.05. For each of them, the viewing zenith angles (VZAs) vary from nadir to 60° and viewing azimuth angles (VAAs) change from 0 to 360° in steps of 1°. Thus, the simulated dataset is composed of 16 cases and every case has 21,600 angle groups.

2.2. Airborne Dataset

Watershed Allied Telemetry Experimental Research (WATER) is a synthetic field experiment aiming to improve the observability, understanding, and predictability of hydrological and related ecological processes at watershed scale [24]. The campaign was conducted in the spring to summer of 2008 on the Heihe River Basin in northwest China (Zhangye City, Gansu Province). The WiDAS system is one of the major airborne sensors used in the WATER campaign, and has four charge-coupled device (CCD) cameras in visible/near-infrared (VNIR) channels and two thermal cameras in the mid infrared (MIR) and TIR channels [21]. These cameras can sequentially observe the surface with a very short time interval (4 s) and a very high overlapping ratio (more than 85% in the MIR and TIR channels) during the flight. Therefore, several sequential images mean that the same ground point can be almost simultaneously observed at different angles. According to the design, there are a total of seven zenith angles (between forward 40° and backward 40°) for the same ground point in the MIR and TIR channels. A multi-angular dataset can be obtained from the collections of the same ground point in the sequential WiDAS images.

Images on 7 July 2008 were chosen to develop the multi-angle dataset owing to their high data quality and clear sky. The solar zenith angle (SZA) and solar azimuth angle (SAA) in the study area were about 22° and 230°, respectively. Because the flight height was approximately 1.5 km above the surface, the spatial resolutions of the MIR/TIR cameras and the VNIR images are approximately 7.9 and 1.25 m, respectively. Six homogeneous sites, including settlement, wheat, maize, orchard, sea buckthorn, and bare soil [25], were selected for further evaluation. The number of images of each site indicates the number of simultaneous observations for the same pixel. Owing to changes in the VZA (and/or VAA), there is a temperature variation (T_diff) for every pixel. The mean T_diff varies from 1.4 K (maize) to 7.7 K (settlement). More detailed information can be found in Figure 1 and

Table 1. Specifications of the six sites.

No.	Long. (°E)	Lat (°N)	Land Cover	LAI	Area (m × m)	UTC (hh:mm:ss)	Number of Images	Mean Tdiff (K)
1	100.264	38.559	Settlement	\	375 × 280	04:13:14–04:13:42	8	7.7
2	100.245	38.515	Wheat	3.4	340 × 315	03:58:10–03:58:38	8	1.7
3	100.247	38.512	Maize	4.2	310 × 260	03:58:10–03:58:34	7	1.4
4	100.239	38.509	Orchard	2.4	330 × 320	03:58:26–03:58:54	8	3.9
5	100.232	38.494	Sea buckthorn	1.9	315 × 245	03:59:14–03:59:42	8	2.9
6	100.193	38.459	Bare soil	\	530 × 375	04:26:58–04:27:18	6	1.8

Note that the symbol\is no data, LAI is leaf area index.

.

2.3. Satellite Dataset

We selected publicly available MODIS and AATSR products as the data sources to develop a multi-angle LST dataset. MODIS LST products (MOD/MYD11_L2, collection-6) were retrieved with the generalized split-window algorithm applied to two TIR bands (31 and 32) [26]. The Level 2 (L2) product is a geophysical product with 1 km spatial resolution at nadir for a swath. Every day, the MODIS LST has a maximum of four clear sky observations. Emissivity values in Bands 31 and 32 were estimated by using the classification-based emissivity method [27] according to ‘static’ land cover types in the pixel. The generalized split-window algorithm is a viewing-angle-dependent algorithm in which six VZA subranges are divided to correct atmospheric and emissivity effects owing to the big scan angle (±65°) [26]. Radiance-based (R-based) validation indicated that the Collection 5 MODIS LST errors are within ±2 K (within ±1 K in most cases) for all the sites except for bare soil sites [28]. By adding two sets of coefficients separately for daytime and nighttime and adjusting the emissivity difference in Bands 31 and 32, the accuracy of bare soil pixels in the Collection 6 was improved [29]. AATSR LST products (ATS_LST_2P) were generated with another split-window algorithm from brightness temperatures measured at 11 and 12 μm for the nadir view of AATSR [30]. Coefficients of the algorithm depend on the land cover type, the fractional vegetation cover (FVC), the atmospheric water vapor, and the viewing angle. Because of the narrow swath width, the temporal resolution of AATSR LST is about three days and the spatial resolution is 1 km. Unlike the case of MODIS LST, the land surface emissivity is implicitly taken into account through the land cover type and FVC-dependent coefficients. According to the results of comparisons between AATSR and in situ LST data over eighteen sites [31], the AATSR LST approaches consistency within the target accuracy (2.5 K) during the day. During the night, the accuracy cannot meet the target (1 K) but the overall error is lower than for the day.

Relevant auxiliary products include global geolocation products (MOD03/MYD03 and ATS_NR_2P) which provide viewing and illumination geometries and a land cover type product (MCD12Q1) which provides yearly land cover classification with a 500 m resolution. In order to be consistent with the canopy scale, only the daytime observations with solar zenith angle (SZA) less than 90° were used in this paper. The establishment of the dataset involved five main steps: (1) LST quality control; (2) spatial collocation; (3) temporal concurrence; (4) spatial homogeneity; and (5) systematic difference removal.

• LST Quality Control

To ensure the availability of the dataset, only pixels with the highest quality were used in this study. For MOD/MYD11_L2 products, the LST value is valid and the quality control (QC) flag is 0 [32]. For ATS_LST_2P products, the LST_uncertainty flag is less than 1 K and the QC flag is 16 [31].

• Spatial Collocation

We identify the same observation location of the two satellites by longitude and latitude. In this paper, the threshold was set at 0.001°. Given the pixel spatial resolution (1 km), the distance is a reasonable value.

• Temporal Concurrence

This step is to check whether the spatially collocated pixels with the best quality are concurrent in time. The interval of the MOD/MYD11_L2 product is five minutes, and that of the ATS_LST_2P product is about 100 min. Therefore, we firstly divided the ATS_LST_2P products into the same interval as the MOD/MYD11_L2 products. Then, the three nearest ATS_LST_2P products were chosen to match for each MOD/MYD11_L2 product. In this way, the biggest possible time difference is 10 min.

• Spatial Homogeneity

To analyze the influence of land cover type on parametric models and further minimize location matching errors, only pixels over relatively homogenous areas were selected. There are two filter criteria, as follows.

(1)

The pixel’s LST should not exceed a certain variability in a 5 × 5 moving window and the standard deviation (STD) should be less than 1 K [33];

(2)

If more than 23 pixels have the same land cover type as the centered pixel inside the 5 × 5 window, the centered pixel is considered to be in a ‘‘pure homogenous or quasi-homogenous area’’ and should be retained [34].

• Systematic Difference Removal

Because of the intrinsic differences in LST retrieval algorithms and satellite sensors, there are still some differences between the two LST products even with the same VZA. The difference is defined as the system error and can be adjusted by a simple linear regression [20]. This step is similar to an inter-comparison between MODIS and AATSR LST products. Firstly, matched pixel pairs with aligned viewing geometries were selected by restricting the maximum relative difference in the secant of the VZA to 0.01 [35]. Based on selected pixel pairs, we built linear regression models between the two LST products for each land cover type. Finally, we chose MODIS LST products as the true value and normalized the AATSR LST.

After the above processing, we firstly established the LST multi-angle dataset at LEO satellite scale. The dataset had 647,788 matched pixel pairs and included 14 International Geosphere- Biosphere Programme (IGBP) land cover types (not including closed shrublands, urban and built-up areas, and snow and ice). Just as Figure 2 shows, the ranges of SZA, VZA, and the relative azimuth angle (RAA) are 20–89°, 0–65°, and 0–360°, respectively. Compared with previous studies [4,20,22], the dataset has higher spatial resolution, wider land cover types, and richer angular sampling. Therefore, it should be able to provide more comprehensive analysis.

2.4. SURFRAD Dataset

The SURFRAD was established in 1993 to support climate research with accurate, continuous, long-term measurements of the surface radiation budget over the United States. A full year (2011) of ground truth LSTs simultaneous with MODIS LSTs were calculated from observed upward ($F_u$) and downward ($F_d$) wideband hemispheric infrared fluxes at the six SURFRAD stations listed in

Table 2. Specifications of six Surface Radiation Budget Network (SURFRAD) stations.

Station Name	Land Cover	Long. (°W)	Lat (°N)	Number of Observations
Bondville, IL	Cropland	88.373	40.051	52
Boulder, CO	Bare soil	105.238	40.126	35
Fort Peck, MT	Grassland	105.102	48.308	47
Goodwin Creek, MS	Grassland	89.873	34.255	171
Penn State, PA	Cropland	77.931	40.720	23
Sioux Falls, SD	Grassland	96.623	43.734	68

. LST is computed using the traditional equation:

$$LST={\left\{\frac{F_u-\left(1-\epsilon \right)\times F_d}{\sigma ·{\epsilon }_{3-\infty }}\right\}}^{0.25}$$

(1)

where $\sigma$ is the Stefan–Boltzmann constant (5.670373 × 10⁻⁸) and ${\epsilon }_{3-\infty }$ is the broadband surface emissivity. The value ${\epsilon }_{3-\infty }$ can be estimated from MODIS narrowband emissivities via the following method [36]:

$${\epsilon }_{3-\infty }=0.0127+0.7852\times {\epsilon }_{29}-0.0151\times {\epsilon }_{31}+0.2139\times {\epsilon }_{32}$$

(2)

where ${\epsilon }_{29}$, ${\epsilon }_{31}$, and ${\epsilon }_{32}$ are the narrowband emissivities at the MODIS 29, 31, and 32 bands, respectively. In this study, MODIS emissivity products (MOD/MYD11B1) [32] were used to obtain the narrowband emissivity. Unlike the MOD/MYD11A1 products, MOD/MYD11B1 products use a physics-based day/night algorithm to dynamically retrieve emissivities at seven spectral bands and LSTs at a coarser spatial resolution (5.6 km). Combining the narrowband emissivities, $F_u$, $F_d$, and Equations (1) and (2), a SURFRAD LST dataset can be established.

3. Parametric Models

3.1. RL Model

Lagouarde and Irvine [8] validated the suitability of a parametric hotspot model proposed by Roujean [37] for the optical domain to the TIR domain, referred to hereafter as the RL model. The RL model describes thermal radiation directional anisotropy requiring two parameters, and is expressed by the following equation:

$$\Delta T\left({\theta }_s,{\theta }_v,\phi \right)=\Delta T_{HS}\frac{\left[e^{-kf}-e^{-k{f}_N}\right]}{\left[e^{-k{f}_{HS}}-e^{-k{f}_N}\right]},$$

(3)

$$\Delta T\left({\theta }_s,{\theta }_v,\phi \right)=T\left({\theta }_s,{\theta }_v,\phi \right)-T_N,$$

(4)

$$\Delta T_{HS}=T_{HS}-T_N,$$

(5)

$$f=\sqrt{{\tan }^2{\theta }_s+{\tan }^2{\theta }_v-2\tan {\theta }_s\tan {\theta }_v\cos \phi },$$

(6)

where ${\theta }_s$, ${\theta }_v$, and $\phi$ are the SZA, VZA, and RAA, respectively. Subscript N stands for nadir observation and $f_N=\tan {\theta }_s$. Subscript HS is a hotspot when the viewing direction coincides with the sun direction and $f_{HS}=0$ [17]. Thus, Equation (3) can be rewritten as

$$\Delta T\left({\theta }_s,{\theta }_v,\phi \right)=\Delta T_{HS}\frac{\left[e^{-kf}-e^{-k\tan {\theta }_s}\right]}{\left[1-e^{-k\tan {\theta }_s}\right]}.$$

(7)

In the RL model, the function $\Delta T\left({\theta }_s,{\theta }_v,\phi \right)$ is used to estimate the anisotropy. The two parameters are $\Delta T_{HS}$, which governs the anisotropy at the hotspot, and $k$, which adjusts the shape of the variations of anisotropy along with the VZAs.

3.2. BRDF Model

The BRDF model was firstly proposed to simulate bidirectional reflectivity under different viewing angles. By replacing reflectivity with surface temperature in the RL model, the BRDF model was updated to the TIR domain [18,19]. In this paper, we use the best-fitting kernel grouping of the RossThick volume kernel and the LiSparseR geometry kernel; the formulation is as follows:

$$T({\theta }_v,{\theta }_s,\phi )=f_{iso}+f_{vol}·k_{vol}({\theta }_v,{\theta }_s,\phi )+f_{geo}·k_{geo}({\theta }_v,{\theta }_s,\phi )$$

(8)

$$k_{vol}({\theta }_v,{\theta }_s,\phi )=\frac{(\pi /2-\xi )\cos \xi +\sin \xi }{\cos {\theta }_v+\cos {\theta }_s}-\frac{\pi }{4}$$

(9)

$$\cos \xi =\cos {\theta }_v\cos {\theta }_s+\sin {\theta }_v\sin {\theta }_s\cos \phi$$

(10)

$$k_{geo}({\theta }_v,{\theta }_s,\phi )=O({\theta }_v,{\theta }_s,\phi )-\sec {\theta }_v^{\prime }-\sec {\theta }_s^{\prime }+\frac{1}{2}(1+\cos {\xi }^{\prime })\sec {\theta }_v^{\prime }\sec {\theta }_s^{\prime }$$

(11)

$$O=\frac{1}{\pi }\left(t-\sin t\cos t\right)\left(\sec {\theta }_v^{\prime }+\sec {\theta }_s^{\prime }\right)$$

(12)

$$\mathrm{cost}=\frac{h}{b}\frac{\sqrt{D^2+{(\tan {\theta }_v^{\prime }\tan {\theta }_s^{\prime }\sin \phi )}^2}}{\sec {\theta }_v^{\prime }+\sec {\theta }_s^{\prime }}$$

(13)

$$D=\sqrt{{\tan }^2{\theta }_v^{\prime }+{\tan }^2{\theta }_s^{\prime }-2\tan {\theta }_v^{\prime }\tan {\theta }_s^{\prime }\cos \phi }$$

(14)

$$\cos {\xi }^{\prime }=\cos {\theta }_v^{\prime }\cos {\theta }_s^{\prime }+\sin {\theta }_v^{\prime }\sin {\theta }_s^{\prime }\cos \phi$$

(15)

$${\theta }_v^{\prime }={\tan }^{-1}(\frac{b}{r}\tan {\theta }_v)$$

(16)

$${\theta }_s^{\prime }={\tan }^{-1}(\frac{b}{r}\tan {\theta }_s)$$

(17)

where $f_{iso}$, $f_{vol}$, $f_{geo}$, $k_{vol}({\theta }_v,{\theta }_s,\phi )$, $k_{geo}({\theta }_v,{\theta }_s,\phi )$ are the isotropic scattering kernel coefficient, the volumetric kernel coefficient, the geometric kernel coefficient, the volumetric kernel, and the geometric kernel, respectively; O is the overlap area between the view and solar shadows; superscript ′ stands for effective angle adjusted for spheroidal shape of the crown; and the ratios $\frac{h}{b}$ and $\frac{b}{r}$ are the dimensionless crown relative height and shape parameters, preselected as 2 and 1, respectively.

3.3. Vinnikov Model

Vinnikov et al. [4] developed an empirical model (hereafter called the Vinnikov model) to normalize surface temperature to the nadir direction by means of a statistical approach. It can be expressed by the following simple equations:

$$\frac{T\left({\theta }_v,{\theta }_s,\phi \right)}{T_0}=1+A·\phi ({\theta }_v)+D·\psi ({\theta }_v,{\theta }_s,\phi )$$

(18)

$$\phi ({\theta }_v)=1-\cos ({\theta }_v)$$

(19)

$$\psi ({\theta }_v,{\theta }_s,\phi )=\sin ({\theta }_v)\cos ({\theta }_s)\sin ({\theta }_s)\cos ({\theta }_s-{\theta }_v)\cos (\phi )$$

(20)

where $T_0$ is the temperature at nadir; and 1, $\phi ({\theta }_v)$, $\psi ({\theta }_v,{\theta }_s,\phi )$, A, and D are the isotropic kernel, emissivity kernel, solar emissivity, emissivity kernel coefficient and solar emissivity coefficient, respectively.

4. Evaluations and Results

4.1. Evaluation Using Simulations

Just as Section 3 shows, the three models do not use the same definition of thermal radiation directional behavior. For convenient interpretation of results, we selected the definition of the RL model which considers anisotropy as the temperature difference between nadir and off-nadir observations. Therefore, the BRDF model becomes

$$\Delta T({\theta }_v,{\theta }_s,\phi )=f_{vol}·[k_{vol}({\theta }_v,{\theta }_s,\phi )-k_{vol}(0,{\theta }_s,0)]+f_{geo}·[k_{geo}({\theta }_v,{\theta }_s,\phi )-k_{geo}(0,{\theta }_s,0)],$$

(21)

and the Vinnokov model becomes

$$\Delta T({\theta }_v,{\theta }_s,\phi )=T_0·\left[A·\phi ({\theta }_v)+D·\psi ({\theta }_v,{\theta }_s,\phi )\right].$$

(22)

For every canopy structure, ordinary least squares linear regression without the intercept was used to fit the BRDF model and the Vinnikov model, and a nonlinear optimization procedure based on the Nelder–Mead method [38] was used to solve the RL model. Figure 3 presents the performance of the three parametric models in describing anisotropy based on the 4SAIL simulated dataset. Errors for most records are in the range of [−0.1 K, 0.2 K] which indicates that the parametric models can estimate directional behaviors well. The Vinnikov model has the best ability, with 94.25% errors within ±0.1 K. Next is the BRDF model, for which 80.80% errors are within ±0.1 K. Then is the RL model, for which the proportion is 75.60%. Moreover, the peak of the histogram is on the interval from 0 to 0.1 K and the trend slant to the right illustrates that the three models slightly overestimate surface temperature. The proportions of positive errors of the RL, BRDF, and Vinnikov models are 73.48%, 62.89%, and 75.11%, respectively.

Figure 4a–d show the anisotropy obtained using parametric models adjusted using 4SAIL under different LIDFs. Obviously, the Vinnikov model can achieve a good accuracy for all LIDFs and presents high robustness. Root mean square errors (RMSEs) for spherical, erectophile, planophile, and plagiophile canopy are 0.09 K, 0.12 K, 0.11 K, and 0.10 K, respectively. For erectophile canopy, the RL and BRDF models have terrible performance, with RMSEs of 0.67 K and 0.38 K, respectively. For the other three LIDFs, there is a clear underestimation of anisotropy for large ΔT. In other words, all parametric models cannot describe anisotropy on hotspots well. Relatively speaking, the BRDF model has the best results.

In order to investigate the influence of LAI on parametric model fitting ability, we added LAIs from 0.5 to 5 with a step size of 0.5. Relative RMSE (RRMSE) was also calculated to eliminate the anisotropy differences between different LAIs. The result is shown in Figure 5. For erectophile and spherical canopy, RMSEs and RRMSEs of the anisotropy fitting error maintained the same trend. As the LAI increases, fitting errors firstly decrease and then increase. At a LAI of 0.5, the RMSE and RRMSE of the Vinnikov model are clearly less than those of the other two models, suggesting that the Vinnikov model has a prominent advantage at low LAI. For planophile and plagiophile canopy, the RMSEs and RRMSEs do not stay the same when the LAI is less than 2. However, the RMSEs are small while the RRMSEs are big, which indicates that the parametric models cannot give a good performance. Therefore, regardless of the LIDFs, fitting effects firstly improve and then decline with increase in the LAI, and the best performance mainly occurs when the LAI varies from 1 to 2. This phenomenon can be explained by the fact that thermal radiation directional signatures are weak at very low and high LAIs [14,16].

The RL model and the BRDF model extended from the VNIR domain can simulate shadowing effects well [37,39]. The Vinnikov model has the emissivity kernel, which is an important difference from the other two models. For low LAI and an erectophile LIDF canopy, the shadowing effect is weak but the emissivity directionality is significant [14]; this may lead to a clearly better performance from the Vinnikov model. On the other hand, the RL and BRDF models are based on the circular crown assumption [37,39]. This feature may also restrict their ability on a low LAI and erectophile LIDF canopy.

Hu et al. [40,41] suggested that two different kernel groups be used for different LAIs to improve the ability of the BRDF model. However, estimating the LAI accurately is not an easy task in practice and would also bring more parameters into the model. Moreover, the RL model cannot be easily inverted because of its exponential formula [41]. In some cases, a multidimensional unconstrained nonlinear optimization procedure can even generate an unrealistic $\Delta T$ which has a clear physical meaning. Based on robustness, simplicity, and the significant advantage for a low LAI and erectophile LIDF canopy, we recommend using the Vinnikov model to simulate directional thermal radiation.

It should be noted that all analyses above using simulations are based on the 4SAIL model. Evaluation results are therefore influenced by the accuracy of the 4SAIL model, which is a limitation of this study. Moreover, the 4SAIL model is suitable for turbid canopy. Evaluation for performance on discrete canopies, such as discontinuous savanna canopy using the modified geometric projection (MGP) model [42] and sparse vegetation using the DART model [13], are intended as further work.

4.2. Evaluation Using Airborne Data

As mentioned in Section 2.2, the least number of simultaneous observations for the same pixel is 6. Therefore, we can fit the three models for every pixel and then evaluate their performances. For two different observations T₁ and T₂, the RL model can be expressed as

$$T_1-T_2=\Delta T_{HS}\left(\frac{\left[e^{-k{f}_1}-e^{-k\tan {\theta }_{s_1}}\right]}{\left[1-e^{-k\tan {\theta }_{{}_{s_1}}}\right]}-\frac{\left[e^{-k{f}_2}-e^{-k\tan {\theta }_{s_2}}\right]}{\left[1-e^{-k\tan {\theta }_{{}_{s_2}}}\right]}\right)$$

(23)

the BRDF model turns out to be

$$T_1-T_2=f_{vol}·\left[k_{vol}({\theta }_{v_1},{\theta }_{s_1},{\phi }_1)-k_{vol}({\theta }_{v_2},{\theta }_{s_2},{\phi }_2)\right]+f_{geo}·[k_{geo}({\theta }_{v_1},{\theta }_{s_1},{\phi }_1)-k_{geo}({\theta }_{v_2},{\theta }_{s_2},{\phi }_2)]$$

(24)

and the Vinnikov model becomes

$$T_1-T_2=A·\left[T_2·\phi ({\theta }_{v_1})-T_1·\phi ({\theta }_{v_2})\right]+D·\left[T_2·\psi ({\theta }_{v_1},{\theta }_{s_1},{\phi }_1)-T_1·\psi ({\theta }_{v_2},{\theta }_{s_2},{\phi }_2)\right]$$

(25)

where subscripts 1 and 2 denote different observations. Similar to in the simulation evaluation, ordinary least squares linear regression without the intercept was used to solve the BRDF model and the Vinnikov model, and a nonlinear optimization procedure was used to calculate the RL model. As Figure 6 shows, the RRMSEs of the six sites are approximately 2–6% and the RMSEs are less than 1 K (except for Site 1), which also illustrates that the three models can simulate directional anisotropy well. Similar to the findings from the simulations, there is a clear tendency of the Vinnikov model to perform better than the BRDF model, and the BRDF model is better than the RL model. Synthesizing the RMSE and RRMSE, we find that the fitting accuracy decreases from herb plants (wheat and maize), to woody plants (sea buckthorn and orchard), to bare soil, to settlement. The performance on bare soil of the parametric models is poor, which is consistent with the low LAI case from the simulations. It is worth noting that wheat and maize have a higher LAI (see Table 1) but a better fitting accuracy than sea buckthorn and orchard, which is inconsistent with the result from simulations. Sea buckthorn and orchard have more complex canopy structures, such as plant height, plant spacing, and even three-dimensional structure, which is more significant for settlement. The canopy parameters can also exert an influence on thermal radiation anisotropy [9,11,42] and therefore restrict the capability of the three parametric models.

4.3. Evaluation Using Satellite Data

Along with at the canopy scale, parametric models were calibrated for each land cover type at the pixel scale. Then, all LSTs were corrected to nadir observation using calibrated models. We selected the root mean square differences (RMSD) between MODIS and AATSR LST as the assessment index for model evaluation. A low RMSD means high comparability between the two products and good performance of the parametric model. The evaluations were conducted in two different ways: verification and validation. The LST multi-angle dataset was randomly divided into two parts in which one was 2/3 of the total and the other was 1/3. We used the 2/3 dataset to verify the models and obtain coefficients. Then, the coefficients are applied as reference and the remaining 1/3 dataset was used to validate the parametric models.

The result is shown in Figure 7. Firstly, removing system error can dramatically decrease the RMSD between two LST products and angular correction can further improve their comparability. Whether verification or validation, models have the same results, which indicates that the three parametric models are stable. For most land cover types, the three models have similar effects. However, in the deciduous broadleaf forest (land cover Class 4), savannas (land cover Classes 8 and 9), grassland (land cover Class 10), cropland (land cover Class 14) and barren or sparsely vegetated (land cover Class 16), the Vinnikov model is better than the other two models. This result is consistent with results from the canopy scale that the Vinnikov model has the best ability in dealing with erectophile canopy (savannas and grassland) and low LAI (barren or sparsely vegetated).

5. Application

From the above analysis, we conclude that the Vinnikov model can obtain more accurate and robust results at both canopy and pixel scales. For this reason, we applied this model to LST angular correction. The LST is extremely dependent on the viewing and illumination geometries, topography, and land cover type [22]. Because of the spatial homogeneity step during the development, the LST multi-angle dataset at LEO satellite scale mainly consists of relatively flat and homogeneous pixels [35]. Therefore, topography is not considered in this study. For each land cover type, the coefficients A and D in Equation (18) can be obtained through an ordinary least squares linear regression method using the LST multi-angle dataset.

From

Table 3. Specification of the Vinnikov model for different land cover classes.

Land Cover Class	A	D	Number of Points	RMSE (K)	Land Cover Class	A	D	Number of Points	RMSE (K)
0	−0.0067	−0.0015	18,073	0.85	8	−0.0178	0.0044	33,697	1.14
1	−0.0068	−0.0002	3284	0.72	9	−0.0175	0.0045	59,706	1.35
2	−0.0173	0.0046	47,542	0.91	10	−0.0228	0.0005	78,303	1.62
3	−0.0102	0.0034	15,485	0.79	11	−0.0115	0.0024	22,325	0.87
4	−0.0214	−0.0023	5602	0.98	12	−0.0184	0.0022	40,985	1.44
5	−0.0093	0.0016	34,211	0.86	14	−0.0279	0.0041	32,568	1.57
7	−0.0045	−0.0078	5612	0.84	16	−0.0209	−0.0005	250,395	1.42

, one can see that the values of A are negative. There is a decrease in emissivity as the VZA increases, and the feature of emissivity will lead to a higher LST value retrieved for fixed emissivity [34]. Therefore, A is negative so that the emissivity kernel can reduce this effect. For land cover Classes 4, 8, 9, 10, 14, and 16 where the Vinnikov model has better performance, there are relatively bigger absolute values of A, suggesting a more significant emissivity effect. Using Spinning Enhanced Visible and Infrared Imager (SEVIRI) and MODIS LST products over four years, Ermida et al. [20] pointed out that positive values of D contribute to a decrease of the LST when VZAs deviate from nadir, except for the hotspot effect. However, coefficient D is negative for land cover Classes 0, 1, 4, 7, and 16. For bare soil (Class 16), there may be no hotspot effects owing to weak shadowing effects. For other land cover classes (Classes 0, 1, 4, and 7), hotspot effects may not have been observed due to insufficient resamples. The two cases can lead to the maximum LST appearing in nadir and a negative D. Moreover, evergreen broadleaf forest (Class 2), woody savannas (Class 8), savannas (Class 9) and cropland/natural vegetation mosaic (Class 14) have a relatively bigger D (≥0.004), and grassland (Class 10) and bare soil (Class 16) have a relatively smaller D (±0.0005). This result is consistent with that of Ermida et al. [20] that there is a more significant shadowing effect on the heterogeneous area characterized by tree (or tall shrubs) coverage than on the homogeneous areas (bare ground or highly dense forests). Ermida et al. also analyzed the relationship between topography and percentage of tree cover (PTC) and D pixel-by-pixel benefit from enormous matched points. We do not have the relative PTC data and cannot conduct a more detailed analysis; this is a limitation in our study.

Firstly, we applied the Vinnikov model and coefficients in Table 3 to the MODIS LST at SURFRAD stations. Before the angular correction to nadir observation, the MODIS LST was lower than the SURFRAD-measured LST, and the RMSE and bias were 2.47 K and −2.05 K, respectively (see Figure 8a). After angular correction, the MODIS LST was improved effectively and the RMSE and bias decreased to 1.58 K and −0.39 K (see Figure 8b).

Then, the MODIS LST at 08:00 on 25 June 2011 (UTC) was used to exemplify the performance of the Vinnikov model. Figure 9a displays the MODIS LST before angular correction and Figure 9b is the result of LST correction to nadir. A comparison of Figure 9b with Figure 9a shows that the spatial variations in the LST after angular correction appeared to be more uniform than those in the LST before angular correction. Figure 9c presents differences between the LST before and after angular correction, which range from 0 to 5 K. Because the effects of SZA and RAA are very small in a swath of MODIS, only VZAs are shown (Figure 9d). A comparison of Figure 9c with Figure 9d indicates that the LST differences are highly related to the VZA.

6. Discussion

By using a SCOPE-simulated anisotropy dataset with six LAIs (0.5, 1, 1.5, 2, 3, 5) and four hotspot parameters (0.01, 0.05, 0.1, 0.5) for spherical LIDF, Duffour et al. [17] compared the RL model and the Vinnikov model. They revealed that the RL model is more efficient than the Vinnikov model at simulating TIR anisotropy, particularly close to the hotspot. Except for the differences between 4SAIL and SCOPE, there are two main reasons for this slight inconsistency. Firstly, the foci of the two studies are different. Duffour et al. [17] paid attention to different hotspot effects while we account for different LIDFs. For the same canopy structure, as Figure 5a shows, the RL model is better than the Vinnikov model when the LAI is greater than 2; this is consistent with the results of Duffour et al. [17]. Secondly, as a purely empirical model, the Vinnikov model does not consider the ‘distance’ between the sun and view directions (Equation (2) for the RL model and Equation (14) for the BRDF model). Failure to consider such a term may severely limit the Vinnikov model’s ability to simulate the hotspot anisotropy [17]. However, the hotspot effect in the TIR domain is not conclusive. Duffour et al. [16] indicate that (1) literature about the hotspot effect essentially deals with the VNIR domain, and values reported display a large range of variation; (2) experimental data are too scarce to be fully confident about the hotspot effect that may be assigned in TIR; (3) the hotspot effect can be set from several definitions. For these reasons, we do not analyze the sensitivity of the parametric models to the hotspot effect by prescribing different hotspot parameters as did Duffour et al. [17].

The quality of reference datasets plays a vital role in the evaluation of parametric models. The 4SAIL model and the WiDAS dataset have been successfully applied to many relevant studies, i.e., directional thermal radiative simulation [14,43], LST angular normalization [44], and component temperature separation [21,45]. For the pixel scale, the LEO satellite LST multi-angle dataset in the study is the first of its kind. On the one hand, similar to the previous GEO satellite LST multi-angle dataset [4,20,22], the LEO satellite LST multi-angle dataset has the following uncertainties: (1) The emissivity uncertainty. The MODIS LST uses a fixed emissivity and the AATSR LST uses the land cover type and FVC to implicitly consider the emissivity. The land emissivity has directional behavior which is also influenced by the FVC and cavity effects [34,46]. Although the retrieval algorithms of MODIS LST and AATSR LST are viewing angle dependent, the emissivity anisotropy has not been fully considered. Moreover, emissivity differences between the two LST products can also affect the coefficients of the parametric models [20]. However, the difference cannot be analyzed quantitatively using the LEO satellite LST dataset because the AATSR LST lacks explicit emissivities. (2) The LST uncertainty. The uncertainty has various sources, i.e., the difference in observation areas at different viewing angles, the LST variation between two different observation times, and atmospheric correction errors. Generally, these uncertainties would become higher and have more significant influence with higher LSTs. On the other hand, unlike in previous studies, the model evaluation is based on the land cover type rather than the pixel. Differences within the same land cover type can also increase the uncertainty of the dataset.

In this study, five steps were set for establishment of the dataset. Among them, spatial collocation and temporal concurrence are used to reduce temporal–spatial differences between the two LST products; spatial homogeneity is used to reduce differences within the same land cover type; and LST quality control and systematic difference removal are used to reduce atmospheric and emissivity effects. Comparison results at the pixel scale are agreeable with the canopy scale. Parameters A and D show features identical to those in the findings of Ermida et al. [20]. All the results suggest that the LEO satellite LST multi-angle dataset is valid for analyzing the parametric models.

Nevertheless, there are still two intrinsic defects owing to data characteristics. The first one is the angular asymmetry between the MODIS LST and the AATSR LST. Unlike MODIS LSTs for which VZAs are from 0 to 65°, the range of VZAs of AATSR LSTs is only from 0 to 22.5°. Therefore, all linear regression models for system error correction are, in theory, only appropriate for the small angle range of AATSR LSTs. This angular asymmetry may affect the results of the above evaluation. However, it is impossible to assess the impact subject to AATSR characteristics. Secondly, MCD12Q1 is a yearly land cover type product. There are different thermal radiation directional signatures during different seasons and growing periods even for the same land cover type [22]. Unfortunately, the seasonal effect cannot be assessed using the MCD12Q1 product, which is another limitation of this study. One solution for solving these problems is to apply new satellite products—for example, using LSTs with large scanning angles and daily surface type products from visible infrared imager radiometer sensors (VIIRSs) or other possible satellites in joint polar satellite systems (JPSSs).

7. Conclusions

This paper evaluated the performance of three parametric directional thermal radiation models—the RL model, the BRDF model, and the Vinnikov model—at canopy and satellite scale using different data sources. Specially, at canopy scale, two different datasets were used: one was an extensive simulated dataset which consisted of 16 different canopy structures by prescribing 4 LIDFs and 4 LAIs; the other was an airborne multi-angle dataset including settlement, wheat, maize, orchard, sea buckthorn, and bare soil derived from WiDAS images. The evaluation showed that the Vinnikov model had the best accuracy and robust ability, particularly for erectophile canopy or low LAI. The accuracies of the three models firstly improved and then declined with increasing LAI, and the best performance occurred for LAI between 1 and 2. The three models could not describe the hotspot anisotropy very well and were also restricted by other parameters, such as plant height, plant spacing, and three-dimensional structure. At pixel scale, the first LEO satellite LST multi-angle dataset was established by combining MODIS and AATSR products. Whether in verification or validation, the three models presented a stable effect. However, the Vinnikov model had the best result in the erectophile canopy (savannas and grassland) and low LAI (barren or sparsely vegetated) areas, which was consistent with the results at the canopy scale. For this reason, we calibrated the Vinnikov model according to different land cover types to correct the MODIS LST product. Validation using SURFRAD-measured LSTs indicated that angular correction could reduce the RMSE of 0.89 K and bias of 1.66 K, respectively. In addition, the corrected result showed better spatial uniformity and higher angular correlation.

The three parametric models are based more or less on the knowledge of directional behavior of the VNIR domain. Indeed, sunlit/shaded elements have different radiations, resulting in different temperatures in the TIR domain. However, thermal radiation has inertia effects and the temperature differences have a certain lag with the moving sun, which is distinct from the VNIR domain. Though the three parametric models each present a good performance in this study, inertia effects still require further research and interpretation. In addition, validations of the three parametric models for discontinuous savanna canopy and sparse vegetation using unmanned aerial vehicle multi-angle observations or using new LEO satellites will be conducted in future work.