On the Derivation of Geometric Optical Kernels for Directional Thermal Radiation

The derivation of widely used geometric optical (GO) kernels in bidirectional reflectance distribution function models, that is, LiSparseReciprocal kernel (KGOLSR) and LiDenseReciprocal kernel, was based on two important assumptions: (1) The shaded components are perfect black and (2) the contributions of two sunlit components are identical. Different from the bidirectional reflectance, thermal radiation directionality effect mainly results from component temperature differences, suggesting the above assumptions are not applicable in most situations. Therefore, this study derived GO kernels for thermal radiation based on temperature differences rather than illumination differences. Specifically, four GO kernels, that is, KGO4 with considering sunlit/shaded vegetation and sunlit/shaded soil, KGO3 with considering sunlit/shaded soil and vegetation, KGO2 with considering vegetation and soil, and KGOg only considering the hottest sunlit soil, have been developed. By using a comprehensive simulated data set, their performances have been thoroughly evaluated and the comparison with KGOLSR has also been analyzed in depth. Results showed that (1) KGO4 had the highest accuracy and KGO3 was the second; in the case of only two available angles, KGOg performed best. (2) Variables such as component temperature, component emissivity, solar zenith angle, and the percentage of tree crown cover mainly affected the comparison result between KGOLSR and KGO2; KGOLSR would have a better performance for a scene with a stronger vegetation effect. Moreover, the best values of two structure characteristics, that is, the crown shape parameter b/r and relative height h/b, for these five kernels have been determined, which can provide instruction for practical application.


