Ziti Jiao
email: jiaozt@bnu.edu.cn

Crystal B Schaaf

Yadong Dong

Miguel Román

Michael J Hill

Jing M Chen

Zhuosen Wang

Hu Zhang

Edward Saenz

Rajesh Poudyal

Charles Gatebe

Francois-Marie Bréon

Xiaowen Li

Alan Strahler

A method for improving hotspot directional signatures in BRDF models used for MODIS

Keywords: BRDF, CAR, MODIS, POLDER, multiangle remote sensing, hotspot signature, hotspot kernel, linear RTLSR model, airborne measurements

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Semi-empirical kernel-driven linear Bidirectional Reflectance Distribution Function (BRDF) models have been widely used to determine the properties of complex heterogeneous environments from multi-angle remote sensing. These models have been used to produce the routine BRDF/Albedo products from the Moderate Resolution Imaging Spectroradiometer (MODIS) [START_REF] Lucht | An algorithm for the retrieval of albedo from space using semiempirical BRDF models[END_REF][START_REF] Schaaf | First operational BRDF, albedo nadir reflectance products from MODIS[END_REF]Schaaf et al., 2011), the Polarization and Directionality of the Earth's Reflectances (POLDER) [START_REF] Bicheron | Bidirectional reflectance distribution function signatures of major biomes observed from space[END_REF][START_REF] Bacour | Variability of land Surface BRDFs[END_REF], the Meteosat Second Generation (MSG) (van Leeuwen and Roujean, 2002;[START_REF] Geiger | Land Surface Analysis Satellite Application Facility (LSA-SAF) Product User Manual (PUM) on Albedo[END_REF], and the Visible/Infrared Imager/Radiometer Suite (VIIRS) on board the platforms of the Suomi National Polar-orbiting Partnership (NPP) [START_REF] Justice | Land and Cryosphere Products from Suomi NPP VIIRS: Overview and Status[END_REF]. They have been also used to retrieve canopy structure parameters (e.g., [START_REF] Chopping | Large area mapping of southwestern forest crown cover, canopy height, and biomass using the NASA Multiangle Imaging Spectro-Radiometer[END_REF][START_REF] Wang | Retrieval of canopy height using moderate-resolution imaging spectroradiometer (MODIS) data[END_REF][START_REF] Hill | Characterizing vegetation cover in global savannas with an annual foliage clumping index derived from the MODIS BRDF product[END_REF][START_REF] He | Global clumping index map derived from the MODIS BRDF product[END_REF], to examine the improved accuracy of land cover classification [START_REF] De Colstoun | Improving global scale land cover classifications with multi-directional POLDER data and a decision tree classifier[END_REF][START_REF] Jiao | Improving MODIS land cover classification by combining MODIS spectral and angular signatures in a Canadian boreal forest[END_REF]Jiao & Li, 2012), to accumulate and apply prior knowledge of BRDF archetypal shapes [START_REF] Li | A priori knowledge accumulation and its application to linear BRDF model inversion[END_REF][START_REF] Jiao | An Anisotropic Flat Index (AFX) to derive BRDF Archetypes from MODIS[END_REF][START_REF] Jiao | An algorithm for retrieval of surface albedo from small view-angle airborne observations through the use of BRDF archetypes as prior knowledge[END_REF], to couple surface reflectance with atmospheric scattering for improving atmospheric correction algorithms (Hu et al., 1999;[START_REF] Wang | Assessment of biases in MODIS surface reflectance due to Lambertian approximation[END_REF][START_REF] Román | Assessing the coupling between surface albedo derived from MODIS and the fraction of diffuse skylight over spatially-characterized landscapes[END_REF][START_REF] Litvinov | Models for surface reflection of radiance and polarized radiance: Comparison with airborne multi-angle photopolarimetric measurements and implications for modeling top-of-atmosphere measurements[END_REF], to correct for the effect of remotely-sensed anisotropic reflectance (e.g., [START_REF] Leroy | Sun and view angle corrections on reflectances derived from NOAA/AVHRR data[END_REF][START_REF] Li | The bidirectional effects of AVHRR 1322[END_REF], and for monitoring ecosystem disturbance and vegetation dynamics (e.g., [START_REF] Friedl | Global land cover mapping from MODIS: algorithms and early results[END_REF][START_REF] Friedl | MODIS Collection 5 global land cover: Algorithm refinements and characterization of new datasets[END_REF][START_REF] Zhang | Monitoring vegetation phenology using MODIS[END_REF].

Despite their wide use, recent studies have recognized that the directional signatures near the Sun illumination direction (aka the hotspot effect) are often underestimated by the semi-empirical BRDF models (e.g., [START_REF] Chen | A hotspot function in a simple bidirectional reflectance model for satellite applications[END_REF][START_REF] Maignan | Bidirectional reflectance of Earth targets: evaluation of analytical mode ls using a large set of sp aceborne measurements with emphasis with the hot spot[END_REF][START_REF] Román | Variability in surface BRDF at different spatial scales (30 m -500 m) ov er a mixed agricultural landscape as retrieved from airborne and satellite spectral measurements[END_REF] such as RTLSR [START_REF] Wanner | On the derivation of kernels for kernel-driven models of bidirectional reflectance[END_REF][START_REF] Lucht | An algorithm for the retrieval of albedo from space using semiempirical BRDF models[END_REF]. In particular, the volumetric scattering (aka Ross) component of RTLSR, originally derived from a horizontally homogeneous plant canopy [START_REF] Ross | The Radiation Regime and Architecture of Plant Stands[END_REF], does not include all possible correlations between the illumination and observation geometries [START_REF] Kuusk | The hotspot effect in plant canopy reflectance, in Photon-Vegetation Interactions[END_REF][START_REF] Jupp | A hot spot model for leaf canopies[END_REF]Qin & Geol, 1995). Although the geometric-optical (Li-Strahler) component derived from geometric optical (GO) models characterizes a significant hotspot effect [START_REF] Li | Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy: effect of crown shape and mutual shadowing[END_REF], the RTLSR model that linearly combines these two (Ross and Li-Strahler) components has difficulties to simulate the both the magnitude and signature of the hotspot effect. This model deficiency, while not significantly impacting an albedo retrieval [START_REF] Huang | Analysis of BRDF and Albedo Retrieved by Kernel-Driven Models Using Field Measurements[END_REF] (based on the integral of the entire view-illumination hemisphere), may still constrain the application of BRDF models in retrieving canopy structure parameters (e.g., clumping index) that need the hotspot amplitude as primary input (e.g., [START_REF] He | Global clumping index map derived from the MODIS BRDF product[END_REF]Zhu et al., 2012).

Various efforts have been made to improve the hotspot effect for such models. [START_REF] Chen | A hotspot function in a simple bidirectional reflectance model for satellite applications[END_REF] enhanced the hotspot effect of the kernel-driven Roujean BRDF model by multiplying the model with an exponential approximation of a physically-based hotspot function. [START_REF] He | Global clumping index map derived from the MODIS BRDF product[END_REF] and Zhu et al. (2012) suggested that the [START_REF] Chen | A hotspot function in a simple bidirectional reflectance model for satellite applications[END_REF] model might still slightly underestimate the reflectance magnitude when it is used to extrapolate the hotspot.

Recent efforts have focused on correcting the RossThick kernel with a hotspot factor [START_REF] Maignan | Bidirectional reflectance of Earth targets: evaluation of analytical mode ls using a large set of sp aceborne measurements with emphasis with the hot spot[END_REF]; hereafter referred to as the Maignan method) based on the geometrical principles of the intersection of viewed and sunlit leaf areas [START_REF] Jupp | A hot spot model for leaf canopies[END_REF]. With a view to improving retrieval of clumping index (CI) from the MODIS BRDF product, [START_REF] He | Global clumping index map derived from the MODIS BRDF product[END_REF] also developed a correction for the MODIS hotspot amplitude by adding the difference between POLDER and MODIS hotspot BRFs, which has also been used to correct the hotspot magnitude of MISR BRFs for CI retrievals (Pisek, Ryu, Sprintsin & He et al., 2013). On the other hand, Zhu et al. (2012) corrected the hotspot amplitude by multiplying both the geometric optical scattering and volumetric scattering items by the exponential approximation of a modified hotspot function; also based on the MODIS RTLSR model.

