EXAMINATION OF THE TERRAIN EFFECT FOR TERRESTRIAL ALBEDO ESTIMATION VIA BRDF MODEL PARAMETERS

: In this paper, we examine the effect of terrain on terrestrial albedo estimation. Terrestrial albedo is one of the most important parameters for understanding the global heat balance. The existing approach for estimating terrestrial albedo involves the estimation of model parameters of the bidirectional reflectance distribution function (BRDF) based on measurements obtained at different geometries. Then, narrowband albedos are estimated from the BRDF model parameters and the broadband albedo is finally estimated by narrowband-to-broadband conversion. Previous studies have not considered the terrain effect for generating the terrestrial albedo. Experiments using in situ measurements showed that the BRDF model, which transforms the geocoordinate of the reflectance of the shadowed terrain, generates the best accuracy. The improvement in the accuracy by the terrain effect correction is limited, and therefore further investigations using more in situ and simulated data are necessary for operational products.


INTRODUCTION
Terrestrial albedo is one of the most important physical parameters for understanding the global circulation of water and heat.Land cover has anisotropic reflectance characteristics that are defined by the bidirectional reflectance distribution function (BRDF), the angular integration of which yields an albedo value.Moderate resolution imaging spectroradiometry (MODIS) terrestrial-albedo products are generated using this BRDF-driven approach (Strahler et al., 1999).Among many proposed BRDF models, kernel-driven models are regarded as robust semiempirical BRDF models that are applicable to any land cover type (Lucht et al., 2000).The kernel function is determined by viewing and illumination geometries.Currently, a common estimation approach for multispectral satellite observations involves a series of processing steps, typically including atmospheric correction for estimating land-surface directional reflectance, angular modeling for calculating spectral albedos, and narrowband-to-broadband albedo conversion (Liang, 2000;Liang, 2004).Such albedo products are acceptably accurate when applied to paddy fields (Susaki et al., 2007).
A BRDF-driven albedo estimation requires a certain number of observations.For example, MODIS operational BRDF/albedo products use at least seven good-quality reflectance observations within a 16-day window.This requirement is challenging for tropical and subtropical climate regions, for which optical images are very often contaminated by clouds.Simulations have revealed that estimates of BRDF model parameters become less stable as more observations are contaminated by noise (Susaki et al., 2004).Consequently, data are often unavailable for MODIS BRDF/albedo products in these regions.Therefore, in this paper, we examine two approaches for stably estimating the terrestrial albedo in areas for which limited observations are available.Narrowband albedos can be estimated using a single observation; however, a more sophisticated approach is to use several observations and select the minimum reflectance.The latter is a BRDF approach, substituting the reflectance for the spectral albedo.Section 2 explains an algorithm to estimate albedo in this way.Section 3 describes the study area and the terrestrial-albedo and bidirectional reflectance factor (BRF) data used in this research.The experimental results are reported in Section 4 and discussed in Section 5. Finally, Section 6 concludes this paper.

ALBEDO ESTIMATION VIA THE BRDF MODEL
Albedo can be defined as the ratio of the outward flux to the inward flux through a hemisphere over a point of interest.The wavelength range depends on the type of albedo.For example, shortwave albedos usually target the range of 0.3-3.0μm. Figure 1 shows a flowchart of an algorithm to generate broadband albedo from spaceborne multi-angular data.Top-of-atmosphere radiance is converted into surface spectral reflectance via atmospheric correction, such as by secondary simulation of the satellite signal in the solar spectrum (6S) (Vermote et al., 1997).Spectral albedo is generated by referring to the BRDF angular model.Finally, broadband albedo is estimated as a linear combination of multiple spectral albedos.This conversion is referred to as narrowband-to-broadband conversion.Hereinafter, we describe the terminology and ideas geared toward estimating broadband albedo.

BRDF Model
A kernel used for estimating BRDF and terrestrial albedo is a function of the bidirectional reflectance determined by the viewing and illumination geometries.In general, we consider two types of optical scattering from an object on the terrestrial surface: volumetric and geometric.
For volumetric scattering, Ross (1981) developed a kernel approximation of the directional reflectance above a horizontally homogeneous plant canopy, and Roujean et al. (1992) derived the Ross-thick kernel as follows: Here, fiso, fvol, and fgeo are unknown coefficients with which users can select any combination of volume-and geometric-scatterings kernels.In this research, we select a popular combination, namely, the Ross-thick kernel for volume scattering and the Li-sparse kernel for geometric scattering.The kernel values are determined once the illumination and viewing geometries are specified.The three unknown coefficients are determined by minimizing the least-squares error function Here, n denotes the number of samples and r denotes the observed reflectance.In operational applications, MODIS BRDF products are generated as follows.If at least seven cloud-free surface observations are available during a 16-day period, a full model inversion is attempted.The available data are first evaluated to discard outliers, and additional checks are performed to guarantee that the kernel weights are positive.If the data pass these evaluations, a full (or normal) inversion is performed to establish the BRDF parameter weights that provide the best rootmean-square error (RMSE) fit (Strahler et al., 1999).