Introduction
Land surface temperature (LST) is a key variable for climate change and energy balance, which has been widely used to a variety of geoscientific studies (Li et al., 2013;Tang, 2018;Tang et al., 2008). To estimate this important parameter, many retrieval algorithms have been proposed and many operational products have also been generated. However, all of them generally have an assumption of isotropic thermal emission (Li et al., 2013). In fact, like the bidirectional reflectance in visible and near-infrared (VNIR) region, thermal infrared (TIR) radiation has been proven to own a significant directional signature (Lagouarde et al., 1995;Lagouarde et al., 2004;Lagouarde et al., 2010;Zhang et al., 2000). This thermal radiation directionality (TRD) effect is mainly resulted from two reasons: (1) the land surface emissivity is angularly dependent Cuenca & Sobrino, 2004;García-Santos et al., 2012;García-Santos et al., 2015;Hu et al., 2019;Sobrino et al., 2005) and (2) various components of heterogeneous surface have different temperatures Bian et al., 2018;Li et al., 2019;Timmermans et al., 2009). As a result, there are different observation values of thermal radiation or LST under different viewing and illumination geometries (Ermida et al., 2014;Ermida et al., 2017). Obviously, the TRD can reduce the intercomparison of LST and therefore limits its wide application . In order to address this problem, lots of thermal radiation © 2020 The Authors. Earth and Space Science published by Wiley Periodicals, Inc. on behalf of American Geophysical Union This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. directional models, including radiative transfer models, geometric optical (GO) models, hybrid models, three-dimensional models, and parametric models, have been developed, and a detailed review can be found in Cao, Liu, et al. (2019).
The thermal radiation of a scene is derived from two parts: One is the thermal emission of individual components and the other one is the scattering radiation within components. The contribution of the component thermal emission is the predominate part , which can be properly described by using GO models. Generally, the GO models assume the scene can be divided into four different components, that is, sunlit/shaded soil and sunlit/shaded vegetation, and the total radiation is the sum of the radiations emitted by these four components weighted by their respective projected fractions. Therefore, the GO model has clear physical meaning and analytical formulations of which the input parameters are relatively available. For this reason, Pinheiro et al. (2004) upgraded the Li-Strahler Geometric-Optical Mutual Shadowing (GOMS) model (Li & Strahler, 1992) in VNIR region to TIR domain, that is, modified geometric projection (MGP) model, by replacing the bidirectional reflectance with thermal radiation. Then, MGP model has been tested against Discrete Anisotropic Radiative Transfer, a rigorous three-dimensional radiative transfer model, and in situ measurements presenting an accurate result (Pinheiro et al., 2006). Based on MGP model, Rasmussen et al. (2010Rasmussen et al. ( , 2011 analyzed the angular effect of LST retrieved from Geostationary Satellites with simulated and field-measured component temperatures. Ermida et al. (2014) proposed to use a numeric method rather than the original analytical method to simplify the estimation of the overlap area, which significantly improves the practicality of MGP model. Meanwhile, owing to its good simulation ability, MGP model has also been used as a reference for evaluating the performance of parametric models (Ermida et al., 2018).
In VNIR domain, semiempirical kernel-driven bidirectional reflectance distribution function (BRDF) models are the basic tool for the correction of bidirectional reflectance (Jupp, 2000). These models are generally a linear combination of isotropic kernel, GO kernel, and volume scattering kernels. Different from bidirectional reflectance, TRD is mainly depended on component temperatures rather than illumination. However, in the consideration of the fact that the temperatures of sunlit components are generally higher than the temperatures of shaded components in a scene , some BRDF models, for example, the RossThick-LiSparseR model (Peng et al., 2011;Ren et al., 2015) and RossThick-LiDenseR model (Hu et al., 2016), have been directly extended to simulate the TRD by replacing the bidirectional reflectance with thermal radiation. Moreover, following the similar methodology of kerneldriven model, Su et al. (2002) proposed LSF-Li model for simulating directional emissivity which is a group of isotropic kernel, GO kernel and two-layer kernel. The GO kernels for these models are LiSparseReciprocal kernel, hereafter named K GOLSR , or LiDenseReciprocal kernel. Actually, the derivation of these two widely used GO kernels was based on two important approximations. The first is shaded components can be treated as ideal blackbody as their contributions are small. The second is the contributions of sunlit components are identical (Wanner et al., 1995). For thermal radiation, the radiances of shaded components cannot equal to zero, and the shaded soil is even hotter than the sunlit vegetation (Ermida et al., 2014;Li et al., 2019;Rasmussen et al., 2011). On the other hand, the radiation contributions of sunlit vegetation and sunlit soil are not identical in the majority of cases (Bian et al., 2018;Bian et al., 2016;Ermida et al., 2014;Rasmussen et al., 2011;Timmermans et al., 2009). Therefore, these two assumptions for reflectance are not applicable to thermal radiation. If one uses BRDF models in VNIR domain to simulate the TRD by just using thermal radiation to replace the reflectance, it may bring some significant errors.
From the inherent characteristic of thermal radiation, this study simplifies the number of components based on the component temperature differences. Correspondingly, four GO kernels for thermal radiation have been developed. Their performances as well as the comparison with K GOLSR have been analyzed by using an extensive simulation data set, which is generated by using MGP model and prescribing 36 component temperatures, 6 component emissivities, 3 solar zenith angles (SZAs), and 6 percentage of tree crown covers (PTCs). Moreover, the best values for their structure parameters have also been determined. The outline of the paper is as follows. Section 2 is a description for the simulated data set. Section 3 illustrates the derivations of GO kernels. The results and discussion are presented in section 4. The main conclusions are summarized in section 5.

Data
MGP model assumes that a scene can be divided into four homogenous and isothermal components, that is, sunlit vegetation (C), shaded vegetation (T), sunlit soil (G) and shaded soil (Z), and the total radiation of this scene can be expressed as (Pinheiro et al., 2004): are the fractions of C, T, G, and Z, respectively; R c , R t , R g , and R z are the emission radiations of C, T, G, and Z, respectively; θ v , θ s , and φ are the viewing zenith angle, SZA and relative azimuth angle, respectively.
According to the Stefan-Boltzmann law, equation (1) can be transferred as follows: is the directional radiative temperature of the scene, T c , T t , T g , and T z are the component temperature of C, T, G, and Z, respectively; ε c , ε t , ε g , and ε z are the component emissivity of C, T, G, and Z, respectively.
The component fractions can be estimated by using the Boolean Scene Model (Serra, 1982), such that (Li & Strahler, 1992) where canopy is assumed as ellipsoid of which horizontal radius, vertical radius and height of crown center are r, b, and h, respectively; A v , A s , and O are the unit projected area from the θ v , θ s and their overlap area, respectively; ζ is the density of ellipsoid centers which can be calculated by using the following equation (Liu et al., 2004): From the equations (2)-(15), one can see that it is necessary to prescribe the component temperature, component emissivity, viewing/illuminating geometries, and structural parameters of canopy to simulate an extensive multiangle radiative temperature data set. The detailed schemes to generate this data set are listed below.
1. Component temperatures. Just like the illumination for the bidirectional reflectance, component temperature is the most important factor for TRD. Following previous studies, the relationships among four component temperatures are (a) the sunlit components are hotter than the shaded components Su et al., 2002); (b) the temperature difference between sunlit and shaded soil is higher than the temperature difference between sunlit and shaded vegetation, and the latter can be even neglected in pixel scale Ermida et al., 2014); (c) the sunlit soil has the highest temperature during daytime (Bian et al., 2018;Bian et al., 2016;Ermida et al., 2014;Su et al., 2002;Timmermans et al., 2009;Verhoef et al., 2007); and (d) the temperature of shaded soil may be higher than the temperature of sunlit vegetation Timmermans et al., 2009), may be between the temperature of sunlit vegetation and shaded vegetation Verhoef et al., 2007), and may be also lower than the temperature of shaded vegetation (Ermida et al., 2014;Verhoef et al., 2007). Therefore, component temperatures are generated according to the following four steps: The differences between T c and T t are from 1 to 5°C by a step of 2°C; c) The differences between T z and T c are from −5 to 10°C by a step of 5°C; d) The differences between T g and T z are from 5 to 15°C by a step of 5°C. As a result, this study has prescribed 36 component temperature profiles.