In this study, we propose a method for improving the hotspot effect of the linear RTLSR BRDF model. The method revises the RossThick kernel using the corrected exponential approximation of the hotspot function (Chen and Cihlar, 1997, thereafter named RTCLSR model here). The principle of the formation of the hotspot is based on a canopy gap size distribution function, but is approximated by using an exponential function with two free parameters (C 1 and C 2 ) characterizing the height and width of the hotspot effect. The hotspot kernel within-crown and between-crown has very similar shapes, and thus can be directly applied to the scenario where a canopy cover is provided with a uniform leaf orientation of horizontally homogeneous plant canopies, as was used in the assumption in deriving this kernel from canopy radiative transfer theory [START_REF] Ross | The Radiation Regime and Architecture of Plant Stands[END_REF]. Such a correction to the RossThick kernel mainly accounts for the correlation between two gap probabilities from sun and view in the hotspot direction, which was not properly considered when deriving the original RossThick kernel.

To validate this model, the study uses a variety of hotspot data to determine appropriate hotspot-related parameters including POLDER, MODIS, airborne multiangle Cloud Absorption Radiometer (CAR) (King et al., 1986[START_REF] Gatebe | Airborne spectral measurements of surface-atmosphere anisotropy for several surfaces and ecosystems over southern Africa[END_REF], Román et al., 2011), and two field-measured data sets [START_REF] Irons | Prediction and measurement of soil bidirectional reflectance[END_REF][START_REF] Deering | Characterization of the Reflectance Anisotropy of Three Boreal Forest Canopies in Spring Summer[END_REF]. We also explore the sensitivity of two hotspot parameters (C 1 and C 2 ) to hotspot-fits. Finally, we examined the influences of this new approach on the retrieval of intrinsic albedos and clumping index.

Hotspot data

POLDER-3 BRDF database

The spaceborne POLDER instrument can sample the land surface for viewing angles up to 60°-70°and for the full azimuth range, at a coarse spatial resolution of approximately 6 km.

Comprehensive BRDF sampling with large spatial coverage enables this instrument to collect observations for the development of BRDF modeling, particularly capturing distinct hotspot signatures. The POLDER-3 sensor onboard the Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar (PARASOL) satellite acquired multiangular measurements for more than 6 years, which were used to create four BRDF databases [START_REF] Bréon | A database of directional reflectance signature (GLC2000) with an analysis tool[END_REF]. The approach uses POLDER data from relatively homogeneous pixels to represent the dominant continental ecosystems. The data are geocoded, calibrated, atmospherically-corrected and cloudy-screened for acquiring the land surface bidirectional reflectance factors (BRFs) for each orbit. This study uses a monthly database containing 14,649 BRFs in 16 IGBP land cover classes. According to phase angle

The spatial distribution of view and Sun geometries for a typical POLDER data set, i.e., IGBP_01_20060609brdf_ndvi08.0824_1671, is exemplified in Figure 1 [START_REF] Bréon | the Cnes PARASOL Team, PARASOL Level-1 Product Data Format and User Manual[END_REF].

MODIS hotspot data

We Since the Terra and Aqua MODIS accumulate multiangular observations through multiple overlapping image swaths, the solar geometry corresponding to each observation will be different during a typical (16-day) retrieval period (e.g., approximately 20 for SZA in the h20v11 tile). This means that MODIS does not capture enough real-time observations in the principal plane (PP) under identical solar illumination conditions. Here, we adopt the assumption used with the MODIS BRDF/Albedo product, whereby the surface BRDF shape does abruptly for the range of SZAs capturing during a 16-day retrieval period. With this assumption, we can accumulate some multiangular observations in the proximity of the hotspot (i.e., )

to constrain the RTCLSR model for acquiring the optimal free parameters describing the hotpot magnitude (C 1 ) and width (C 2 ). Figure 1 (d) presents an extreme example that uses the RTCLSR model to reconstruct BRDF shapes as a function of phase angle for IGBP savanna in the NIR band. Red solid points represent observations extracted by using this method in the approximate PP. The black curve is the reconstructed BRDF shape using the RTCLSR model with C 1 =0.4, and C 2 =4.5°as the optimal values for the 2,275 MODIS data sets in this tile. The SZA was set to 30.58 , and the observation in the closest proximity to the hotspot direction (red point).

The BRDF shape reconstructed using the Maignan method (green curve) is provided here for visual comparison.

Figure 1 near here 2.3. Finer resolution measurements [START_REF] Gatebe | Airborne spectral measurements of surface-atmosphere anisotropy for several surfaces and ecosystems over southern Africa[END_REF], which captures hotspot signatures over 2 field sites. To acquire multiangular measurements, the CAR instrument is flown using a clockwise circular pattern above the surface repeatedly, and observes the reflected solar radiation at a fine angular resolution (i.e., 0.5 intervals through its 190°aperture at a rate of 100 scans per minute). This sampling scheme results in a BRDF retrieval that is based on 76,400 and 114,600 BRFs measurements per channel per complete orbit, which corresponds to a representative sample of the landscape-level (~5 km) reflected surface [START_REF] Gatebe | Airborne spectral measurements of surface-atmosphere anisotropy for several surfaces and ecosystems over southern Africa[END_REF]. At an altitude of 600m above the targeted surface area and 1 accuracy of CAR measurements is within an error margin of 0.3% (~2.0-3.5 meters as derived from the high resolution scene across the entire scan track). This accuracy holds well particularly with off-nadir looking observations [START_REF] Gatebe | Characterization of errors in the use of integrating-sphere systems in the calibration of scanning radiometers[END_REF]. These CAR measurements used in this study are averaged at an angular resolution of 1°in the viewing hemisphere and are taken in the red (0.682 µm) and NIR (0.870 µm) bands.