Albedo Estimation
The estimated BRDF model parameters then allow black-sky and white-sky albedos to be determined.The black-sky albedo is a virtual albedo in the absence of a diffuse component.The whitesky albedo is also a virtual albedo in the absence of a direct component when the diffuse component is isotropic.The blacksky and white-sky albedos are expressed by ( 11) and ( 12), respectively: The coefficients in ( 13) and ( 14) are given in Table I.As shown in (15), the actual albedo a (λ) at wavelength λ is expressed as a linear combination of black-sky and white-sky albedos using the atmospheric optical depth τ (λ):  where S denotes the fraction of diffuse skylight.
Broadband albedo is estimated using the linear regression of several narrowband albedos, called narrowband-to-broadband (NTB) conversion (Liang, 2000).Liang et al. (2002) introduced conversion formulas for various optical sensors.The broadband albedos produced by these formulas have an uncertainty of approximately 0.02.For example, the formula for estimating shortwave albedo from MODIS data is expressed by Equation ( 16): Here, αi denotes narrowband albedo derived from the reflectance of band i.

Reflectance of Geometric Scattering
Reflectance of geometric scattering is expressed by Equation ( 17) (Wanner and Strahler, 1995): Here,  ,  ,  , and  are the proportions of sunlit crown, sunlit ground, shadowed crown, and shadowed ground, respectively. ,  ,  , and  are the reflectances of sunlit crown, sunlit ground, shadowed crown, and shadowed ground, respectively.

Figure 2. Modeling of shadow effects on geometric scattering
Assume that discrete objects are available in three-dimensional coordinates.The gap probability (, ) can be expressed by Equation ( 18).

𝑞(𝜃, 𝜙) = 𝑒𝑥𝑝 −𝜆 • 𝐴(𝜃, 𝜙) . (18)
Here,  and  are the viewing angle and azimuth angle, respectively. is the density of object centers on a plane at the base of the layer.(, ) is the average area of an object projected at angle  and . can be expressed by the number of objects n and (, ): As a result, the proportional area of the crown from a viewpoint is expressed as follows: Here,  and  are the viewing zenith and azimuth angles, respectively.Similarly, the proportion of non-crown from the viewpoint is expressed as follows: The proportion of both sunlit and viewed ground is calculated by multiplying the proportions of ( ,  ) and ( ,  ) and excluding the proportion of overlapping area: Here,  and  are the illumination zenith and azimuth angles, respectively.Oarea is the mean area of overlap between the projected areas of objects onto the ground at angle ( ,  ) with its projected area at angle ( ,  ).
Now, we consider the proportion of sunlit crown relative to sunlit and shadowed crown.This relative proportion of sunlit crown can be expressed as follows: Here, g′ is the scattering phase angle between the direction vectors for the illumination and viewing directions.Then, pc can be rewritten using Equations ( 20) and ( 23): In the Li-sparse model, it is assumed that Rc is equal to Rg, and that RT and Rz are equal to 0. Using Equations ( 17), ( 22), and (24), Equation ( 25) is derived: Because  is a small value, () 1 +  .Equation ( 25) is rewritten as Equation ( 26): − ( ,  ) + ( ,  ) −  . (26) (, ) can be approximated as an area of a sphere with radius r: The geometric kernel Kgeo can be defined from Equation ( 27): Here, O of Equation ( 28) is equivalent to that of Equation (4).c1 is equivalent to fgeo in Equation ( 9), and c2 is a component of fiso in Equation ( 9).While Equation ( 29) is independent of the solar zenith angle, Equation ( 3) is dependent on this angle (Wanner et al., 1995).In this research, we use Equation (3) as the geometric kernel Kgeo.
The traditional kernel-driven BRDF model used for estimating terrestrial albedo assumes that the ground is flat; a model considering the effect of terrain is necessary for the accurate estimation of the BRDF model and terrestrial albedo.Figure 3 shows that the area of crown projected to the ground depends on the slope.Accordingly, the proportion of the total area and the observed reflectance from the scene may change.In this research, we investigate the change in area by applying the coordinate transform or area conversion.

