# Determination of Optimum Viewing Angles for the Angular

## Transcript Of Determination of Optimum Viewing Angles for the Angular

Sensors 2015, 15, 7537-7570; doi:10.3390/s150407537 Article

OPEN ACCESS

sensors

ISSN 1424-8220 www.mdpi.com/journal/sensors

Determination of Optimum Viewing Angles for the Angular Normalization of Land Surface Temperature over Vegetated Surface

Huazhong Ren 1,2,4, Guangjian Yan 1,*, Rongyuan Liu 1,4,*, Zhao-Liang Li 3,4, Qiming Qin 2, Françoise Nerry 4 and Qiang Liu 5

1 State Key Laboratory of Remote Sensing Science, School of Geography, Beijing Normal University, Beijing 100875, China; E-Mail: [email protected]

2 Institute of Remote Sensing and Geographic Information System, Peking University, Beijing 100871, China; E-Mail: [email protected]

3 Key Laboratory of Agri-Informatics, Ministry of Agriculture/Institute of Agricultural Resources and Regional Planning, Chinese Academy of Agricultural Sciences, Beijing 100081, China; E-Mail: [email protected]

4 ICube Laboratory, Universitéde Strasbourg, 67412 Illkirch, France; E-Mail: [email protected] 5 College of Global Change and Earth System Science, Beijing Normal University, Beijing 100875,

China; E-Mail: [email protected]

* Authors to whom correspondence should be addressed; E-Mails: [email protected] (G.Y.); [email protected] (R.L.); Tel.: +86-10-5880-2085 (G.Y.).

Academic Editor: Fabrizio Lamberti

Received: 30 January 2015 / Accepted: 23 March 2015 / Published: 27 March 2015

Abstract: Multi-angular observation of land surface thermal radiation is considered to be a promising method of performing the angular normalization of land surface temperature (LST) retrieved from remote sensing data. This paper focuses on an investigation of the minimum requirements of viewing angles to perform such normalizations on LST. The normally kernel-driven bi-directional reflectance distribution function (BRDF) is first extended to the thermal infrared (TIR) domain as TIR-BRDF model, and its uncertainty is shown to be less than 0.3 K when used to fit the hemispheric directional thermal radiation. A local optimum three-angle combination is found and verified using the TIR-BRDF model based on two patterns: the single-point pattern and the linear-array pattern. The TIR-BRDF is applied to an airborne multi-angular dataset to retrieve LST at nadir (Te-nadir) from different viewing directions, and the results show that this model can

Sensors 2015, 15

7538

obtain reliable Te-nadir from 3 to 4 directional observations with large angle intervals, thus corresponding to large temperature angular variations. The Te-nadir is generally larger than temperature of the slant direction, with a difference of approximately 0.5~2.0 K for vegetated pixels and up to several Kelvins for non-vegetated pixels. The findings of this paper will facilitate the future development of multi-angular thermal infrared sensors.

Keywords: BRDF; angular normalization; land surface temperature; multi-angular; WiDAS

1. Introduction

Land surface temperature (LST) is required in many applications, including agrometeorology, climate and environmental studies [1,2]. Thermal infrared images from aircraft and spaceborne satellites provide a unique opportunity to map this parameter at regional and even global scales. Many methods have been proposed to retrieve LST data from remotely sensed data using different specifications of thermal infrared sensors and the atmospheric and emissivity data situations [3–5]. These methods can been roughly grouped into three categories: single-channel algorithms [6,7], multi-channel methods (e.g., the split-window algorithm [8,9] and the temperature and emissivity separation method [10,11]), and multi-time methods (e.g., the temperature-independent spectral index (TISI) method [12,13], the two-temperature method (TTM) [14], and the physical day/night algorithm [15]). For operational purposes, these methods often take the observed pixel as a homogeneous and isothermal target; this assumption is reasonable for pure or quasi-pure pixels such as bare soil, sand, snow and dense vegetated surfaces. However, for mixed pixels including two or more components at different temperatures and emissivities, their pixel temperature exhibits spectral and angular variations. As a result, the above assumption will be incorrect, and the retrieved LST will only present, at the least, the effective temperature at its corresponding viewing direction and cannot be directly taken as the temperature at nadir or under other directions in theory.

The angular behavior of LST has been investigated in many previous studies [16–22], and this angular variation results primarily from the angular variation in the pixel emissivity for three-dimensional surfaces and the relative weights of more than one component (e.g., leaf and background soil) with different temperatures included in the scene. Certain ground measurements have indicated that the LST difference at nadir and off-nadir observations can be as large as 5 K for bare soils and up to 10 K for urban areas [23]. For satellite images, the pixels also face a similar situation because the pixels in the same image are sometimes observed at different viewing angles, especially by those sensors with a large field of view (FOV). For example, the Moderate Resolution Imaging Spectroradiometer (MODIS) scans the land surface in the cross-track direction with a viewing zenith angle (VZA) varying from −65°to +65°[21,24]. Thus, angle-dependent variations in the retrieved LST are inevitable, which make the LSTs of different pixels in the same image incomparable and eventually lead to an error up to about 2.5 K [21]. Similar cases can be observed in other satellite sensors such as the Advanced Very High Resolution Radiometer (AVHRR) and Spinning Enhanced Visible and InfraRed Imager (SEVIRI) instruments. Therefore, it is crucial to perform angular

Sensors 2015, 15

7539

normalization of the LST data, i.e., to translate LST data obtained in different directions to a specified direction such as the nadir direction [3,25].