CAR data source was mainly from the early Smoke Clouds and Radiation-Brazil (SCAR-B) field campaign on August, 1995. Two kinds of data sets were collected from the well-defined surfaces of cerrado and dense forest [START_REF] Tsay | Airborne spectral measurements of surface anisotropy during SCAR-B[END_REF] To compare these models being explored at a field scale, we also analyzed two high-quality multiangular field data sets reported in previous studies (e.g. [START_REF] Li | A priori knowledge accumulation and its application to linear BRDF model inversion[END_REF][START_REF] Strugnell | A global albedo data set derived from AVHRR data for use in climate simula tions[END_REF]Huang & Jiao 2012;[START_REF] Jiao | An Anisotropic Flat Index (AFX) to derive BRDF Archetypes from MODIS[END_REF]. These include soil multiangular measurements (Iron et al., 1992) [START_REF] Deering | Characterization of the Reflectance Anisotropy of Three Boreal Forest Canopies in Spring Summer[END_REF] were also acquired with PARABOLA instrument with a 15°IFOV at a black spruce site that was mainly made up of old black spruce (picea mariana) with scattered emergent tamarack. The tree height was less than 10 m and the total stem density was 8040 live stems/ha, with a basal area of 40 m 2 /ha. Canopy closure averaged about 55 %, and the leaf area index (LAI) measured by an LAI-2000 in spring 1994 was 3.7 (Chen et al., 1997).

Model and method

RTLSR model

The semi-empirical, kernel-driven, linear BRDF model is a linear combination of three basic scattering components: isotropic scattering, volumetric scattering, and geometric-optical (GO) scattering. This model adopted a general form [START_REF] Roujean | A Bidirectional Reflectance Model of the Earth's Surface for the Correction of Remote Sensing Data[END_REF][START_REF] Lucht | An algorithm for the retrieval of albedo from space using semiempirical BRDF models[END_REF]:

( , , , ) ( ) ( ) ( , , ) ( ) ( , ,
)

v s iso vol vol v s geo geo v s R f f K f K (1)
Where f iso ( ), f vol ( ) and f geo ( ) are the spectrally dependent model parameters. K vol ( v , s ) and K geo ( v , s ) are kernel functions of view zenith v , illumination zenith s and relative azimuth and provide shapes for volumetric scattering and geometric-optical scattering BRDFs; f iso ( ) is a spectral constant for isotropic scattering; f vol ( ) and f geo ( ) are spectral constants, i.e., model anisotropic parameters that weight the two BRDFs; R( v , s , ) is BRDF in waveband .

K geo and K vol have been derived from physical approximation of the radiative transfer at the surface. K geo is derived from the GO model [START_REF] Li | Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy: effect of crown shape and mutual shadowing[END_REF][START_REF] Roujean | A Bidirectional Reflectance Model of the Earth's Surface for the Correction of Remote Sensing Data[END_REF] and characterizes dome-shaped BRDF curves. It is a function that describes the shadowing and surface scattering from the canopy. The operational MODIS BRDF/Albedo algorithm adopted the LiSparseReciprocal kernel (K LSR ) that was derived from a sparsely-vegetated canopy surface.

' ' ' ' ' sec sec cos 1 2 1 sec sec , , s v s v s v LSR O K (2)
Notably, where O( v , s ) is overlap function of view and illumination shadows on the ground, and hence the hotspot effect (i.e., the probability of seeing the sunlit background from the same gap) is included in this kernel. Further expressions and several intermediate variables are detailed to several papers (e.g., [START_REF] Wanner | On the derivation of kernels for kernel-driven models of bidirectional reflectance[END_REF][START_REF] Lucht | An algorithm for the retrieval of albedo from space using semiempirical BRDF models[END_REF].

Notably, two ratios (h/b and b/r) for describing vegetation structure are used to parameterize overlap function in K LSR . Here, h is the mean height where a crown center is located, b is the mean vertical half axis of the modeled ellipsoid, and r is the mean horizontal radius. These two ratios are related to the hotspot effects in K LSR on vegetation canopy scale [START_REF] Li | Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy: effect of crown shape and mutual shadowing[END_REF].

The operational MODIS RTLSR algorithm adopts h/b=2 and b/r=1. In Figure 2 (top), red dashed curve is for h/b=2 and b/r=1.2, and magenta dashed curve is for h/b=2.5 and b/r=1. Figure 2 demonstrates that, although these two ratios were in theory related to the general hotspot effect in K LSR , they do not seem to be very sensitive to the changes of hotspot effect at a reasonable range.

K vol is a similar function that describes the volumetric scattering component from canopy, based on an assumption of a single-scattering approximation of the radiative transfer (RT) theory by [START_REF] Ross | The Radiation Regime and Architecture of Plant Stands[END_REF]. The operational RTLSR model adopted RossThick kernel (K RT ) derived from a homogenously layered canopy with a large leaf area index (LAI) [START_REF] Wanner | On the derivation of kernels for kernel-driven models of bidirectional reflectance[END_REF] 4 cos cos sin cos 2

s v RT K (3)
This kernel characterizes bowl-shaped BRDF curves, but does not consider the correlation between illumination and observation processes for the observed hotspot. In theory, the probability of observing a sunlit component can be taken as the product of two probabilities: one for the view line from the observer and the other for illumination beam from Sun. When the view line and an illumination beam coincide, the observer either sees the sunlit foliage or the sunlit background, while the shadows of foliage and background are hidden from the view. As a result, the correlation tends to be 1:1, as the illumination and view directions get closer.

RTLSR hotspot correction

To consider the correlation between the view and illumination processes, [START_REF] Maignan | Bidirectional reflectance of Earth targets: evaluation of analytical mode ls using a large set of sp aceborne measurements with emphasis with the hot spot[END_REF] corrected the K RT with a hotspot factor:

4 / 1 1 cos cos sin cos 2 1 0 s v RTM K (4) 1 0 / 1 1
is a simplified hotspot factor derived by [START_REF] Bréon | Analysis of hot spot directional signatures measured from space[END_REF] from the theory on the calculation of an overlay function of the intersection of viewed and sunlit leaf areas [START_REF] Jupp | A hot spot model for leaf canopies[END_REF]. This hotspot factor was used to correct K RT , thereafter named K RTM . Here is phase angle, 0 is a characteristic angle in relation to the ratio of scattering element size and the canopy vertical density, which follows the range 0 =[1°,2°]. A 0 =1.5°has been suggested as a typical value representing a wide range of landscape conditions [START_REF] Bréon | Analysis of hot spot directional signatures measured from space[END_REF][START_REF] Maignan | Bidirectional reflectance of Earth targets: evaluation of analytical mode ls using a large set of sp aceborne measurements with emphasis with the hot spot[END_REF].

Note that the K RTM value near the hotspot direction will double rapidly (as much as K RT ) as the phase angle approaches 0°. To consider hotspot variability, we specified a 0 =3.0°for

K RTM for a comparison with the standard 0 =1.5°(Figure 2). The results indicate that an increase in 0 can increase hotspot width, but diverge from K RT beyond the hotspot region compared with K RTC . Therefore, we do not further examine the 0 influence on hotspot-fits in this study, rather than adopt 0 =1.5°as the most appropriate value. Here, the K RTM has a simplified feature whereby the variation in hotspot height and width is determined by phase angle .

For volumetric scattering kernel (K vol ), it is frequently required to meet the empirical requirement of K vol = 0 when both the viewing and illumination geometries point to nadir. This requirement aims to ensure that the isotropic parameter is physically reflectance for the model retrieval [START_REF] Roujean | A Bidirectional Reflectance Model of the Earth's Surface for the Correction of Remote Sensing Data[END_REF]. Conversely, using the Maignan method, the nadir-view and nadir-sun reflectance are specified as R(0,0,0,

) = f iso ( ) + f vol ( ) /4.
Here, we introduce the modified exponential function of the hotspot kernel function to K RT , hereby termed the RossThickChen kernel (K RTC ), which includes two free parameters characterizing hotspot variations [START_REF] Chen | A hotspot function in a simple bidirectional reflectance model for satellite applications[END_REF]:

2 1 cos cos sin cos 2 2 1 C s v RTC e C K (5) Where 2 1 1 C e C
is the modified hotspot function. The two free parameters, C 1 and C 2 , allow a large dynamic range of hotspot variation, and thus facilitate the analysis of the variation of hotspot height and width in fitting hotspot BRFs. For this retrieval, we also make an empirical adjustment to retain the isotropic BRDF model parameter, f iso , as the corresponding nadir-view and nadir-sun BRF.

Originally, the K RTC hotspot function was derived from a theoretical expression that accounts for the hotspot effect for a forest stand, based on a geometric-optical model [START_REF] Chen | A four-scale bidirectional reflectance model based on canopy architecture[END_REF]. In this specific scenario, C 1 , is linearly related to the difference between the reflectances of the foliage and background at the wavelengths of interest [START_REF] Chen | A hotspot function in a simple bidirectional reflectance model for satellite applications[END_REF].

In deriving this hotspot kernel, [START_REF] Chen | A four-scale bidirectional reflectance model based on canopy architecture[END_REF] demonstrated that this modification had very similar shapes at within-crown and between-crown scales because the gap size is scaled by the gap depth. For plane parallel layers of leaves, the average gap depth is taken as the mean distance between two layers of leaves and is related to the leaf area density. In this hotspot function (K RTC ), C 1 is related to the magnitude of the hotspot peak, while C 2 defines the half width of the hotspot in relation to the ratio of canopy height to the size of the predominant canopy structure. Notably for this study, to keep consistent with the original kernel form, the unit of C 2 is radians in terms of Eq. ( 5), but it is converted into degrees in subsequent use to describe the hotspot width.

Here, the K RTC kernel is empirically adjusted to meet K RTC = 0 for the nadir-view and nadir-sun geometries, which results in a downward shift of kernel shape (Figure 2, bottom). Note that this adjustment does not affect the fitting ability of this model since the kernel shape is retained. Figure 2

Hotspot parameter retrieval and analysis

The inversion strategy for retrieving the three parameters of this linear BRDF model is to minimize the root mean square error (RMSE) between model predicted and observed BRFs. The full-inversion is a simple matrix inversion that is independent of land surface types and is performed pixel by pixel using all high-quality measurements available. Details are referred to the papers [START_REF] Roujean | A Bidirectional Reflectance Model of the Earth's Surface for the Correction of Remote Sensing Data[END_REF][START_REF] Lucht | An algorithm for the retrieval of albedo from space using semiempirical BRDF models[END_REF][START_REF] Shuai | Quality assessment of BRDF/albedo retrievals in MODIS operational system[END_REF].

To derive the C 1 and C 2 values, we calculated fit-RMSEs to minimize the model-observation fits using accumulated measurements in a close vicinity to the hotspot .

Results

C 1 and C 2 values for POLDER

In this section, we assess the performance of the two hotspot parameters, based on the underlying IGBP class and two vegetation phenological phases, based on a comparison of the model predicted hotspot BRFs with observed BRFs derived from the POLDER BRDF database.

For this study, our analysis focuses on the RTCLSR and Maignan method. The C 1 values are significantly different between some IGBP classes; particularly in the red, and exhibits a high correlation between the red and the NIR (R 2 = 0.62). This reveals that a higher hotspot in the red band is frequently accompanied by a higher hotspot in the NIR band.

However, the C 2 parameter presents a low correlation between the two bands (R 2 = 0.21),

indicating that hotspot width is not band-dependent for these BRDF data.

Figure 3 (top right) shows the near-hotspot fit-RMSEs ( < 5°) between models in 6 bands.

Results indicate that the fit-RMSEs are lower in the visible than in NIR and shortwave infrared (SWIR). As expected, the RTLSR model had the largest fit-RMSEs in the vicinity of hotspot direction. We also derived globally optimized fit-RMSEs (RTCLSR_GLOB) for the entire POLDER data and the IGBP-bounded fit-RMSEs (RTCLSR_IGBP). The absolute average difference (AAD) of the fit-RMSEs (left ordinate) between Maignan and RTCLSR_IGBP models was 0.0023 per band, while the relative average difference (RAD) of fit-RMSEs (right ordinate) between these two models was ~20-30% in the blue, red and NIR bands [490, 670, 765nm]. This reveals that the RTCLSR model with the two free parameters can provide further improved hotspot-fits. Figure 3 (bottom) presents the fit-RMSEs for 15 IGBP classes between models in the red and NIR bands. Some improvements occur for several IGBP classes using the RTCLSR.

The improvements seem more pronounced for forest than for herbaceous classes in the red band, while the opposite appeared to be the case in the NIR band. are chiefly higher in the NIR than in the red, but do not show significant differences between these two seasons. This comparison reveals that NDVI alone doesn t distinctly capture this type hotspot variation as was detailed by using two hotspot parameters.

Figure 4 near here

Hotspot BRFs

The scatterplots in Figure 5 (top panel) show the agreement between model-predicted and observed hotspot BRFs in the red and NIR bands for the entire POLDER database. The correlations between the two are very high (R 2 > ~0.90) in both the red and NIR bands. In the red band, the Maignan method slightly underestimates the hotspot BRF; particularly at a range of low reflectances (i.e., bias = -0.005 for BRFs < 0.1), but overestimates hotspot BRFs at a range of high reflectances (i.e., bias = 0.011 for . In the NIR band, the Maignan method overestimates hotspot BRFs by 0.014 units at a range of 3, and by 0.008 at range of BRFs > 0.4. These biases are slightly corrected by using the RTCLSR model. Although the RTCLSR model with two free parameters fits the hotspot BRFs a little better than the Maignan method, the differences between them in fitting observed BRFs in the close proximity of hotspot were minor.

To further examine these two models in fitting hotspot BRFs as a function of surface type, we selected four IGBP classes to reflect a forest, mixed, grass, and sparsely vegetated gradient representing different canopy physiognomies and structures. In the red band, the Maignan method tends to underestimate the hotspot BRFs by 0.014 units for class. The underestimation is reduced to 0.011 units for (WSa) class, and the model fits well (underestimates by 0.003) for GrL) class;

but somewhat overestimates by 0.009 units for class.

The relative average difference (RAD) between these two model reaches to 12% for ENF class.

These biases are generally corrected by using the RTCLSR model (Figure 5). In the NIR band, the biases in reconstructing hotspot BRFs by Maignan method are -0.007 for ENF, 0.001 for WSa, -0.017 for GrL, and 0.018 for BSV. These biases are reduced to the range from -0.001 to -0.003 by using the RTCLSR model. This demonstrates that the RTCLSR model can provide further improvement for hotspot-fits for a subset of IGBP classes.

Hotspot parameter values for one tile of MODIS hotspot data

We retrieve two hotspot parameter values using one tile of MODIS hotspot data (h20v11) (Table 2). In general, C 1 values are larger in lower reflectance, which indicates relatively prominent hotspot effects due to lack of multiple scatterings from vegetation cover, which shows a consistent variation tendency with two hotspot parameter values for POLDER data. C 1 and C 2 values do not present significant differences between three major land cover types in this tile, dominated by savannas. However, the optimized C 1 and C 2 values present some differences in the corresponding bands between MODIS and POLDER, likely because MODIS hotspot data are derived only from one tile. From Figure 1 (d), we can also see that the BRDF shapes reconstructed using the Maignan (green curve) and RTCLSR models (black curve) in the NIR band are very close except for the hotspot peak. The optimal C 1 =0.4 for 2,275 MODIS data sets in this tile characterizes a low hotspot height, compared with the Maignan method for this extreme example. Note that the MODIS has difficulty in acquiring sufficient hotspot signatures on a global scale. As such, analysis of potential scale inconsistencies between MODIS and POLDER was not comprehensively performed using current hotspot data.

Results with in-situ and airborne data

We used airborne CAR SCAR-B cerrado and forest measurements [START_REF] Tsay | Airborne spectral measurements of surface anisotropy during SCAR-B[END_REF], as well as two high-quality field BRDF data sets (Iron et al., 1992;[START_REF] Deering | Characterization of the Reflectance Anisotropy of Three Boreal Forest Canopies in Spring Summer[END_REF] to evaluate the RTCLSR approach. Measurements are limited in the proximity of the principal plane, and according to our sampling design (Figure 1, a). The CAR cerrado dataset was collected on a forest-grass vegetation system that is known for its distinct hotspot signature in the red and NIR bands. Dense forest covered by tall trees with a close canopy had two-flight measurements. For these data sets, we only make use of the band where the hotspot signature was the most prominent.

Using CAR SCAR-B cerrado measurements, Figure 7 presents the model predicted and observed BRFs in the red and NIR bands (not shown for the RTLSR model). We provide three specific cases to identify the potential difference between these two methods. Figure 7 (a-b)

present the optimalized C 1 and C 2 for the RTCLSR model for this specific data set. Figure 7 (c-d) use the C 1 = 1.0 and C 2 = 3°as default values because K RTC with the default C values captures a very close kernel shape with K RTM (Figure 2). Figure 7 (e-f) present two-model consistencies by adjusting two hotspot parameters of the RTCLSR model to approach to Maignan predicted shapes for this specific data set (observations not shown). This comparison shows that the RTCLSR model with two free hotspot parameters can approach to Maignan method to a very great degree. Using CAR SCAR-B forest measurements, we compared the difference between the model-predicted and observed BRFs in the principal plane, and over the entire viewing geometries between two models (Figure 8). The difference in reconstructing the hotspot BRFs between these two models are distinctly exhibited over the entire viewing hemisphere. We further performed the statistical analysis for these two CAR datasets for BRFs near the hotspot hemisphere of VZA While these two approaches do not present difference for the entire CAR datasets, the significant improvement occurs in recreating hotspot signatures using the optimized hotspot parameters of the RTCLSR model for the RMSEs, biases and correlation coefficients (Figure 9). Since the BRDF sampling can have an effect on the sensitivity of model parameter retrievals for kernel-driven models [START_REF] Lucht | An algorithm for the retrieval of albedo from space using semiempirical BRDF models[END_REF], we also examined the sensitivity of the two hotspot parameters to the model parameters. To do this, we evaluated the modeled hotspot reflectance as a function of C 1 (taking C 2 = 3°as the default value) for three sample sizes (i.e., 12, 60 and 161 measurements) in principal plane. These BRDF observations are sampled from the airborne CAR cerrado measurements, which contained 29,160 BRF samples. In each case, 6

uniformly-distributed observations in the vicinity of the hotspot region (i.e., 5°) were used, and the rest of the observations outside of this hotspot region were randomly sampled. This investigation demonstrates that the sensitivity of hotspot parameters is somewhat related to varying BRDF sample numbers, which in turn tend to have an effect on model parameter retrievals. This effect holds true even when using the same observations in the vicinity of the hotspot 5°). This implies that the model predicted hotspot BRFs using a set of optimalized prior C 1 and C 2 values in RTCLSR model (combining K LSR and K RTC ) would be sensitive to the three model parameters (i.e., f iso , f geo and f vol ) to a certain extent that the adjustment of hotspot BRF dynamics by three model parameters is no longer effective. In such a situation, the two free hotspot parameters in RTCLSR model provide an improved capacity to capture more accurate hotspot signatures.

Figure 10 near here

Influence on retrieval of intrinsic albedo and clumping index

To evaluate the performance of the RTCLSR method for potential applications, we evaluated the retrieval of intrinsic albedos and clumping index (CI) from available POLDER and MODIS data. The BSA values were calculated for SZAs ranging from 0°to 75°at the interval of 15°.In general, the intrinsic albedo retrievals between the RCLSR and RTCLSR models are in a very high correlation for both POLDER (R 2 = ~1.0) and MODIS (R 2 > ~0.9) data; with negligible biases (< ~0.0003 in red and < ~0.0006 in NIR), while the difference in CI values between models was rather significant (Table 1 and Figure 12). As compared with the RTCLSR model results, a major overestimation of CI values occurs with the RTLSR model, but a slight underestimation of CI values occurs with the Maignan method. This is attributed to the difference in reconstructing hotspot BRFs between models. Figure 12 exhibits the scatter plots between models for ~11,632 CI values for 15 IGBP classes in the NIR (Figure 12, a-b) and the red bands (Figure 12, c-d) using POLDER data. According to the CI inversion algorithm, the CI biases between models (despite a high correlation, R 2 = ~0.9) result from biases in modeling hotspot BRFs. For most IGBP classes from POLDER data, about 20% of the relative average difference (RAD) in retrieving CIs occurs between the RTLSR model and the RTCLSR model, while about 10% of RAD occurs between the Maignan method and the RTCLSR in the NIR (Table 1). In the red band (Figure 12, c-d), the RAD in CI retrievals between the Maignan method and the RTCLSR model is generally small, while it still remains about 20% between the RTLSR and the RTCLSR. For the MODIS data, the RAD in CI retrievals arrive at about 6% and 12% in the red and the NIR bands between the RTCLSR and the Maignan (Figure 12, e-f).

29 Notably, we use more measurements in CI estimates because we can, in theory, reconstruct the hotspot BRFs using two prior IGBP-based hotspot parameter values that are derived from the high-quality observations in Section 4.1.1. In addition, we merely compare the retrieved CI values between models following the objective of this paper. To extensively validate these CI estimates with ground measurements would require the expanded use of additional measurements beyond the ones used in this study.

Conclusions and discussion

In this study, we proposed a new method to correct the hotspot effects for the RTLSR model known for its use within the operational MODIS BRDF/Albedo product. The method makes use of the exponential function with two free parameters (C 1 and C 2 ) to characterize hotspot height and width. Our results indicate that, although the Maignan method with no free hotspot parameters characterizes the hotspot effect in a relatively high accuracy (particularly for POLDER), the RTCLSR model can provide a further improvement in hotspot-fits. This is mainly attributed to the two free hotspot parameters that can be adjusted to reach their optimal values for the near-hotspot measurements available. Further analysis of the two hotspot parameters using hotspot data reveals that the hotspot signatures are somewhat related to surface type and vegetation phenology for available POLDER data. We found that the hotspot height (C 1 ) value is bigger in the red band than in the NIR band for most vegetation types, indicating a more prominent hotspot in the red band possibly due to the strong chlorophyll absorption of vegetation foliage.

The RTCLSR model is quite consistent with the other two models in the intrinsic albedos retrievals, but is somewhat different in CI retrievals through the use of the inversion algorithm of [START_REF] Chen | Global mapping of foliage clumping index using multi-angular satellite data[END_REF]. This reveals that the albedo retrievals using the RTCLSR would be consistent with the archived albedo product using the routine RTLSR algorithm; however, the RTCLSR model would provide more accurate hotspot signatures that may lead to an improved understanding of vegetation biophysical parameter retrievals in relation to hotspot BRFs. A similar analysis of the RTCLSR with the MODIS hotspot data, airborne CAR data and field measurements reveals the broad adaptability of this new method for different spatial resolutions; however, further investigation into the scale issue of hotspot effect using the RTCLSR model is still a challenge mainly due to the lack of sufficient hotspot data at different spatial resolutions.

Investigation into the sensitivity of these two hotspot parameters shows that C 1 is more sensitive to hotspot-fits than C 2 . The total BRDF sampling design can also play a role on the sensitivity of the hotspot parameters to hotspot-fits, because the total sample numbers (even having a good distribution) can have an effect on the model parameter retrievals, in particular the sparse BRDF sampling cases typified by the MODIS retrieval.

The original design of the kernel-driven linear Ross-Li BRDF model comprised a collection of kernels for different scenarios of land surface types. Thus, kernel functions are derived from different assumptions for vegetation canopy structures, and view and illumination geometries.

These assumptions may result in differences in modeling the radiation field [START_REF] Wanner | On the derivation of kernels for kernel-driven models of bidirectional reflectance[END_REF] especially in the hotspot direction (Huang & Jiao et al., 2012). Methodologically, the exponential approximation of this hotspot kernel function may also be used to correct the other volumetric scattering kernel (i.e., RossThin), because this exponential hotspot function makes use of two free parameters (C 1 and C 2 ) to fit hotspot signatures, rather than use surface biophysical parameters as inputs to drive a theoretical hotspot model. The exponential function form is rooted in gap probability theory and has been one of major contributions to hotspot modeling [START_REF] Qin | An evaluation of hotspot models for vegetation canopies[END_REF]. The corrected volumetric kernel can be combined with various geometric optical kernels for potential applications. However, the use of multi-kernel combination models would require recalibration of the C 1 and C 2 parameters.

In the situation where hotspot signatures are not available, the C 1 and C 2 values would need to be estimated on an a priori basis for it to initialize an RTCLSR model inversion. This has been attempted by fitting the hotspot data of several spatial resolutions via two means. Globally optimized C 1 and C 2 values can characterize the overall accuracy of the hotspot effect for certain satellite hotspot data, while variable-related optimized C 1 and C 2 values can help explore hotspot variation as a function of some underlying variables such as surface type and vegetation phenology. We found that a prior C 1 =1 and C 2 = 3°provide a stable initialized value in the RTCLSR model that is comparable in performance to the Maignan method.

Because the hotspot effect is very sensitive to the phase angle between the view and illumination in the retro-solar direction [START_REF] Bréon | Analysis of hot spot directional signatures measured from space[END_REF], its analysis for remote sensing applications requires a very high geometric accuracy. This implies that it is especially difficult to capture accurate hotspot observations in field measurements due to the shadows of the goniometer and sensors that must have a small enough IFOV (e.g., . Use of the field measurements with an IFOV= 15°in this study aims to stress that the reciprocal nature allows the RTCLSR model to acquire hotspot-fits with the least RMSEs for this data set. However, with two free hotspot parameters determined a priori by using enough hotspot data, this new method provides an improved understanding of the hotspot effect, and thus has potentials for certain ecological applications in regard to the hotspot BRFs for complex heterogeneous environments.

Appendix A

The shaded area in Figure 1

  extracted 2,275 MODIS hotspot data sets from one 10 × 10 tile (h20v11) of MODIS surface reflectance products (MOD09 and MYD09) at a spatial resolution of approximately 500-m, for 7 solar reflective bands, and using high-quality BRDF sampling distributions, represented by a phase angle range . There are only 80 high-quality hotspot data sets within a phase angle range of 1.5°; indicating that MODIS seldom acquires exact hotspot observations. Each data set consists of at least 10 observations with a proper directional sampling in approximate principal plane (PP), and includes fewer hotspot measurements within the defined area (Figure 1, a). The corresponding MODIS tile h20v11 captures a wide range of grass-shrub-savanna vegetation types in the Southern Hemisphere. The 329 th Day of Year (DOY) in 2010 represents a maturity season (Zhang et al., 2003). Solar zenith angles (SZAs) mainly ranged from 20°to 40°during the study period. The full-model inversion quality of the MODIS BRDF/Albedo product (MCD43A2) accounts for approximately 70 % of total pixels. The spatial distribution of MODIS observations and solar geometry from an extreme data set in this tile is presented in Figure 1 (c), where solid points are for MODIS observations, and circles are for solar geometries. The concentric circles are at 20°intervals.

  aligns K RTC with K RT and K RTM for a convenient comparison. The results indicate that K RTC coincides with K RT in the scattering direction beyond hotspot region, and adjusts the hotspot height and width by changing C 1 and C 2 (Figure 2 bottom). A larger C 1 value also indicates a higher hotspot, while a larger C 2 value indicates a wider hotspot.

Figure

  Figure 2 near here

  Where R obs ( v , s , ) and R model ( v , s , ) are the model predicted and observed BRFs near the hotspot region in the viewing and solar geometries, v , s at wavelength , as in Equation (1). We derived the optimized C 1 and C 2 values from the least RMSEs using the iteration step size of 0.1 for C 1 = [0.3,1.2] and C 2 = [1.0°,6.0°]. This range of C 1 and C 2 can ensure realistic hotspot variations while optimizing processing time. Following this method, we retrieved the C 1 and C 2 values for the entire POLDER database; including parameterizations by IGBP cover type and phenological phase, respectively. To derive the optimal hotspot parameter values, we constrain this model by accumulating all measurements near hotspot direction ( ) for these two parameter types, including all 15 IGBP classes, excluding and 2 phenological stages: (1) maturity and (2) dormancy. We then compared the model predicted and observed BRFs near hotspot direction by using scatterplots with regression lines. To display model-observation fits, we examined the model predicted and observed BRFs as a function of phase angle for the selected observations for Evergreen Needleleaf Forest type (ENF) . Details on how to derive the shaded area (Figure 1, a) are provided in Appendix A. We also performed a comparison against finer resolution in-situ data. To analyze the sensitivity of the hotspot parameters to hotspot-fits in conjunction with the model parameters, we used simulated BRDF shapes generated from the airborne CAR data set. Finally, we examined the intrinsic albedos (i.e., White Sky Albedo (WSA) and Black Sky Albedo (BSA)) and clumping index (CI) retrievals between models. To retrieve CI values, we made use of the algorithm by Chen et al. (2005) and assigned a simplified geometry of 45°view zenith angles in both the backward-and forward-scattering directions for the corresponding principal plane hotspot and dark spot, respectively (Zhu et al. 2012). Such a simplification is somewhat different from Chen et al. (2005), which focused on locating the optimal view geometry to capture the dark spot. Since the three models are generally consistent in fitting observations outside of the hotspot region, this simplification would not affect the inter-comparison of CI retrievals between models.

Figure 3 (

 3 Figure 3 (top left) presents the globally optimized C 1 and C 2 values (RTCLSR_GLOB) for 6

Figure 3

 3 Figure 3 near here

Figure 5

 5 Figure 5 near here

Figure 1

 1 Figure 1 (a) to approximate the principal plane. Data sources are marked on each subplot with a

Figure

  Figure 6 near here

Figure 7 ,

 7 Figure 7, 8 and 9 near here

Figure

  Figure 10 near here

Figure 11 (

 11 Figure 11 (bottom) presents the sensitivity of modeled hotspot reflectance as a function of

  (a) contains four parts. As the 1 st constraint, all observations falling within the phase angle having the radius of R should be included in set A. Here, R but R can be adjusted properly for different hotspot data.R ATo select the observed BRFs close to principal plane (PP) for visually comparing with modeled BRFs in PP, we define the distance, H, which is perpendicular to the principal plane, as the 2 nd constraint. Obviously, a less H value can select observations in the shaded area to get closer to PP. This generates set B as follows. Here, As the 3 rd constraint, we should consider that, in a 2-D plot that exhibits the observed BRFs in approximate PP, spurious visualizations, particularly in the proximity of hotspot direction in PP, possibly result from some observed BRFs that are not actually near hotspot region in PP, e.g., a H taking Sun as the center (Figure1, a) to be exhibited as near-hotspot BRFs in PP. To remove these observations in 2-D plots that present observed and modeled BRFs in approximately PP (e.g., Figure 1, d), we define an intersection angle ( ) between PP line and the line passing through Sun (Figure 1, a), and derive set C and set D as follows. Here, these hotspot data used in this study. shaded area in Figure 1 (a) can be derived by implementing set operation for four sets above: A B C D

Figure 1 .Figure 2 .

 12 Figure 1. (a) The method to select observations near principal plane and in the proximity of hotspot

Figure 3 .

 3 Figure 3. The globally optimized C 1 and C 2 values, and the fit-RMSEs derived from the entire POLDER

Figure 4 .

 4 Figure 4. Two hotspot parameters (top) as a function of IGBP class in terms of the timing of Jun-Jul-Aug/

Figure 5 .

 5 Figure 5. Comparison of the model predicted and observed hotspot BRFs for the entire POLDER BRDF

Figure 6 .Figure 7 .

 67 Figure6. POLDER observatoins (red points) and the reconstructed BRDF shapes byMaignan (green)

Figure 8 .

 8 Figure 8. Using CAR SCAR-B forests to examine the difference between model predicted and observed

Figure 9 .Figure 10 .Figure 11 .

 91011 Figure 9. Scatterplots showing the difference between modeled and observed BRFs near the hotspot

Figure 12 .

 12 Figure 12. Comparison of CI retrievals between models using POLDER data in the NIR band (top) and in

Figure 1 .Figure 2 .

 12 Figure 1. (a) the method to select observations near principal plane and in the proximity of hotspot

Figure 3 .

 3 Figure 3. The globally optimized C 1 and C 2 values, and the fit-RMSEs derived from the entire POLDER

Figure 4 .

 4 Figure 4. Two hotspot parameters (top) as a function of IGBP class in terms of the timing of

Figure 5 .

 5 Figure 5. Comparison of the model predicted and observed BRFs for the entire POLDER BRDF data and two IGBP classes, i.e., ENF and WSa in the red band.

Figure 6 .

 6 Figure 6. POLDER observatoins (red points) and the reconstructed BRDF shapes by Maignan (green)

Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .

 7891011 Figure 7. Comparison of Maignan with RTCLSR model using CAR/SCAR-B cerrado measurements in

  (b). Since Snow and Ice class tends to have maximum reflectance in forward scattering direction that current surface BRF models cannot characterize, we exclude this class from the database. The BRDF signatures are based on inversions of approximately 20 orbits worth of spectral BRF data, for each of the six POLDER bands. Up to 16 different multiangular measurements for a given POLDER pixel are included in each orbit. BRFs measured at two typical bands (red, centered 670 nm and NIR, centered 865 nm) are mainly selected to assess the hotspot effect reconstructed from the

proposed model. Notably, POLDER-3 spectral measurements are not simultaneous, meaning that each channel is acquired with a slightly different viewing geometry. As such, to analyze the variability of hotspot signatures, this requires a careful assessment of measurement differences, including normalization of the POLDER viewing geometry

Table 1 .

 1 Statistics of the clumping index values retrieved by three models

	Class	Number	RTLSR Mean Std. Mean Std. Mean Std. Maignan RTCLSR (RTLSR-RTCLSR) / RTCLSR (%)	(Maignan-RTCLSR) /RTCLSR (%)
	ENF	793	0.57 0.04 0.49 0.03 0.51 0.04	13.09	4.32
	EBF	898	0.76 0.08 0.56 0.06 0.64 0.07	19.02	11.85
	DNF	219	0.66 0.02 0.57 0.02 0.59 0.02	10.54	4.29
	DBF	682	0.80 0.06 0.63 0.05 0.67 0.05	19.59	6.89
	MiF	690	0.76 0.10 0.57 0.08 0.63 0.09	21.55	9.02
	CSh	549	0.77 0.08 0.58 0.06 0.63 0.07	22.83	7.22
	OSh	1247	0.86 0.10 0.68 0.10 0.74 0.10	16.31	8.68
	Wsa	1035	0.83 0.06 0.62 0.06 0.68 0.06	22.49	8.32
	Sav	734	0.87 0.06 0.67 0.06 0.74 0.06	17.64	10.11
	GrL	1075	0.88 0.09 0.68 0.10 0.76 0.10	16.21	9.56
	Pwe	21	0.95 0.06 0.65 0.05 0.77 0.05	22.29	15.68
	CrL	960	0.89 0.07 0.67 0.06 0.75 0.06	19.55	10.48
	Ubu	835	0.84 0.06 0.65 0.06 0.72 0.06	17.20	10.13
	CNVM 727	0.89 0.06 0.66 0.05 0.75 0.05	19.00	10.94
	BSV	1167	0.99 0.16 0.85 0.19 0.92 0.17	8.00	7.44
	Table 2. The globally optimized C 1 and C 2 values derived from one 10 × 10 tile (h20v11) of MODIS
		surface reflectance products (MOD09 and MYD09) in 7 reflected solar bands
	Band	Band1	Band2	Band3	Band4	Band5		Band6	Band7
	(nm)	(620-670)	(841-876)	(459-479)	(545-565)	(1230-1250)	(1628-1652)	(2105-2155)
	C 1	0.7	0.4	0.7	0.7	0.6		0.7	0.7
	C 2 (°)	5.2	4.5	5.2	5.2	3.5		5.2	5.2

Acknowledgments

Francois-Marie Bréon (LSCE, France) for his kindly providing a piece of code to reprocess the POLDER data, and for his great comments and suggestions on this paper.


Mich Ele Sebag
email: michele.sebag@polytechnique.fr

Antoine Ducoulombier
email: antoine.ducoulombier@lri.fr

Extending Population-Based Incremental Learning to Continuous Search Spaces

An alternative to Darwinian-like arti cial evolution is o ered by Population-Based Incremental Learning (PBIL): this algorithm memorizes the best past individuals and uses this memory as a distribution, to generate the next population from scratch. This paper extends PBIL from boolean to continuous search spaces. A Gaussian model is used for the distribution of the population. The center of this model is constructed as in boolean PBIL. Several ways of de ning and adjusting the variance of the model are investigated. The approach is validated on several large-sized problems.

Introduction

Evolutionary algorithms (EAs) [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF][START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF][START_REF] Fogel | Evolutionary Computation. Toward a New Philosophy of Machine Intelligence[END_REF] are mostly used to nd the optima of some tness function F de ned on a search space .

F : ! IR From a machine learning (ML) perspective 9], evolution is similar to learning by query: Learning by query starts with a void hypothesis and gradually re nes the current hypothesis through asking questions to some oracle.

In ML, the sought hypothesis is the description of the target concept; the system generates examples and asks the oracle (the expert) whether these examples belong to the target concept. In EA, the sought "hypothesis" is the distribution of the optima of F; the system generates individuals and asks the oracle (a routine or the user) what their tness is. In all cases, the system alternatively generates questions (examples or individuals) depending on its current hypothesis, and re nes this hypothesis depending on the oracle's answers.

One core di erence between ML and evolution is that ML, in the arti cial intelligence vein, manipulates high-level, or intensional description of the hypothesis sought. Conversely, evolution deals with a low-level, or extensional description of the sought distribution: the distribution of the optima is represented by a collection of individuals (the current population).

The Population Based Incremental Learning (PBIL) approach bridges the gap between ML and EAs: it explicitly constructs an intensional description of the optima of F, expressed as a distribution on 2, 3]. This distribution is alternatively used to generate the current population, and updated from the best individuals of the current population. The advantage of the approach is that, as claimed throughout arti cial intelligence 12], the higher level the information, the more explicit and simple the information processing can be. And indeed, PBIL involves much less parameters than even the canonical GAs 6].