Component emissivities.
Previous studies indicate that the vegetation emissivity can be treated as a constant with the variation of view angles (Göttsche et al., 2013) whereas the soil emissivity presents an angular dependency Ren et al., 2011). García-Santos et al. (2012) revealed that soil emissivity decreases with the increase of θ v , and their relationship can be approximated as follows: where α is the parameter to express the strength of angular dependency. Two emissivity groups are set in this study. For high emissivity group, vegetation emissivity ε v and soil emissivity at nadir ε s (0) are set as 0.99 and 0.97, respectively; for low emissivity group, these two variables are set as 0.97 and 0.93, respectively (Bian et al., 2018). It is worth noting that ε c = ε t = ε v and ε g = ε z = ε s . Meanwhile, α is prescribed as three values, that is, 3.5, 3.3, and 3.1, to represent strong, middle, and weak angular effect, respectively (Ermida et al., 2018). Eventually, the total of six cases of component emissivity has been developed.
3. Viewing and illumination geometries. Limited by the valid swath angle of current satellite sensors such as Visible Infrared Imaging Radiometer Suite or Moderate Resolution Imaging Spectroradiometer (MODIS), θ v is set from 0°to 60°at 3°step and φ is set from 0°to 360°at 10°step. Meanwhile, this study considers three θ s , that is, 15°, 30°, and 45°. 4. Structure parameters. Same as the MODIS BRDF product , the ratios b/r representing crown shape parameter and h/b representing crown relative height, are preselected to 1 and 2, respectively. The values of PTC refer to Rasmussen et al. (2011) Earth and Space Science viewing zenith angles and 37 relative azimuth angles. By combining these parameters and MGP model, an extensive multiangle radiative temperature data set has been developed.