Multi-angular observation of the same target is considered to be the most promising method of solving this problem, and such observations can be achieved using ground goniometers: the Portable Apparatus for Rapid Acquisitions of Bi-directional Observations of Land and Atmosphere (PARABOLA) [26], the field goniometer system (FIGOS) [27], the portable Multi-Angle Observation System (MAOS) [23,28], and others [29,30]. However, only a small number of studies have been published on LST angular normalization or on the simultaneous retrieval of directional emissivity and temperature from multi-angular TIR images mainly due to a lack of multi-angular observations at satellite level. To date, only ATSR (Advanced Along Track Scanning Radiometer) series satellites [31,32] have provided multi-angular observations of thermal infrared (TIR) data because of the difficulty and the complexity of manufacturing highly controlled multi-angular sensors and also because no studies that discuss requirements in terms of the number of direction specifications of angular observations have been published. To address this situation, we investigate the minimum requirements of viewing observations and directions for achieving angular normalization of LSTs using multi-angular TIR data to provide technique support for the future design of such sensors. This paper is organized as follows: Section 2 will introduce a new model to simulate the directional thermal radiance and discuss the performance of the kernel-driven bi-directional reflectance distribution function (BRDF) model in the TIR domain as a TIR-BRDF model. Section 3 will investigate the local optimum three-angle combinations for the TIR-BRDF using a single-point pattern and discuss the three-angle combinations using a linear array pattern that is always used in airborne and satellite sensors. Section 4 will apply the TIR-BRDF model to an airborne multi-angular TIR image from the WiDAS (Wide-angle infrared Dual-mode line/area Array Scanner) system [22,33] in the WATER (Watershed Allied Telemetry Experimental Research) campaign [34]. Sections 5 and 6 will provide discussions and conclusions of this paper. The full names and the corresponding abbreviations of some terms are listed in the Table A1 of Appendix.

2. Modeling of Directional Thermal Radiation

Modeling the directional thermal radiance (DTR) of homogenous or heterogeneous canopies has attracted a substantial amount of attention. The main reason for the DTR is that the fractions of different components with different temperature and emissivity vary with the changes of viewing angles. Based on this concept, considering N different components in a pixel under clear-sky conditions, a generalized parameterization of the DTR can be expressed as [22]:

L(v , s , ) B[DBT (v , s , )] (v , v )B[Te (v , s , )]

N

fi (v , s , ) i Bi (Ti ) Lmulti

(1)

i1

where θv and θs are the viewing zenith angle (VZA) and the solar zenith angle (SZA), respectively; φ is the relative azimuth angle (RAA) between the viewing azimuth angle (VAA) φv and the solar azimuth angle (SAA) φs; L(θv, θs, φ) is the directional thermal radiance (the inversion of this term according to the Planck’s law B[] will produce the directional brightness temperature (DBT)); Meanwhile, ε(θv,φv) is the pixel ensemble directional emissivity and Te(θv, θs, φ) is the directional

Sensors 2015, 15

7540

effective temperature; fi are the fractions of various components, such as the commonly used sunlit and shadowed leaves and sunlit and shadowed soil; Ti and εi are the temperatures (unit: K) and emissivities of those components, and note that εi are assumed to be independent of viewing angle; and Lmulti is the multi-scattering between soil and leaves and between leaves inside the canopy.

We consider the homogenous canopy and estimate the component fractions fi from a parameterization [35] of the Scattering by Arbitrarily Inclined Leaves (SAIL) Hotspot (SAILH) model [36,37] using the solar and viewing directions, canopy structure and model parameters. In addition, only three components (leaves and sunlit and shadowed soils) are considered because the temperature difference between sunlit and shadowed leaves is very small compared to the temperature difference between sunlit and shadowed soils. Furthermore, we calculate the multi-scattering term Lmulti with the following equation:

Lmulti (1 M )Lleaf (1 g ) (1 )[1 b()M ][1 b()](1 v) Lleaf

(2)

The first part of the right-hand side of Equation (2) is the downward leaves’ radiation reflected by the soil, and the second part is the multi-scattering interior to the canopy; εg and εv are the emissivities of ground soil and leaf, respectively; b(θ) is the directional gap frequency in the θ direction; and M is the hemispheric gap probability of the canopy. For a canopy with a spherical leaf angle distribution and random dispersion, b(θ) and M can be expressed as Equation (3), in which LAI is the leaf area index. α denotes the cavity effect that accounts for multiple scatterings inside the canopy, the details of which can be found in [38].

b() exp[ 0.5

LAI ] , and M 1

2

b()d exp(0.825LAI )

(3)

cos

2

Table 1 lists the main input variables of the simulation DBT, including the solar zenith and azimuth angles; LAI of the canopy; the temperatures of the leaves (Tleaf), sunlit soil (Tsun_soil) and shadowed soil (Tshd_soil); and the emissivities of leaf (εv) and soil (εg). Note that the simulated DBT is assumed to be atmospherically corrected.

Table 1. The main input parameters for the directional brightness temperature (DBT) simulation.

SAA SZA

LAI

Tleaf

Tsun_soil

Tshd_soil

εv

εg

120° 30° 0.5~5.0 with a step of 0.5 305 K 320 K 315 K 0.985

0.95

In addition, we also used the kernel-driven BRDF model, which is always used to link the bi-directional reflectance to the viewing geometry in the visible/near-infrared wavelengths and to link the DTR and DBT to the viewing geometry by replacing the bi-directional reflectivity with the DTR [22,39]. As a result, the kernel-driven BRDF in the TIR field is written as Equation (4). For simplicity, the new kernel-driven BRDF model is henceforth called TIR-BRDF.

DTR(v,s ,) B[DBT(v,s ,)] fiso fvol kvol (v, s , ) fgeo kgeo (v,s ,)

(4)

where fiso is the isotropic scattering term, fvol is the coefficient of the volumetric kernel kvol, and fgeo is the coefficient of the geometric kernel kgeo. The Ross-Thick kernel developed in [40] and the Li-SparseR kernel derived in [41] on the basis of the geometric-optical mutual shadowing BRDF model [42] are used as the volumetric and geometrical scattering kernels, respectively. A suitable

Sensors 2015, 15

7541

expression of kvol was derived by Roujean et al. [40], called the Ross-Thick kernel for its assumption of a dense leaf canopy (see Equation (5)). It is a single-scattering solution of the radiative transfer equation for plane-parallel dense leaf canopy with uniform leaf angle distribution, a Lambertian background and equal leaf transmittance and reflectivity, but it does not account for the hotspot effect.

( / 2 ) cos sin

kvol (v , s , ) cos cos 4

(5)

v

s

where, ξ is the phase angle, related to the sun-target-observer position as:

cos v cos s sin v sin s cos

(6)