PBIL was designed for binary search spaces. It actually constructs a distribution on = f0; 1g N represented as an element of 0; 1] N . The basics of this scheme are rst brie y recalled in order for this paper to be self contained (section 2). Our goal here is to extend this scheme to a continuous search space IR N . Continuous PBIL, noted PBIL C , evolves a Gaussian distribution on noted N(X; ). The center X of the distribution is evolved much like in the binary case; evolving the standard deviation of this distribution is more critical, and several heuristics to this aim are proposed (section 3). PBIL C is nally validated and compared to evolution strategies on several large-sized problems (section 4). The paper ends with some perspectives for further research. Let denote a population of individuals in = f0; 1g N . An element h of H = 0; 1] N can be associated to , by de ning h i as the fraction of individuals in having their i-th bit set to 1. Conversely, an element h in H de nes a distribution over : one draws an element X = (X 1 ; : : :; X N ) in by setting X i to 1 with probability h i .

PBIL relies on the following premises 2]: a) if evolution succeeds, the population converges toward a single1 optimum of F; b) the more converged the population , the better it is represented by h. Assuming these, PBIL discards all information in the population not contained in h: The population is simply considered as a manifestation of h. The attention is thus shifted from evolving by means of mutation and recombination, to evolving h (Fig. 1). To this aim, PBIL uses the information contained in the current population t : h is evolved, or rather updated, by relaxation from the best individual X max in t : h t+1 = (1 ) : h t + : X max ; in ]0; 1 Distribution h t can be viewed as the memory of the best individuals generated by evolution. Relaxation factor corresponds to the fading of the memory: the higher , the faster h t moves toward the current local optimum.