GO Kernels
Currently, the widely used GO kernel for BRDF model is K GOLSR , which is derived from GOMS model Wanner et al., 1995). Similarly, MGP model is also based on GOMS model. For these reasons, we develop the GO kernels for thermal radiation from MGP model.
For the sparse scene, one can approximate e x ≈ 1 + x. Then, equations (3)-(6) can be expressed as follows: If we substitute equations (17)-(20) into equation (1), there is In order to derivate the K GOLSR , there are two assumptions: Shaded components are perfectly black and the radiances of sunlit components are identical. Then, equation (21) can be expressed as follows: In consideration of the reciprocal property in the viewing and illumination geometries, the K GOLSR is Obviously, it is necessary to use observations from at least two angles to estimate R (θ v , θ s , φ) with K GOLSR .
Here, we infer the GO kernels for thermal radiation from the perspective of the component temperature differences. If one do not simplify the number of components, equation (21) can be described as follows: Then, the four-component GO kernel (K GO4 ) is where C 1 , C 2 , and C 3 are constants and C 1 ≥ 0,C 2 ≤ 0. In this case, there need measurements from at least four angles to estimate R (θ v , θ s , φ) with K GO4 .
Account for the fact that the temperature difference between sunlit and shaded soil is higher than the temperature difference between sunlit and shaded vegetation Ermida et al., 2014;Ermida et al., 2018), we can assume that R c = R t and equation (21) can be simplified as follows: 10.1029/2019EA000895

Earth and Space Science
Then, the three-component GO kernel (K GO3 ) is as follows: where C 1 and C 2 are constants and C 1 ≤ 0. In this case, the least requirement of angle number is 3.
We can further simplify components, such that only soil and vegetation components are considered, that is, R c = R t and R z = R g . Correspondingly, equation (21) becomes as follows: Then, the two-component GO kernel (K GO2 ) is as follows: where C is constant. Similar to K GOLSR , observations from at least two angles are needed.
Finally, based on the fact that the temperature of sunlit soil is highest, only R g is considered. Then, equation (21) is as follows: The empirical GO kernel with only considering sunlit soil component (K GOg ) is as follows: where C is constant and no more than 0. Clearly, the angular requirement is also 2.

Overall Performance
As is described in section 2, the generated data set has 3,888 cases and every case has 777 angle groups. Here, the difference between radiative temperature at off-nadir and radiative temperature at nadir is used to express the TRD . For every case, the root-mean-square errors (RMSEs) of fitting TRD for K GO4 , K GO3 , K GO2 , and K Gog have been separately calculated. In order to present their improvements clearly, the accuracies of K GOLSR have been estimated for comparison. Besides the RMSE, relative RMSE (RRMSE), the ratio between RMSE and TRD range, is also used to remove the influence of original magnitude of TRD.
The overall performance is shown in Figure 1. The average and maximum TRDs of these 3,888 cases are 2.3 and 7.5°C, respectively. In terms of K GO4 , the average and maximum RMSEs are 0.01 and 0.04°C, respectively; the average and maximum RRMSEs are 0.38% and 1.67%, respectively. Such high accuracies indicate that the approximation e x ≈ 1 + x makes hardly any influence. All red scatter points locate below the 1:1 line suggesting that the performances of K GO4 are better than the performances of K GO3 . This result can be expected since K GO4 consider all components. However, the simplification of sunlit and shaded vegetation cannot bring a significant influence. As the X axis shows, the maximum RMSE and RRMSE of K GO3 are 0.22°C and 8.99%, respectively. Further checking original data reveals that the average RMSE and RRMSE are 0.03°C and 1.24%, respectively. These results are satisfactory illustrating that the simplification of vegetation components is desirable. All of the blue points are above the 1:1 line, and most of them are below the green as well as magenta points, which indicates that K GOg is worse than K GO3 but performs best under the same angular requirements, that is, two angles. The average and maximum RMSEs of K GOg are 0.06 and 0.37°C, respectively; the average and maximum RRMSEs are 3.34% and 18.05%, respectively. These results are also satisfactory indicating that the way only considering of the hottest component is feasible. With regard to K GO2 , its performances are lower than K GOg , whereas their changes are relatively stable. The average and maximum RMSEs of K GO2 are 0.25 and 0.81°C, respectively; the average and maximum RRMSEs are 11.28% and 21.05%, respectively. By contrast, RMSEs of K GOLSR have a significant change range of which they can be better than K GO3 but can be also worse than K GO2 . The average and maximum RMSEs of K GOLSR are 0.25 and 1.08°C, respectively; the average and maximum RRMSEs are 10.66% and 23.41%, respectively.
Further examination reveals that there are six ranking types in the performances of these five GO kernels. As the Table 1 shows, K GO4 has the highest accuracies for all 3,888 cases. Overall, in addition to the case of B3G2, K GO3 ranks second with a proportion of 94.6%. K GOg performs the third best and the proportion, that is, 3GB2 and 3G2B, is 90.25%. The performance of K GOLSR fluctuates greatly, it can rank second with a ratio of 5.40%; it can also rank third with a ratio of 3.85%; but in most cases it ranks fourth with a ratio of 45.37%. K GO2 performs worst with a fifth-place ratio of 54.47%.