Coordinate Transformation for Terrain Effect Correction
To estimate BRDF for various forest structures, Nagamine and Murakami (1981) proposed a formula to transform the coordinate system so that the geolocation of the sun and satellite relative to the terrain surface is calculated.Figure 4 shows the concept of the relation between the coordinates.Equation ( 31) represents the generation of the unit vector of the transformed coordinate system from the slope elevation and aspect.
Here, s and t are the slope elevation and azimuth angles, respectively.( , ,  ) is the original coordinate system, and (, , ) is the transformed coordinate system.
Next, we calculate the geolocation of the sun and satellite in the transformed coordinate system.In the original coordinate system, a unit vector n, expressed by Equation (32), is calculated.Here, ZE and AZ are the zenith and azimuth angles of an object in the original coordinate system, respectively.The vector of Equation ( 32) is calculated for the sun and satellite.This vector n is substituted into (, , ) of Equation ( 31), and the unit vector in the transformed coordinate system, n′, is generated.By (33) We can correct for the terrain effect by substituting the transformed zenith and azimuth angles of the sun and sensor.
Previous studies have applied a Kgeo correction.However, in this research, we investigate two models to correct for the terrain effect by transforming the whole Kgeo and transforming the term related to  , which is impacted most substantially by the terrain effect and is important for clarifying the performance of the model.

Area Conversion for Terrain Effect Correction
Wu and Strahler (1994) proposed a formula expressed by Equation ( 33) that converts the area of the shadowed background  into the adjusted area  to correct for the terrain effect: Here, β is the slope elevation angle, θ is the solar zenith angle, and  is the relative azimuth angle between the normal slope and an object, sun, or satellite.c is a coefficient.In Equation ( 29),  has terms representing the shadowed areas from illumination direction sec  and viewing directions sec  .We corrected for the terrain effect in kernel-driven BRDF models by applying c of Equation ( 33) to those terms.
Previous research has reported a terrain effect correction in BRDF models, but the effect on the albedo estimated via BRDF model parameters has not been reported.Therefore, in this research, we assessed the accuracy of the albedo estimated using BRDF model parameters.

DATA
We selected as a study area a tower at the Yamashiro Experimental Forest in southern Kyoto Prefecture, Japan (Yamashiro experimental forest, 2018).The soils are Regosols with sandy loam or loamy sand texture and contain fine gravel composed of residual quartz crystals from granite parent material (Kaneko et al., 2007).Deciduous broad-leaved, evergreen broadleaved, and coniferous tree species account for 66%, 28%, and 6% of the living tree biomass, respectively (Goto et al., 2003).
Table 2 shows the details of the in situ measurements.We measured BRF for wavelengths of 300-2500 nm at the top of the tower at 26.5 m using a FieldSpec 3 spectroradiometer (ASD, 2018).For BRF measurements, we set the azimuth angle between the solar and viewing directions to 15°, 45°, 90°, 135°, and 180° to minimize measurement time.We set the angle to 15° instead of 0° to avoid contamination by the equipment shadow.The viewing angle was set to 0°, 15°, 30°, and 45°, as shown in Figure 5.We used a 5° field-of-view lens for the measurement.We also measured the broadband albedo for wavelengths of 300-3,000 nm using an MR-60 spectrometer (EKO, 2018).The leaf area index was also measured using an LAI-2200C plant canopy analyzer (LI-COR, 2018).
For the terrain data, we downloaded a digital elevation model with a spatial grid of 10 m from the Geospatial Information Authority of Japan.The calculated slope elevation and aspect were 6.8° and 228.6°, respectively.

EXPREMENTS
We derived the reflectance of MODIS bands 1 to 7 from the measured BRF data.The measured digital numbers of the samples were converted into reflectances.Calculating the spectral albedo using Equation ( 15) requires the atmospheric optical depth τ (λ) at wavelength λ.In the experiments, we used τ (λ = 550 nm) for all bands to simplify the calculation.
In the experiments, we examined five approaches to geometric kernels in the kernel-driven BRDF model.In Model (a), we applied a coordinate transformation to the surface reflectance term in  and in Model (b), we applied it to the whole  .In Model (c), we applied an area conversion to the surface reflectance term in  and in Model (d), we applied it to the whole  .For reference, in Model (e), we calculated the albedo by applying the existing method reported by Strahler et al. (1999), which does not consider the terrain effect.
Table 3 shows the albedos estimated using Equation ( 16) for NTB conversion.Some albedos in Table 3 had large residuals.Therefore, we derived another NTB conversion formula from in situ measurements, shown in Equation ( 34): The non-negative coefficients in Equation ( 34) were derived from 6 measurements in bare soil, 1 in grassland, and 3 in the forest.Seven bands of narrowband albedos were used as explanatory variables and broadband albedos were used as explained variables.Equation ( 34) was obtained when Akaike's Information Criterion (Akaike, 1973) was minimized.The validation of Equation (34) using another data set shows that the RMSE was 0.014, and the formula was effectively accurate.Table 4 shows the albedos estimated using Equation (34) for NTB conversion.