In contrast to standard evolution, PBIL explicitly explores the space H of distributions on . And, as noted already, this higher level representation allows for a simpler information processing: besides the population size, PBIL involves a single key parameter, , to be compared to the various parameters controlling mutation and recombination. Further, the exploration is deterministic, in the sense that h t is deterministically updated from the current population2 .

Discussion

Let us reformulate PBIL as a learning-by-query algorithm, by de ning a partial generality order on the set of distributions H. The generality of a distribution h is clearly related to the diversity of the population generated from h, and the diversity of the population with regard to bit i is inversely proportional to jh i :5j. Accordingly, a distribution h is more speci c than h 0 , if, for each bit i, either 0 h i h 0 i :5, or :5 h 0 i h i 1.

PBIL initializes h to the most general distribution h 0 = (:5 : : : ; :5), and gradually specializes it along generations. Let X h denote the (boolean) individual most similar to h t ; then, h t is specialized on all bits i such that X h i = X max i .

The complete convergence of the scheme is avoided as h t i never reaches 0 or 1;

in theory, PBIL can generate any individual at any time.

In practice, PBIL can su er from premature convergence. This happens when h t gets too speci c3 , and no new good individual is discovered. PBIL o ers two heuristics to resist premature convergence 2]: Using the average of the two best individuals in t , rather than the single best one. This way, h t is generalized on all bits discriminating these individuals. Perturbing h t with a Gaussian noise: with a given probability (5%), a Gaussian variable with a low standard deviation is added to h t i . This way, the center of the distribution is durably perturbed, which helps escaping from local minima.