10.1029/2019EA000895
Earth and Space Science Figure 2 is the fitting results of five GO kernels under different component temperature differences. Clearly, the Group 3GB2 and the Group 3G2B account for most proportions. An obvious phenomenon is that K GOLSR is better than K GO2 when T c ≥ T z and vice versa. In terms of T z − T c = −5 and 0°C, the accuracies of K GOLSR become higher along with the increase of T c − T t or the decrease of T g − T z . In particular, for T z − T c = −5°C and T g − T z = 5°C, K GOLSR is better than K GOg if T c − T t < 5°C and K GOLSR is even better than K GO3 if T c − T t = 5°C. These results indicate that the more significant the effect of sunlit vegetation, the better the performance of K GOLSR . From the perspective of thermal radiation, the contribution of sunlit soil is greatest. If T c ≥ T z , as the increases of T c , the distribution of component temperatures is more close to the assumptions of K GOLSR ; that is, only two sunlit components make contributions and their contributions are identical. Therefore, the above result can be expected. In the case of T z − T c = 5 and 10°C, the performances of K GO2 become better with the decrease of T g − T z or T c − T t . This phenomenon can also be expected since a small temperature difference between sunlit and shaded components means they can be simplified as one component which is in accordance with the hypothesis of K GO2 . Figure 3 presents the result of five GO kernels for different component emissivities. The type of emissivity mainly has an influence on the mutual relationship between the performances K GOLSR and K GO2 . First, as the increase of angular effect, the counts of 3GB2 would fall and the counts of 3G2B would grow. That is to say, K GO2 becomes better, whereas K GOLSR becomes worse. The reason may be that, in the case of more significant angular dependency of component emissivity, the effect of soil becomes more important or the effect of vegetation become slighter. As a result, it is more difficult to satisfy the assumption of K GOLSR that the contributions of sunlit soil and sunlit vegetation are identical resulting in a worse performance. On the other hand, the ε v − ε s (0) is 0.02 for high emissivity group and is 0.04 for low emissivity group. Therefore, there is a relatively strong soil effect for high emissivity group and a relatively strong vegetation effect for low emissivity group. For the same reason, the performance of K GO2 is better than the performance of K GOLSR for the high emissivity group whereas the accuracy of K GOLSR is better than the accuracy of K GO2 for the low emissivity group.

Earth and Space Science
Similarly, the main influence of SZA is also the mutual relationship between the performances of K GOLSR and K GO2 . As Figure 4 shows, with the increase of SZA, the accuracy of K GO2 improves, whereas the accuracy of K GOLSR decreases. Specifically, as SZA ranges from 15°to 30°to 45°, the number of 3G2B ranges from 372 to 679 to 694, whereas the number of 3GB2 ranges from 804 to 503 to 457.
The fitting results of five GO kernels for different PTCs are shown in Figure 5. Along with the increase of PTC, the counts of 3G2B are ever decreasing.

10.1029/2019EA000895
Earth and Space Science