The Li-Sparse kernel derived from the geometric-optical mutual shadowing BRDF model [42] was reported to work well with the observed data. The original form of this kernel is not reciprocal in the viewing and solar directions, a property that is expected from homogeneous natural surface viewed at coarser spatial resolution, but then was refined to be reciprocal by assuming the sunlit component in the viewed scene simply varies as 1/cosθs. As a result, the reciprocal model was the Li-SparseR kernel [43,44]:

kgeo (v , s , ) O(v , s , ) sec 'v sec 's 1 (1 cos ')sec 'v sec 's (7) 2

with,

O 1 (t sin t cos t)(sec 'v sec 's )

cos t h b

D2 (tan 'v tan 's sin )2 sec 'v sec 's

D tan2 'v tan2 's 2 tan 'v tan 's cos

cos ' cos 'v cos 's sin 'v sin 's cos

' tan1(b tan ) , ' tan1(b tan ) ,

v

r

v

s

r

s

where, O is the overlap area between the view and solar shadows. The term cost should be constrained to the range (−1, 1), as values outsides of this range imply no overlap and should be disregarded. The ratio h/b and b/r are the dimensionless crown relative height and shape parameters, respectively, and should be preselected. In this thesis, h/b = 2 and b/r = 1, i.e., the spherical crowns are separated from ground by half their diameter.

Figure 1 shows the simulated DBT distribution (the first column), the fitted DBT from the TIR-BRDF model (the second column), and their temperature difference (the third column) at (a) LAI = 0.5; (b) LAI = 1.0 and (c) LAI = 2.0 based on the input variables shown in Table 1 for the DBT simulation using Equation (1) and the SAILH parameterization model [35]. The maximum zenith angle is constrained to 60°because the kernel-driven BRDF model was reported to obtain unacceptable results for larger zenith angles [38,45].

Sensors 2015, 15

7542

(a) LAI = 0.5

(b) LAI = 1.0

(c) LAI = 2.0 Figure 1. The hemispheric distribution of the simulated DBT from Equation (1) (the first column), the fitted DBT from the TIR-BRDF model (the second column) and their temperature difference (the third column) at (a) LAI = 0.5; (b) LAI = 1.0 and (c) LAI = 2.0, respectively. As observed from Figure 1, the angular variation in the DBT can be as high as approximately 3.0 K under the simulation conditions (Table 1). However, this angular variation is based on the components’ temperature difference and canopy structure: larger component temperature differences generally lead to larger angular variations. The TIR-BRDF model produces a good fit of the DBT distribution, with errors lower than 0.3 K. The hotspot effect in the solar direction is significant; however, we also find that the maximum temperature difference occurs around this direction and that the DBT of the TIR-BRDF model is generally larger than that of the simulation result in the directions around the solar beam. The increase in LAI decreases the DBT because the fraction of leaves increases and the leaves’ temperature is lower than that of the soils. Furthermore, the LAI also influences the temperature difference caused by the regression used in the TIR-BRDF model. Figure 2 shows the RMSE (root-mean-square error) and the maximum temperature difference as a function of LAI. The figure

Sensors 2015, 15

7543

indicates that the RMSE is smaller than 0.1 K and that the maximum temperature is smaller than 0.3 K; both obtain their global maximums at approximately LAI = 2.

Figure 2. The influence of LAI on RMSE and maximum temperature difference from the TIR-BRDF model.

Furthermore, we also implemented the TIR-BRDF model using different solar positions, component temperatures and emissivities, and canopy structures, and the results further confirmed the performance of the TIR-BRDF model for calculating the DBT distribution. Note that because the component emissivities are considered to be angle independent and because, consequently, the DBT of a bare surface is isotropic, we only discuss at the level of the canopy rather than for the case of a bare surface.

3. Viewing Angle Specification for Angular Normalization of LST

The above result is obtained using the TIR-BRDF model with all viewing zenith and azimuth angles. However, in practice, it is almost impossible to observe the same target at numerous directions. On the contrary, only some angular observations are usually conducted. Consequently, a problem arises: how many observations and under what viewing angles are needed for the TIR-BRDF model to accurately fit the DBT and consequently enable the angular normalization of the LST? To determine the optimum observation groups, this paper used two different patterns in the following discussions: single-point patterns and linear-array patterns.

3.1. Single-Point Pattern Analysis for Multiple-Viewing Angle Specification

According to Equation (4), at least three angular observations are needed to drive the TIR-BRDF model. Additional observations can improve the reliability of the fitted coefficients in theory, but the errors in those observations may produce greater uncertainty in the fitted result. Therefore, we start with only three viewing angles to determine the local optimum viewing angle combination that produces the smallest error in the TIR-BRDF model and to provide various suggestions for the future development of multi-angular airborne or space-borne TIR sensors. Therefore, to better facilitate the

Sensors 2015, 15

7544

development of a mechanical design of such a sensor, we herein assume that all observations are performed with the azimuths in the same plane, i.e., their azimuths equal φ or φ + 180°.

Taking φ = 0° for example, the three observations’ zenith angles vary in the azimuth plane (0°~180°) at a step of 10°to a maximum of 60°in the zenith direction. Different observations are required to have different value of zenith or azimuth angles. As a result, there are 126 groups with three angles. To determine the optimum group, we first use Equation (1) to simulate hemispheric DBTs under varying solar positions (SAA from 0°to 330°with a step of 30°, SZA = 10°, 30°and 50°), LAIs (from 0.5 to 5.0 with a step of 0.5) and the component temperatures and emissivities in Table 1, thus resulting in a total of 396 cases. Consequently, the three coefficients in Equation (4) are obtained by regressing the DBT with the known kernels (kvol and kgeo) calculated in different viewing geometry, and this process is additionally controlled by an optimization algorithm on the basis of that the brightness temperature at the hotspot direction (i.e., solar beam direction) has the large value because this direction corresponds to the largest proportions of sunlit soil and leaves which have a higher temperature. Table 2 shows the frequency of the temperature RMSEs in the ranges of [0.0, 0.5] K and [0.5, 1.0] K and the number of times the maximum and minimum RMSEs are obtained for different angle groups. Furthermore, because it was found that the viewing angle group with nadir viewing direction (VZA = 0° and VAA = 0°) had better results than those without such nadir observation, and also because of space limitations of this paper, we did not display those results without the nadir observation.

Furthermore, four criteria are proposed here to filter the final optimum angle combinations:

(1) Most of the RMSEs in the TIR-BRDF model should fall in the range of 0.0~0.5 K. (2) The maximum difference should not or seldom occur for acceptable combinations of angles. (3) Most of the minimum differences can be obtained using the combinations. (4) Angle difference of the adjacent viewing angles should be large enough in order to cause large

difference in the component fractions and then lead a large angular variation in the observed brightness temperature. But a very large VZA are not recommended in order to avoid large pixel size differences in the angle combination.

Table 2. The frequency of the root-mean-square error (RMSE) in different three-direction combinations. Columns 1 and 2 are the VAA and VZA for the first direction observation, and columns 3 and 4 are for the second direction observation, and columns 5 and 6 are for the third direction observation.

1st VAA

0 0 0 0 0 0 0

1st VZA

0 0 0 0 0 0 0

2nd VAA

0 0 0 0 0 0 0

2nd VZA

10 10 10 10 10 10 20

3rd VAA 180 180 180 180 180 180 180

3rd VZA

10 20 30 40 50 60 10

RMSE [0.0~0.5] K

195 241 299 316 337 323 243

RMSE [0.5~1.0] K

44 43 40 33 7 4 42

Max *

122 49 16 15 1 5 45

Min #

0 1 2 0 12 11 1

Sensors 2015, 15

7545

Table 2. Cont.

1st 1st 2nd 2nd 3rd 3rd VAA VZA VAA VZA VAA VZA

RMSE (0.0~0.5) K

RMSE (0.5~1.0) K

Max *

Min #

0

0

0 20 180 20

295

18

11

0

0

0

0 20 180 30

315

34

14

4

0

0

0 20 180 40

337

20

2

3

0

0

0 20 180 50

359

3

1

15

0

0

0 20 180 60

354

7

19

11

0

0

0 30 180 10

299

39

17

2

0

0

0 40 180 10

316

35

16

0

0

0

0 40 180 20

337

20

2

3

0

0

0 40 180 30

361

2

0

6

0

0

0 40 180 40

364

1

0

3

0

0

0 40 180 50

373

0

0

26

0

0

0 40 180 60

373

0

1

47

0

0

0 50 180 10

336

8

1

14

0

0

0 50 180 20

359

4

1

16

0

0

0 50 180 30

360

0

0

14

0

0

0 50 180 40

373

0

0

26

0

0

0 50 180 50

353

2

2

14

0

0

0 50 180 60

365

8

5

69

0

0

0 60 180 10

322

5

4

10

0

0

0 60 180 20

354

8

22

12

0

0

0 60 180 30

373

0

1

10

0

0

0 60 180 40

373

0

1

47

0

0

0 60 180 50

365

8

3

69

0

0

0 60 180 60

349

0

21

0

* Max stands for the frequency of the maximum temperature difference, while # Min stands for the frequency

of the minimum temperature difference.

According to the above four criteria and the results shown in Table 2, we find that large VZAs with respect to the nadir direction can reduce the fitting error, and several angle combinations of VAA and VZA can be considered as candidates for the optimum groups: ① [(0°, 0°), (0°, 30°), (180°, 50°)]; ② [(0°, 0°), (0°, 40°), (180°, 50°)]; and ③ [(0°, 0°), (0°, 40°), (180°, 60°)] (see highlight in Table 2). Although the combination [(0°, 0°), (0°, 50°), (180°, 60°)] produces the minimum error most often, the number of times the maximum error is obtained (five times) and the large viewing zenith angles make this group impractical.

Among the three candidates, group ③ cannot be recommended because of its large VZA in the direction (180°, 60°), whose pixel sizes are up to four times the pixel sizes at nadir observation if the sensor’s IFOV (instantaneous field of view) remains constant in every direction. Group ② performs slightly better than group ①. However, the final decision cannot be made without a sensitivity analysis of the errors in the observed data. To investigate the sensitivity of the different angle groups to the DBT errors, artificial temperature noises within [−0.5, 0.5] K and [−1.0, 1.0] K with a uniform distribution are added to the simulated DBT of Equation (1). The three coefficients of the TIR-BRDF model are recalculated and subsequently used to obtain the hemispheric DBT. The output of the

Sensors 2015, 15

7546

TIR-BRDF model is constrained whereby the DBT must reach its maximum value in the solar direction because the sunlit components in this direction have the largest percentage in this hotspot direction.

Figure 3 shows the histograms of the temperature RMSEs for the three-angle TIR-BRDF model using the angle groups ①–③ with respect to simulated random measurement error in [−0.5, 0.5] K and [−1.0, 1.0] K included in the input DBT data. The figure illustrates that group ③ produces the smallest errors under the two noise conditions, followed by groups ② and ①. The RMSE percentages of groups ①–③ in the range 0.0~1.0 K are approximately 92.4%, 96.6% and 96.3% for the noise level [−0.5, 0.5] K and 80.1%, 88.5% and 91.1% for the noise level [−1.0, 1.0] K. These results indicate that the three groups can produce a temperature error of less than 1.0 K for most cases if the noise involved in the observed DBT data is no more than 1.0 K. However, it is obvious that groups ② and ③ are less sensitive to the noise in the observed data. As stated above, group ③ cannot be used as the optimum viewing angle combination because of the large viewing zenith angle and the occurrence of maximum errors. As a result, group ② should be the optimum group in theory. However, because the pixel size increases with increasing VZA, the ratios of pixel area at VZA = 30°, 40°and 50°to that of nadir are 1.32, 1.70 and 2.40, respectively. It is necessary to use a pixel with a middle area to connect the observations from nadir to large VZAs, especially for heterogeneous surfaces. In this case, the VZAs in group ① produce a more continuous series of pixel sizes (1, 1.3 and 2.4 times the nadir pixel size). In addition, the angle intervals in the slant direction of group ① (i.e., 20°) are larger than those of group ② (i.e., 10°), which makes group ① less sensitive to the errors in the observation angle in theory. Therefore, according to the single-pattern analysis, we prefer the angle combination of group ① and take this group as the local optimum angle combination for the DBT regression using the three-angle TIR-BRDF model.

Figure 3. Histograms of the temperature RMSE caused by the three-angle TIR-BRDF model with respect to temperature noise about (a) [−0.5, 0.5] K and (b) [−1.0, 1.0] K included in the input DBT data, respectively.

OPEN ACCESS

sensors

ISSN 1424-8220 www.mdpi.com/journal/sensors

Determination of Optimum Viewing Angles for the Angular Normalization of Land Surface Temperature over Vegetated Surface

Huazhong Ren 1,2,4, Guangjian Yan 1,*, Rongyuan Liu 1,4,*, Zhao-Liang Li 3,4, Qiming Qin 2, Françoise Nerry 4 and Qiang Liu 5

1 State Key Laboratory of Remote Sensing Science, School of Geography, Beijing Normal University, Beijing 100875, China; E-Mail: [email protected]

2 Institute of Remote Sensing and Geographic Information System, Peking University, Beijing 100871, China; E-Mail: [email protected]

3 Key Laboratory of Agri-Informatics, Ministry of Agriculture/Institute of Agricultural Resources and Regional Planning, Chinese Academy of Agricultural Sciences, Beijing 100081, China; E-Mail: [email protected]

4 ICube Laboratory, Universitéde Strasbourg, 67412 Illkirch, France; E-Mail: [email protected] 5 College of Global Change and Earth System Science, Beijing Normal University, Beijing 100875,

China; E-Mail: [email protected]

* Authors to whom correspondence should be addressed; E-Mails: [email protected] (G.Y.); [email protected] (R.L.); Tel.: +86-10-5880-2085 (G.Y.).

Academic Editor: Fabrizio Lamberti

Received: 30 January 2015 / Accepted: 23 March 2015 / Published: 27 March 2015

Abstract: Multi-angular observation of land surface thermal radiation is considered to be a promising method of performing the angular normalization of land surface temperature (LST) retrieved from remote sensing data. This paper focuses on an investigation of the minimum requirements of viewing angles to perform such normalizations on LST. The normally kernel-driven bi-directional reflectance distribution function (BRDF) is first extended to the thermal infrared (TIR) domain as TIR-BRDF model, and its uncertainty is shown to be less than 0.3 K when used to fit the hemispheric directional thermal radiation. A local optimum three-angle combination is found and verified using the TIR-BRDF model based on two patterns: the single-point pattern and the linear-array pattern. The TIR-BRDF is applied to an airborne multi-angular dataset to retrieve LST at nadir (Te-nadir) from different viewing directions, and the results show that this model can

Sensors 2015, 15

7538

obtain reliable Te-nadir from 3 to 4 directional observations with large angle intervals, thus corresponding to large temperature angular variations. The Te-nadir is generally larger than temperature of the slant direction, with a difference of approximately 0.5~2.0 K for vegetated pixels and up to several Kelvins for non-vegetated pixels. The findings of this paper will facilitate the future development of multi-angular thermal infrared sensors.

Keywords: BRDF; angular normalization; land surface temperature; multi-angular; WiDAS

1. Introduction

Land surface temperature (LST) is required in many applications, including agrometeorology, climate and environmental studies [1,2]. Thermal infrared images from aircraft and spaceborne satellites provide a unique opportunity to map this parameter at regional and even global scales. Many methods have been proposed to retrieve LST data from remotely sensed data using different specifications of thermal infrared sensors and the atmospheric and emissivity data situations [3–5]. These methods can been roughly grouped into three categories: single-channel algorithms [6,7], multi-channel methods (e.g., the split-window algorithm [8,9] and the temperature and emissivity separation method [10,11]), and multi-time methods (e.g., the temperature-independent spectral index (TISI) method [12,13], the two-temperature method (TTM) [14], and the physical day/night algorithm [15]). For operational purposes, these methods often take the observed pixel as a homogeneous and isothermal target; this assumption is reasonable for pure or quasi-pure pixels such as bare soil, sand, snow and dense vegetated surfaces. However, for mixed pixels including two or more components at different temperatures and emissivities, their pixel temperature exhibits spectral and angular variations. As a result, the above assumption will be incorrect, and the retrieved LST will only present, at the least, the effective temperature at its corresponding viewing direction and cannot be directly taken as the temperature at nadir or under other directions in theory.

The angular behavior of LST has been investigated in many previous studies [16–22], and this angular variation results primarily from the angular variation in the pixel emissivity for three-dimensional surfaces and the relative weights of more than one component (e.g., leaf and background soil) with different temperatures included in the scene. Certain ground measurements have indicated that the LST difference at nadir and off-nadir observations can be as large as 5 K for bare soils and up to 10 K for urban areas [23]. For satellite images, the pixels also face a similar situation because the pixels in the same image are sometimes observed at different viewing angles, especially by those sensors with a large field of view (FOV). For example, the Moderate Resolution Imaging Spectroradiometer (MODIS) scans the land surface in the cross-track direction with a viewing zenith angle (VZA) varying from −65°to +65°[21,24]. Thus, angle-dependent variations in the retrieved LST are inevitable, which make the LSTs of different pixels in the same image incomparable and eventually lead to an error up to about 2.5 K [21]. Similar cases can be observed in other satellite sensors such as the Advanced Very High Resolution Radiometer (AVHRR) and Spinning Enhanced Visible and InfraRed Imager (SEVIRI) instruments. Therefore, it is crucial to perform angular

Sensors 2015, 15

7539

normalization of the LST data, i.e., to translate LST data obtained in different directions to a specified direction such as the nadir direction [3,25].

Multi-angular observation of the same target is considered to be the most promising method of solving this problem, and such observations can be achieved using ground goniometers: the Portable Apparatus for Rapid Acquisitions of Bi-directional Observations of Land and Atmosphere (PARABOLA) [26], the field goniometer system (FIGOS) [27], the portable Multi-Angle Observation System (MAOS) [23,28], and others [29,30]. However, only a small number of studies have been published on LST angular normalization or on the simultaneous retrieval of directional emissivity and temperature from multi-angular TIR images mainly due to a lack of multi-angular observations at satellite level. To date, only ATSR (Advanced Along Track Scanning Radiometer) series satellites [31,32] have provided multi-angular observations of thermal infrared (TIR) data because of the difficulty and the complexity of manufacturing highly controlled multi-angular sensors and also because no studies that discuss requirements in terms of the number of direction specifications of angular observations have been published. To address this situation, we investigate the minimum requirements of viewing observations and directions for achieving angular normalization of LSTs using multi-angular TIR data to provide technique support for the future design of such sensors. This paper is organized as follows: Section 2 will introduce a new model to simulate the directional thermal radiance and discuss the performance of the kernel-driven bi-directional reflectance distribution function (BRDF) model in the TIR domain as a TIR-BRDF model. Section 3 will investigate the local optimum three-angle combinations for the TIR-BRDF using a single-point pattern and discuss the three-angle combinations using a linear array pattern that is always used in airborne and satellite sensors. Section 4 will apply the TIR-BRDF model to an airborne multi-angular TIR image from the WiDAS (Wide-angle infrared Dual-mode line/area Array Scanner) system [22,33] in the WATER (Watershed Allied Telemetry Experimental Research) campaign [34]. Sections 5 and 6 will provide discussions and conclusions of this paper. The full names and the corresponding abbreviations of some terms are listed in the Table A1 of Appendix.

2. Modeling of Directional Thermal Radiation

Modeling the directional thermal radiance (DTR) of homogenous or heterogeneous canopies has attracted a substantial amount of attention. The main reason for the DTR is that the fractions of different components with different temperature and emissivity vary with the changes of viewing angles. Based on this concept, considering N different components in a pixel under clear-sky conditions, a generalized parameterization of the DTR can be expressed as [22]:

L(v , s , ) B[DBT (v , s , )] (v , v )B[Te (v , s , )]

N

fi (v , s , ) i Bi (Ti ) Lmulti

(1)

i1

where θv and θs are the viewing zenith angle (VZA) and the solar zenith angle (SZA), respectively; φ is the relative azimuth angle (RAA) between the viewing azimuth angle (VAA) φv and the solar azimuth angle (SAA) φs; L(θv, θs, φ) is the directional thermal radiance (the inversion of this term according to the Planck’s law B[] will produce the directional brightness temperature (DBT)); Meanwhile, ε(θv,φv) is the pixel ensemble directional emissivity and Te(θv, θs, φ) is the directional

Sensors 2015, 15

7540

effective temperature; fi are the fractions of various components, such as the commonly used sunlit and shadowed leaves and sunlit and shadowed soil; Ti and εi are the temperatures (unit: K) and emissivities of those components, and note that εi are assumed to be independent of viewing angle; and Lmulti is the multi-scattering between soil and leaves and between leaves inside the canopy.

We consider the homogenous canopy and estimate the component fractions fi from a parameterization [35] of the Scattering by Arbitrarily Inclined Leaves (SAIL) Hotspot (SAILH) model [36,37] using the solar and viewing directions, canopy structure and model parameters. In addition, only three components (leaves and sunlit and shadowed soils) are considered because the temperature difference between sunlit and shadowed leaves is very small compared to the temperature difference between sunlit and shadowed soils. Furthermore, we calculate the multi-scattering term Lmulti with the following equation:

Lmulti (1 M )Lleaf (1 g ) (1 )[1 b()M ][1 b()](1 v) Lleaf

(2)

The first part of the right-hand side of Equation (2) is the downward leaves’ radiation reflected by the soil, and the second part is the multi-scattering interior to the canopy; εg and εv are the emissivities of ground soil and leaf, respectively; b(θ) is the directional gap frequency in the θ direction; and M is the hemispheric gap probability of the canopy. For a canopy with a spherical leaf angle distribution and random dispersion, b(θ) and M can be expressed as Equation (3), in which LAI is the leaf area index. α denotes the cavity effect that accounts for multiple scatterings inside the canopy, the details of which can be found in [38].

b() exp[ 0.5

LAI ] , and M 1

2

b()d exp(0.825LAI )

(3)

cos

2

Table 1 lists the main input variables of the simulation DBT, including the solar zenith and azimuth angles; LAI of the canopy; the temperatures of the leaves (Tleaf), sunlit soil (Tsun_soil) and shadowed soil (Tshd_soil); and the emissivities of leaf (εv) and soil (εg). Note that the simulated DBT is assumed to be atmospherically corrected.

Table 1. The main input parameters for the directional brightness temperature (DBT) simulation.

SAA SZA

LAI

Tleaf

Tsun_soil

Tshd_soil

εv

εg

120° 30° 0.5~5.0 with a step of 0.5 305 K 320 K 315 K 0.985

0.95

In addition, we also used the kernel-driven BRDF model, which is always used to link the bi-directional reflectance to the viewing geometry in the visible/near-infrared wavelengths and to link the DTR and DBT to the viewing geometry by replacing the bi-directional reflectivity with the DTR [22,39]. As a result, the kernel-driven BRDF in the TIR field is written as Equation (4). For simplicity, the new kernel-driven BRDF model is henceforth called TIR-BRDF.

DTR(v,s ,) B[DBT(v,s ,)] fiso fvol kvol (v, s , ) fgeo kgeo (v,s ,)

(4)

where fiso is the isotropic scattering term, fvol is the coefficient of the volumetric kernel kvol, and fgeo is the coefficient of the geometric kernel kgeo. The Ross-Thick kernel developed in [40] and the Li-SparseR kernel derived in [41] on the basis of the geometric-optical mutual shadowing BRDF model [42] are used as the volumetric and geometrical scattering kernels, respectively. A suitable

Sensors 2015, 15

7541

expression of kvol was derived by Roujean et al. [40], called the Ross-Thick kernel for its assumption of a dense leaf canopy (see Equation (5)). It is a single-scattering solution of the radiative transfer equation for plane-parallel dense leaf canopy with uniform leaf angle distribution, a Lambertian background and equal leaf transmittance and reflectivity, but it does not account for the hotspot effect.

( / 2 ) cos sin

kvol (v , s , ) cos cos 4

(5)

v

s

where, ξ is the phase angle, related to the sun-target-observer position as:

cos v cos s sin v sin s cos

(6)

The Li-Sparse kernel derived from the geometric-optical mutual shadowing BRDF model [42] was reported to work well with the observed data. The original form of this kernel is not reciprocal in the viewing and solar directions, a property that is expected from homogeneous natural surface viewed at coarser spatial resolution, but then was refined to be reciprocal by assuming the sunlit component in the viewed scene simply varies as 1/cosθs. As a result, the reciprocal model was the Li-SparseR kernel [43,44]:

kgeo (v , s , ) O(v , s , ) sec 'v sec 's 1 (1 cos ')sec 'v sec 's (7) 2

with,

O 1 (t sin t cos t)(sec 'v sec 's )

cos t h b

D2 (tan 'v tan 's sin )2 sec 'v sec 's

D tan2 'v tan2 's 2 tan 'v tan 's cos

cos ' cos 'v cos 's sin 'v sin 's cos

' tan1(b tan ) , ' tan1(b tan ) ,

v

r

v

s

r

s

where, O is the overlap area between the view and solar shadows. The term cost should be constrained to the range (−1, 1), as values outsides of this range imply no overlap and should be disregarded. The ratio h/b and b/r are the dimensionless crown relative height and shape parameters, respectively, and should be preselected. In this thesis, h/b = 2 and b/r = 1, i.e., the spherical crowns are separated from ground by half their diameter.

Figure 1 shows the simulated DBT distribution (the first column), the fitted DBT from the TIR-BRDF model (the second column), and their temperature difference (the third column) at (a) LAI = 0.5; (b) LAI = 1.0 and (c) LAI = 2.0 based on the input variables shown in Table 1 for the DBT simulation using Equation (1) and the SAILH parameterization model [35]. The maximum zenith angle is constrained to 60°because the kernel-driven BRDF model was reported to obtain unacceptable results for larger zenith angles [38,45].

Sensors 2015, 15

7542

(a) LAI = 0.5

(b) LAI = 1.0

(c) LAI = 2.0 Figure 1. The hemispheric distribution of the simulated DBT from Equation (1) (the first column), the fitted DBT from the TIR-BRDF model (the second column) and their temperature difference (the third column) at (a) LAI = 0.5; (b) LAI = 1.0 and (c) LAI = 2.0, respectively. As observed from Figure 1, the angular variation in the DBT can be as high as approximately 3.0 K under the simulation conditions (Table 1). However, this angular variation is based on the components’ temperature difference and canopy structure: larger component temperature differences generally lead to larger angular variations. The TIR-BRDF model produces a good fit of the DBT distribution, with errors lower than 0.3 K. The hotspot effect in the solar direction is significant; however, we also find that the maximum temperature difference occurs around this direction and that the DBT of the TIR-BRDF model is generally larger than that of the simulation result in the directions around the solar beam. The increase in LAI decreases the DBT because the fraction of leaves increases and the leaves’ temperature is lower than that of the soils. Furthermore, the LAI also influences the temperature difference caused by the regression used in the TIR-BRDF model. Figure 2 shows the RMSE (root-mean-square error) and the maximum temperature difference as a function of LAI. The figure

Sensors 2015, 15

7543

indicates that the RMSE is smaller than 0.1 K and that the maximum temperature is smaller than 0.3 K; both obtain their global maximums at approximately LAI = 2.

Figure 2. The influence of LAI on RMSE and maximum temperature difference from the TIR-BRDF model.

Furthermore, we also implemented the TIR-BRDF model using different solar positions, component temperatures and emissivities, and canopy structures, and the results further confirmed the performance of the TIR-BRDF model for calculating the DBT distribution. Note that because the component emissivities are considered to be angle independent and because, consequently, the DBT of a bare surface is isotropic, we only discuss at the level of the canopy rather than for the case of a bare surface.

3. Viewing Angle Specification for Angular Normalization of LST

The above result is obtained using the TIR-BRDF model with all viewing zenith and azimuth angles. However, in practice, it is almost impossible to observe the same target at numerous directions. On the contrary, only some angular observations are usually conducted. Consequently, a problem arises: how many observations and under what viewing angles are needed for the TIR-BRDF model to accurately fit the DBT and consequently enable the angular normalization of the LST? To determine the optimum observation groups, this paper used two different patterns in the following discussions: single-point patterns and linear-array patterns.

3.1. Single-Point Pattern Analysis for Multiple-Viewing Angle Specification

According to Equation (4), at least three angular observations are needed to drive the TIR-BRDF model. Additional observations can improve the reliability of the fitted coefficients in theory, but the errors in those observations may produce greater uncertainty in the fitted result. Therefore, we start with only three viewing angles to determine the local optimum viewing angle combination that produces the smallest error in the TIR-BRDF model and to provide various suggestions for the future development of multi-angular airborne or space-borne TIR sensors. Therefore, to better facilitate the

Sensors 2015, 15

7544

development of a mechanical design of such a sensor, we herein assume that all observations are performed with the azimuths in the same plane, i.e., their azimuths equal φ or φ + 180°.

Taking φ = 0° for example, the three observations’ zenith angles vary in the azimuth plane (0°~180°) at a step of 10°to a maximum of 60°in the zenith direction. Different observations are required to have different value of zenith or azimuth angles. As a result, there are 126 groups with three angles. To determine the optimum group, we first use Equation (1) to simulate hemispheric DBTs under varying solar positions (SAA from 0°to 330°with a step of 30°, SZA = 10°, 30°and 50°), LAIs (from 0.5 to 5.0 with a step of 0.5) and the component temperatures and emissivities in Table 1, thus resulting in a total of 396 cases. Consequently, the three coefficients in Equation (4) are obtained by regressing the DBT with the known kernels (kvol and kgeo) calculated in different viewing geometry, and this process is additionally controlled by an optimization algorithm on the basis of that the brightness temperature at the hotspot direction (i.e., solar beam direction) has the large value because this direction corresponds to the largest proportions of sunlit soil and leaves which have a higher temperature. Table 2 shows the frequency of the temperature RMSEs in the ranges of [0.0, 0.5] K and [0.5, 1.0] K and the number of times the maximum and minimum RMSEs are obtained for different angle groups. Furthermore, because it was found that the viewing angle group with nadir viewing direction (VZA = 0° and VAA = 0°) had better results than those without such nadir observation, and also because of space limitations of this paper, we did not display those results without the nadir observation.

Furthermore, four criteria are proposed here to filter the final optimum angle combinations:

(1) Most of the RMSEs in the TIR-BRDF model should fall in the range of 0.0~0.5 K. (2) The maximum difference should not or seldom occur for acceptable combinations of angles. (3) Most of the minimum differences can be obtained using the combinations. (4) Angle difference of the adjacent viewing angles should be large enough in order to cause large

difference in the component fractions and then lead a large angular variation in the observed brightness temperature. But a very large VZA are not recommended in order to avoid large pixel size differences in the angle combination.

Table 2. The frequency of the root-mean-square error (RMSE) in different three-direction combinations. Columns 1 and 2 are the VAA and VZA for the first direction observation, and columns 3 and 4 are for the second direction observation, and columns 5 and 6 are for the third direction observation.

1st VAA

0 0 0 0 0 0 0

1st VZA

0 0 0 0 0 0 0

2nd VAA

0 0 0 0 0 0 0

2nd VZA

10 10 10 10 10 10 20

3rd VAA 180 180 180 180 180 180 180

3rd VZA

10 20 30 40 50 60 10

RMSE [0.0~0.5] K

195 241 299 316 337 323 243

RMSE [0.5~1.0] K

44 43 40 33 7 4 42

Max *

122 49 16 15 1 5 45

Min #

0 1 2 0 12 11 1

Sensors 2015, 15

7545

Table 2. Cont.

1st 1st 2nd 2nd 3rd 3rd VAA VZA VAA VZA VAA VZA

RMSE (0.0~0.5) K

RMSE (0.5~1.0) K

Max *

Min #

0

0

0 20 180 20

295

18

11

0

0

0

0 20 180 30

315

34

14

4

0

0

0 20 180 40

337

20

2

3

0

0

0 20 180 50

359

3

1

15

0

0

0 20 180 60

354

7

19

11

0

0

0 30 180 10

299

39

17

2

0

0

0 40 180 10

316

35

16

0

0

0

0 40 180 20

337

20

2

3

0

0

0 40 180 30

361

2

0

6

0

0

0 40 180 40

364

1

0

3

0

0

0 40 180 50

373

0

0

26

0

0

0 40 180 60

373

0

1

47

0

0

0 50 180 10

336

8

1

14

0

0

0 50 180 20

359

4

1

16

0

0

0 50 180 30

360

0

0

14

0

0

0 50 180 40

373

0

0

26

0

0

0 50 180 50

353

2

2

14

0

0

0 50 180 60

365

8

5

69

0

0

0 60 180 10

322

5

4

10

0

0

0 60 180 20

354

8

22

12

0

0

0 60 180 30

373

0

1

10

0

0

0 60 180 40

373

0

1

47

0

0

0 60 180 50

365

8

3

69

0

0

0 60 180 60

349

0

21

0

* Max stands for the frequency of the maximum temperature difference, while # Min stands for the frequency

of the minimum temperature difference.

According to the above four criteria and the results shown in Table 2, we find that large VZAs with respect to the nadir direction can reduce the fitting error, and several angle combinations of VAA and VZA can be considered as candidates for the optimum groups: ① [(0°, 0°), (0°, 30°), (180°, 50°)]; ② [(0°, 0°), (0°, 40°), (180°, 50°)]; and ③ [(0°, 0°), (0°, 40°), (180°, 60°)] (see highlight in Table 2). Although the combination [(0°, 0°), (0°, 50°), (180°, 60°)] produces the minimum error most often, the number of times the maximum error is obtained (five times) and the large viewing zenith angles make this group impractical.

Among the three candidates, group ③ cannot be recommended because of its large VZA in the direction (180°, 60°), whose pixel sizes are up to four times the pixel sizes at nadir observation if the sensor’s IFOV (instantaneous field of view) remains constant in every direction. Group ② performs slightly better than group ①. However, the final decision cannot be made without a sensitivity analysis of the errors in the observed data. To investigate the sensitivity of the different angle groups to the DBT errors, artificial temperature noises within [−0.5, 0.5] K and [−1.0, 1.0] K with a uniform distribution are added to the simulated DBT of Equation (1). The three coefficients of the TIR-BRDF model are recalculated and subsequently used to obtain the hemispheric DBT. The output of the

Sensors 2015, 15

7546

TIR-BRDF model is constrained whereby the DBT must reach its maximum value in the solar direction because the sunlit components in this direction have the largest percentage in this hotspot direction.

Figure 3 shows the histograms of the temperature RMSEs for the three-angle TIR-BRDF model using the angle groups ①–③ with respect to simulated random measurement error in [−0.5, 0.5] K and [−1.0, 1.0] K included in the input DBT data. The figure illustrates that group ③ produces the smallest errors under the two noise conditions, followed by groups ② and ①. The RMSE percentages of groups ①–③ in the range 0.0~1.0 K are approximately 92.4%, 96.6% and 96.3% for the noise level [−0.5, 0.5] K and 80.1%, 88.5% and 91.1% for the noise level [−1.0, 1.0] K. These results indicate that the three groups can produce a temperature error of less than 1.0 K for most cases if the noise involved in the observed DBT data is no more than 1.0 K. However, it is obvious that groups ② and ③ are less sensitive to the noise in the observed data. As stated above, group ③ cannot be used as the optimum viewing angle combination because of the large viewing zenith angle and the occurrence of maximum errors. As a result, group ② should be the optimum group in theory. However, because the pixel size increases with increasing VZA, the ratios of pixel area at VZA = 30°, 40°and 50°to that of nadir are 1.32, 1.70 and 2.40, respectively. It is necessary to use a pixel with a middle area to connect the observations from nadir to large VZAs, especially for heterogeneous surfaces. In this case, the VZAs in group ① produce a more continuous series of pixel sizes (1, 1.3 and 2.4 times the nadir pixel size). In addition, the angle intervals in the slant direction of group ① (i.e., 20°) are larger than those of group ② (i.e., 10°), which makes group ① less sensitive to the errors in the observation angle in theory. Therefore, according to the single-pattern analysis, we prefer the angle combination of group ① and take this group as the local optimum angle combination for the DBT regression using the three-angle TIR-BRDF model.

Figure 3. Histograms of the temperature RMSE caused by the three-angle TIR-BRDF model with respect to temperature noise about (a) [−0.5, 0.5] K and (b) [−1.0, 1.0] K included in the input DBT data, respectively.