A more fundamental limitation of PBIL comes from the distribution space, which implicitly assumes the linear separability of the problem (genes are considered independent). This distribution space appears too poor to t complex tness landscapes, such as the Long Path problem 7]. Previous experiments show that distributions used in PBIL have di culties to overlap the narrow path 14]. Recent extensions to PBIL have considered richer distribution spaces 4].

3 Continuous PBIL This section rst brie y discusses a previous attempt to extend PBIL to continuous search spaces, then details the proposed method and outlines PBIL C .

Continuous PBIL with dichotomic distributions

To the best of our knowledge, the only extension of PBIL to continuous search spaces has been proposed in 15]. This algorithm explores the search space much like the delta-coding approach 17]. The domain of each gene is divided into two intervals ("low" and "high" values); the current distribution h (h in 0; 1] N ) is used to determine which interval an individual belongs to:

X i 2 a; b] Probability(X i > a + b
2 ) = h i X i is then drawn with uniform probability in the selected interval.

At each generation, h is updated like in the boolean case, by memorizing whether the best individual takes low or high values for each gene:

h t+1 i = (1 ) : h t i + : (X max i > a + b 2 )
When h i gets speci c enough (h i < :1 or h i > :9), the population gets concentrated in a single interval (resp. a; a+b 2 ] or a+b 2 ; b]). The search is then focused: the domain of the gene is set to the interval considered and h i is reinitialized to :5.