Determination of Best Structure Parameters
In above subsections, the b/r and h/b are set as 1 and 2, which equal to the values of simulated data set. Actually, the b/r and h/b have different values for different scenes resulting in district performances even for the same kernels (Wanner et al., 1995). That is to say, if one uses the uniform selection, for example, 1 and 2, it may bring some errors. Therefore, the determination of their best values is very meaningful for practical application of GO kernels. Combing the schemes of previous studies (Table 2), this study set 32 groups of b/r and h/b representing a variety of situations. Specifically, the values of b/r are 0.5, 0.75, 1, 1.25, 2, and 2.5 and the h/b is from 1.5 to 3.5 by a step of 0.5. Moreover, two special combinations, that is, one is 2 + 1 for shrub scene and the other one is 3.5 + 2 for conifer scene, have also been developed. By using the same parameters described in section 2 and upgraded b/r and h/b, a more comprehensive data set with 124,416 cases, that is, every group of b/r and h/b has 3,888 cases, has been generated. Then, the accuracies of 30 settings except for shrub and conifer scenes have been separately estimated based on the above data set. The best preselection means that using it can generate the most robust performance. For this reason, the average RMSE reflecting the overall accuracy and the largest RMSE reflecting the extreme accuracy of 124,416 cases are chosen as the evaluation indicators.  The results are presented in Figure 6. From the perspective of RMSE mean, if one uses uniform parameters, the performance ranking is K GO4 , K GO3 , K GOg , K GO2 , and K GOLSR . If b/r and h/b are known, as analyzed in section 4.1, the order is K GO4 , K GO3 , K GOg , K GOLSR , and K GO2 . This slight difference indicates that the K GOLSR is more sensitive about the values of b/r and h/b. With respect to RMSE maximum, the rank list would become K GO4 , K GO3 , K GO2 , K GOg , and K GOLSR , suggesting that K GOg is less sensitive than K GOLSR but more sensitive than other three kernels. Overall, the influence mode of these two structure characteristics is different for different GO kernels. For K GO4 and K GO3 , the main influence factor of RMSE mean is h/b, whereas RMSE maximum is more affected by b/r. K GO2 only contains b/r, and therefore, its performances have no connection with h/b. In terms of K GOg and K GOLSR , both b/r and h/b have an important influence but the effect of b/r is more significant not only for RMSE mean but also for RMSE maximum.
Regarding to K GO4 , the best choices of b/r and h/b are 2.5 and 2, respectively. This group would bring a relatively small overall error and extreme error of which RMSE mean is 0.08°C and RMSE maximum is 0.74°C. Concerning with K GO3 , the best values of b/r and h/b are 2.5 and 2, respectively. In this case, RMSE mean is 0.11°C and RMSE maximum is 0.92°C. For K GOg , its overall accuracies would improve with the increases of b/r or the decreases of h/b. When b/r ≥ 2, its extreme errors would be smaller with the increases of h/b. For these reasons, in aggregate, the best preselections of b/r and h/b are 2.5 and 1.5, respectively. Correspondingly, RMSE mean and RMSE maximum are 0.16 and 1.14°C, respectively. As the increases of b/r, the performances of K GO2 become ever better. Therefore, the best scheme of K GO2 is that b/r takes the maximum, that is, 2.5. Correspondingly, the smallest RMSE mean and RMSE maximum are 0.24 and 1.15°C, respectively. With reference to K GOLSR , RMSE mean and maximum would enhance with the increase of b/r and h/b for almost all situations. In general, the best combination of b/r and h/b is 0.5 and 2 of which the RMSE mean is 0.32°C and RMSE maximum is 2.96°C.

Conclusion
Based on GO model, this study derived GO kernels for thermal radiation from the perspective of the component temperature difference rather than the illumination difference. As a result, four GO kernels have been developed. Their performances as well as the comparison with a widely used GO kernel in BRDF model, that is, K GOLSR , have been analyzed in detail by using an extensive data set. The main conclusions were summarized below.
1. In the case of sufficient viewing angles, K GO4 performed best and then was K GO3 . If only two angles were available, the overall accuracy of K GOg was highest; though the result of K GOLSR was slightly better the result of K GO2 , it could bring a bigger extreme error.