DISCUSSION
In the experiments, we examined five approaches to geometric kernels in the kernel-driven BRDF model using in situ measurements in a forest.As a result, Model (a), which transforms the coordinate of the surface reflectance term in  , showed the best accuracy for broadband albedos.As shown in both Tables 3 and 4, Model (a) generated the best results for the data observed on April 26, May 16, and July 10, 2018.For the data observed on March 29, 2018, the results for Model (a) were the second best, but the differences from the best models in Tables 3 and 4 were 0.0002 and 0.0012, respectively.For the data observed on September 25, 2017, the estimates were not close to the actual values.Now, we discuss the factors that contributed to the results.We separate the reflectance of the terrain and that of the canopy.The reflectance of the terrain is deeply dependent on the terrain effect; the relative geometry between the illumination and the slope is affected by the terrain elevation and azimuth angles.On the other hand, the canopy reflectance term may be independent of the terrain effect because the canopy arrangement is not determined by only the terrain elevation and azimuth angles.In this regard, the Model (a), which transforms the coordinate of the surface reflectance term in  , generates more reasonable results than those of Model (b), which transforms the coordinate of the whole  .Model (b), which also transforms the reflectance of the canopy, is expected to overcorrect for the terrain effect.
The area conversion may have a similar performance to that of the coordinate transformation."Abs dif2," defined as the difference between the estimate of interest and the value estimated by the existing model (Tables 3 and 4) indicates that the results obtained by Model (c) are similar to those by Model (a), and the results obtained by Model (d) are similar to those by Model (b).Note that the area conversion requires the canopy shape parameters.The area conversion in Models (c) and (d) may work under an appropriate model of the canopy.However, the results may depend on the canopy shape model.It may be difficult to model the shape of the canopy.Therefore, the coordinate transformation, independent of the canopy shape model, is more reasonable than the area conversion.
In addition, we compared the albedos estimated using the NTB conversion formula used for MODIS products with the albedos estimated using the calibrated NTB conversion formula, expressed by Equation ( 34).The results are shown in Tables 3  and 4, respectively.It is not the case that the results in Table 4 exceed those in Table 3.The performance of the calibration of the NTB conversion formula should be investigated using more in situ data.

CONCLUSIONS
We investigated the effects of the slope inclination angle on the terrestrial albedo via a kernel-driven BRDF model.The terrain effects can be expressed in a geometric scattering kernel and are derived by transforming the local coordinate or by changing the shadow area.We compared five models using in situ BRF data measured at the top of a tower in a forest.As a result, Model (a), which transforms the coordinate of the surface reflectance term in  , showed the best accuracy for broadband albedos.
As a future task, we will investigate the performance of the models using more in situ BRF data obtained under different LAI values.In addition, we will use simulated data obtained using our simulator (Jin and Susaki, 2017).Finally, we evaluate the increase in computation time by considering the terrain effect and the improvement in accuracy and determine which is more suitable for operational terrestrial albedo products, BRDF models without or with terrain effects.
Table 3. Broadband albedo obtained by applying the NTB conversion formula reported by Liang (2000).Gray column denotes the best result."Abs dif1" and "abs dif2" denote the differences from the actual value and from the value estimated using the existing model, respectively.

Figure 1 .
Figure 1.Algorithm for estimating broadband albedo using the kernel-driven BRDF model from spaceborne sensor data

Figure 3 .
Figure 3. Change in the crown area projected to the ground caused by the terrain effect.θ is the solar zenith angle, b is half the height of the crown, and r is the crown radius.

Figure 4 .
Figure 4. Transformation of the local coordinate.s and t are the slope elevation and azimuth angles, respectively.
substituting n′ into Equation (33), the transformed zenith and azimuth angles of the object, ZE′ and AZ′, are generated, respectively:

Figure 5 .
Figure 5. Measurement of BRF at the top of a tower in a forest (a) from an oblique viewpoint and (b) at the nadir viewpoint.

Table 2 .
Measurements in the study area

Table 4 .
Broadband albedo obtained by applying the NTB conversion formula generated by calibration with in situ measurements.Gray column denotes the best result."Abs dif1" and "abs dif2" denote the differences from the actual value and from the value estimated using the existing model, respectively.