In this scheme, evolution gradually focuses on the region most often containing the best individuals. One limitation is that a region which has been discarded at some point is hardly explored ever after, and this violates the ergodicity requirement. Furthermore, the search might be insu ciently focused, given the poor (uniform) distribution used within the selected interval.

Continuous PBIL with Gaussian distributions

Our approach rather explores Gaussian distributions N(X; ) on the search space , given as products of Gaussian distributions N(X i ; i ) on each gene domain. With no loss of generality, is set to 0; 1] N in the following.

Like PBIL, PBIL C starts with a rather general distribution; then it alternatively uses this distribution to draw the population, and uses the population to update the distribution. The center of the distribution X t is initialized to the center of the search space (:5; : : : ; :5). At each generation, X t is updated from a linear combination of the two best and the worst individuals in the current population, inspired from PBIL and Di erential Evolution 16]: X t+1 = (1 ) : X t + : (X best; 1 + X best; 2 X worst ) The diversity of the population, controlling the convergence of evolution, depends on the variance = ( 1 ; : : : N ) of the distribution. Several heuristics have been investigated to adjust parameters i . A The simplest possibility is to use a constant value. The trade-o between exploration and exploitation is thus settled once for all: the search cannot become too speci c and it cannot be speeded up either. B A second possibility is to make evolution itself adjust . PBIL C here proceeds exactly as a self-adaptive (1; )-evolution strategy (ES) 4 where stands for the size of the population, except that the parent is replaced by the center X t of the distribution. C A third possibility is to adjust depending on the diversity of the current best o spring; t is then set to the variance of the K best current o spring:

i = s P K j=1 (X j i X i ) 2 K
where X denotes the average of the best K o spring X 1 ; : : : X K . D Last, can be learned in the same way as X itself, by memorizing the diversity of the K best o spring:

t+1 i = (1 ) t i + s P K j=1 (X j i X i ) 2 K 3.

Discussion

At rst sight, PBIL C is quite similar to a (1; )-ES, the o spring being generated from the single parent (X t ; t ). The di erence is twofold.

In (1; )-ES, the parent is simply replaced by the best o spring, whereas PBIL C updates X t by relaxation. Let any o spring X k be written X t + Z k , with Z k being a random vector drawn according to N(0; t ). Then it comes: X t+1 = (1 )X t + (X best;1 +X best;2 X worst ) = X t + (Z best;1 +Z best;2 Z worst )

The evolution of X t can be viewed as a particular case of weighted recombination as studied by Rudolph 11]; a theoretical analysis shows that weighted recombination with optimal weights should be preferred to the simple replacement of the parents. Interestingly, the heuristic recombination used in PBIL C is intermediate between two particular cases with good theoretical properties (for F(X) = P X 2 i ): the half sum of the two best o spring, and the di erence of the best and the worst o spring. PBIL C uses xed, hence non-optimal, weights; but note that intervenes as an additional scaling factor, controlling the variance of X t .

Independently, the variance of X t is also controlled from t . PBIL C uses global mechanisms (options A, B and D) to adjust t , by opposition to the local adjustment of achieved by self-adaptive mutation. Actually, the adjustment of (option D) much resembles the 1/5th rule used to globally adjust in early evolution strategies 10]. The di erence is that the 1/5th rule criterion compares the o spring to the parents, and considers whether a su cient fraction of o spring is more t than the parents. In opposition, PBIL C only examines the diversity of the best t o spring: it does not need to restrict the exploration, even if the o spring are less t than the parent, because the center of the explored region moves more slowly than in standard ES.

To sum up, PBIL C controls the exploration-exploitation tradeo in a way rather di erent from that of (1; )-ES. First of all, the single parent does not jump directly to a desirable location (the best o spring, or some weighted combination of the remarkable o spring), but rather makes a very small step toward this desirable location (e.g. is set to 10 2 in the experiments). Variance is adjusted in a similarly cautious way. It appears that ES takes instant decisions, on the basis of the instant information. On the opposite, PBIL C maintains a long-term memory, slowly updated from the instant information, and bases its cautious decisions on this long-term memory.

Validation

This section describes the goal of the experiments and the problems considered. We then report and discuss the results obtained.

Experiment Goals and Problems

Our goal is to study the respective advantages of evolving extensional vs intensional information about the tness landscape. Practically, PBIL C , evolving an intensional information represented as a distribution, is compared to selfadaptive evolution strategy, evolving an extensional information represented as usual as a population. We deliberately consider large-sized search spaces (N = 100) for the following reason. In low or middle-sized spaces, populations or distributions might convey similarly accurate information about the tness landscape. This is not true in large-sized spaces: any reasonable number of point s can only convey a very poor information about IR 100 . Experimenting PBIL C in IR 100 will show how intensional evolution stands the curse of dimensionality. Functions and search spaces considered are displayed in Table 1. Functions F 1 to F 3 have been used to evaluated binary PBIL 2]. Besides the size of the search space, F 1 and F 2 su er from an additional di culty, epistasis (the genes are linked via the y i ). Functions F 6 to F 8 have been extensively studied in the literature, for lower-sized search spaces (N 30).

Notation

Experimental setting

We used two reference algorithms: boolean PBIL working on a discretization of the continuous problem (each continuous variable is coded through 9 binary variables), using either a binary or a Gray coding; and a (10 + 50)-ES with self adaptive mutation 1]. In the PBIL case, the size of the population is set to 50 and the relaxation factor is set to :01.

PBIL C involves the same setting as PBIL ( = 50 and = :01). Four options regarding the variance of the distributions have been considered (section 3.2): A Constant variance. B Self-adapted variance: PBIL C here behaves like a self-adaptive (1; )-ES, except that the parent is replaced by X t . C Instant variance: i is set to the variance of the best K o spring in the population. Several values of K were considered: =2; =3; =5. D Relaxed variance: i is the variance of the best K o spring relaxed over the past generations; the relaxation factor is again set to = :01.

Results

Algorithm F1 F2 F3 Table 2 displays the results obtained on functions F 1 ; F 2 and F 3 . Results obtained by boolean PBIL are taken from 2]; additional results not reported here, show that boolean PBIL signi cantly outperforms several variants of GAs and Hill-Climbers on these functions. Note that all algorithms end rather far from the actual optimum (10 7 ). Still, PBIL C signi cantly outperforms standard ES on these problems | provided that the variance of the distribution is adequately set. Note also that PBIL C outperforms PBIL itself, working on a binary or Gray discretization of these continuous problems. This might be due either to the loss of information entailed by discretization, or because PBIL, as already mentioned, explores a too restricted distribution space.

The worst results of PBIL C are obtained when is self-adapted or set to the diversity of the current best o spring (options B and C); they are due to a fast decreasing of . And, in retrospect, a vicious circle occurs when tightly depends on the diversity of the o spring: the less diverse the o spring, the smaller , hence the less diverse the o spring... Setting to a constant value (option A; the particular values were chosen after 10,000 evaluations preliminary runs) leads to satisfactory results, even outperforming those of standard ES. Further experiments will show whether this is rather due to the superiority of weighted recombination (replacing a parent by a combination of o spring) over replacement | or to the "long-term memory" e ect, as the parent slowly moves toward the weighted combination of the o spring instead of jumping there.

The best option appears to learn the variance in the same way as the center of the distribution X t (option D). Further, the fraction K of the o spring considered to update apparently is not a critical parameter 5 .

Algorithm

F6 F7 F8 These trends are con rmed by preliminary experiments on F 6 ; F 7 and F 8 (Table 3): PBIL C signi cantly outperforms self-adaptive ES on two out of the three problems, the best option for adjusting being the relaxation from a small fraction of the best o spring.

Conclusion

The main originality of PBIL is to reformulate evolution into new, higher-level, terms: rather than specifying all operations needed to transform a population into another population (selection, recombination, mutation, replacement), one only speci es how to evolve or update a distribution given the additional information supplied by the current population. At this level, many core traits of evolution (e.g. diversity, speed of changes) are explicit and can be directly controlled.

Overall, evolution shifts from the stochastic exploration of the search space , to learning a distribution on by reinforcement from the current population. This paper extends PBIL from boolean to continuous search spaces, by learning Gaussian distributions N(X; One argument for learning distributions is that it expectedly scales up more easily than evolving populations: a reasonable size population gives little information on large-sized search space. Experimental results on large-sized problems show that PBIL C actually outperforms standard ES on ve out of six problems (with one or two orders of magnitude) and also outperforms the original PBIL working on a discretized version of the continuous problems considered.

Nevertheless, given the size of the search space, PBIL C ends rather far from the optimum on four out of six problems. Further experiments will consider other problems, and study how PBIL C behaves in the last stages of exploitation. Another perspective of research is to evolve several distributions rather than a single one. This would relax the main limitation of the PBIL scheme, that is, the fact that it can only discover a single optimum. Indeed, learning simultaneously several distributions is very comparable to evolving several species. The advantage is that comparing an individual to a few distributions might be less expensive and again more transparent, than clustering the population, adjusting the selection or the tness function to ensure the co-evolution of species.

Figure 1 .

 1 Figure 1. Comparing the generation steps in standard evolution and PBIL

Table 1 :

 1 Fitness functions considered,

	F 1 F 2 F 3 F 6 F 7 F 8	De nition with yi = xi + yi 1,i y1 = x1 with y1 = x1 yi = xi + sinyi 1, i with yi = :024 (i + 1) i (x 2 i Acos(2 x i )) + 100A i jyij P 100 10 5 + 100 10 5 + P i jyij 100 10 5 + P i jyij P P i x i sin p x i P i x 2 i Q i cos( xi p i+1 )	2	2 xi	Domain 3; 3] 100 3; 3] 100 3; 3] 100 5; 5] 100 30; 30] 100 100; 100] 100

i = 1 : : : 100

Table 2 :

 2 Best Fitness (averaged on 20 runs) for 200,000 evaluations

	Best results indicated in bold	Exact optimum of F1 ; F2 and F3 = 10 7

Table 3 :

 3 Best Fitness (averaged on 20 runs) for 200,000 evaluations

	Best results indicated in bold	Exact optimum of F6 and F8 = 0

  ). The resulting PBIL C algorithm can be thought of as a (1; )-ES, with the following di erences. ES takes instant decisions, on the basis of the instant information. PBIL C maintains a long-term memory, takes its decisions on the basis of this long-term memory, and slowly updates the memory from the instant information. Practically, the parent of a (1; )-ES jumps toward the best o spring; in opposition, the center of the distribution in PBIL C cautiously moves toward a weighted combination of the o spring. Similarly, self-adaptive ES locally adjusts the variance of mutation by means of instant decisions; in opposition, PBIL C cautiously updates the variance from the global diversity of the best o spring.

This claim obviously makes no room for diversity preserving mechanisms 8].

This raises the question of whether PBIL still pertains to the eld of "Parallel problem solving from nature": is "nature" allowed to explicitly manipulate distributions? Still, a de nition of "nature" is clearly beyond the scope of the paper.

[START_REF] Baluja | Removing the genetics from the standard genetic algorithms[END_REF] Parameter partly controls the speci city of h t , and plays the same role as selection in GAs: the diversity decreases, everything else being equal, as goes to 1.

In self-adaptive ES, besides the Xi an individual X carries the variance i of the mutation to be applied on the Xi 13, 1]: Mutation rst evolves the i, then uses the new i to perturb the Xi. Evolution thus hopefully adjusts the i "for free", at the individual level.

This holds for all problems except F3, which is the problem with most diversity in the tness of the o spring. This might be an indication for choosing K adaptively: e.g. retain the o spring whose tness is greater than a given function of the average tness and deviation of the tness in the current population.

Acknowledgments

Many thanks to Marc Schoenauer, for many valuable comments, and to the second anonymous referee, for very insightful comments and suggestions.