Evaluation of Medium Spatial Resolution BRDF-Adjustment Techniques Using Multi-Angular SPOT4 (Take5) Acquisitions

Abstract

High-resolution sensor Surface Reflectance (SR) data are affected by surface anisotropy but are difficult to adjust because of the low temporal frequency of the acquisitions and the low angular sampling. This paper evaluates five high spatial resolution Bidirectional Reflectance Distribution Function (BRDF) adjustment techniques. The evaluation is based on the noise level of the SR Time Series (TS) corrected to a normalized geometry (nadir view, 45° sun zenith angle) extracted from the multi-angular acquisitions of SPOT4 over three study areas (one in Arizona, two in France) during the five-month SPOT4 (Take5) experiment. Two uniform techniques (Cst, for Constant, and Av, for Average), relying on the Vermote–Justice–Bréon (VJB) BRDF method, assume no variation in space of the BRDF shape. Two methods (VI-dis, for NDVI-based disaggregation and LC-dis, for Land-Cover based disaggregation) are based on disaggregation of the MODIS-derived BRDF VJB parameters using vegetation index and land cover, respectively. The last technique (LUM, for Look-Up Map) relies on the MCD43 MODIS BRDF products and a crop type data layer. The VI-dis technique produced the lowest level of noise corresponding to the most effective adjustment: reduction from directional to normalized SR TS noises by 40% and 50% on average, for red and near-infrared bands, respectively. The uniform techniques displayed very good results, suggesting that a simple and uniform BRDF-shape assumption is good enough to adjust the BRDF in such geometric configuration (the view zenith angle varies from nadir to 25°). The most complex techniques relying on land cover (LC-dis and LUM) displayed contrasting results depending on the land cover.

1. Introduction

Land surface anisotropy arising from non-constant observation and illumination geometries can cause significant variations in the surface reflectance unrelated to surface changes. These variations, which cause time series noise, can affect land surface analysis. Based on POLDER observations, Bréon et al. [1] have determined that the ratio of minimum to maximum surface reflectance due to sun-view geometry differences can be up to five in the visible domain. During the past 20 years, these variations have been deeply analyzed and characterized using coarse resolution optical sensors that are usually associated with a large field of view (higher than 110°, e.g., Multi-angle Imaging SpectroRadiometer (MISR), Moderate Resolution Imaging Spectroradiometer (MODIS), MEdium Resolution Imaging Spectrometer (MERIS), and Visible/Infrared Imager/Radiometer Suite (VIIRS)).

The anisotropic reflectance variations are quantified using Bidirectional Reflectance Distribution Function (BRDF) models, which depend on coefficients describing the surface (isotropic and anisotropic components) and on the observation and illumination geometry, described by the view zenith, the sun zenith and the sun-view relative azimuth angles. Most of the BRDF-retrieval methods rely on the inversion of BRDF model coefficients using several observations over a limited time period (e.g., 16 days for MODIS BRDF products, MCD43) of the same area from various sun-view geometries: examples include MODIS [2,3], MISR [4] and POLDER [5]. In 2009, a new approach to adjust BRDF effects from MODIS SR data was introduced Vermote Justice Bréon (VJB) model, [6]. Compared to previously cited methods, the VJB method supports some simplification hypotheses that relate BRDF variation to Normalized Difference Vegetation Index (NDVI) variation.

The launch of the Sentinel-2 satellite with the MultiSpectral Instrument (MSI) sensor will initiate a new era in terms of land surface analysis relying on high spatial and temporal resolution data [7]. Time series derived from MSI, which can be complemented by other high spatial resolution sensors data, bring to users new capabilities for vegetation monitoring at a fine scale [8]. However, few studies have addressed the question of compensating for the BRDF observed with high spatial resolution sensors (10–30 meters). Indeed, these sensors do not have sufficient angular and temporal dynamics to characterize the BRDF directly as it has been done with large field of view sensors. The BRDF effects remain quite limited for near-nadir and narrow field of view sensors (e.g., Landsat (15° field of view)), but they can cause higher SR variations for wider swath sensors (e.g., Disaster Monitoring Constellation (DMC), MSI) or with off-nadir pointing capability sensors (e.g., sensors on-board SPOT (Satellites Pour l'Observation de la Terre) constellation). Retrieving a spatially explicit BRDF from such sensors has to rely on a-priori estimations of surface anisotropy.

As of 2015, there are several existing approaches to model the BRDF effect at high spatial resolution. Roy et al. [9] used directly MODIS BRDF parameters at 500 m to adjust the BRDF effects for Landsat imagery. Li et al. [10] utilized regional BRDF information from MODIS BRDF products and corrected angular effects for Landsat scenes. Shuai et al. [11] developed an approach to extract BRDF parameters for different surface types based on MODIS pure pixels and used this information to correct angular effects for the corresponding surface types. Based on [11], Franch et al. [12] have adapted an alternative approach to invert the BRDF parameters without relying on pure pixel and using the 0.05° MODIS BRDF information retrieved with the VJB model [6]. Other approaches rely on the use of pre-defined BRDF shapes, often called archetypes that capture anisotropic scattering information (see [13] for complete reviews of these approaches). Derived from the coarse resolution BRDF products and using land cover stratification and/or vegetation indices, these archetypes are used to retrieve the BRDF shape at high spatial resolution. Finally, Gao et al. [14] retrieved per-crop BRDF parameters using 500 m MODIS BRDF products and a detailed crop type map for the United States. They built an 8-day 1-degree Look-Up Map (LUM) to account for spatial and temporal differences of the per-crop BRDF parameters.

Quantifying the uncertainty related to these BRDF-adjustment techniques is a critical part of producing consistent SR data set [15]. Several approaches to perform such assessment have been explored. The most traditional validation approach relies on albedo validation using in situ measurement [3,11,12,16]. Albedo is derived from the angular and spectral integration of the BRDF among the full sun-view angle range (i.e., white sky albedo) and the solar spectral domain. The albedo uncertainty includes these uncertainties related to various integrations: (i) Angular integration of BRDF; (ii) spectral integration of the spectral albedo [17]; and (iii) spatial integration of the measurement footprint [16]. These integrations make the BRDF uncertainty retrieval more difficult to estimate. Moreover, for any in situ approaches, the method is limited by the number of sites and land cover. Another approach consists of analysis of the time series noise of a BRDF-adjusted SR. In that way, it corresponds more to an evaluation approach rather than a validation approach [15]. The directional effect contributes significantly to the time series noise as the target reflectance varies with the observation geometry [1]. Bréon et al. [18] and Vermote et al. [6] have demonstrated the relevance of this evaluation approach on vegetated surfaces related to smooth surface changes. Nonetheless, it requires temporally dense time series associated with high observation geometry variations, such as MODIS.

We showed above that high spatial resolution sensors usually have a narrow scan angle associated with less temporally dense time series. However, some sensors, such as the ones on board the SPOT Constellation can be pointed as much as 27° on either side of nadir. This pointing capability allows the sensors to acquire imagery of the same target with a relative high frequency and under various view angles. In the early 1990s, Pinter et al. [19] showed already that the analysis of SPOT1 and SPOT2 (launched in 1986 and 1990, respectively) time series requires extensive knowledge of the surface BRDF. Launched in 1998, SPOT4 was located during more than 15 years on a 26-day cycle orbit, but from February to June 2013 before the SPOT4 decommissioning date, the platform was put in a five-day orbit. This six-month experiment was called SPOT4 (take5) [20]. Thanks to the SPOT4 pointing capability, we requested two sets of five-day image time series with two main observation geometries of the same site near Maricopa (Arizona, USA). The acquisitions from the two angles are delayed by only one day. This configuration offers a unique data set to evaluate BRDF-adjustment techniques through the quantification of the high spatial resolution time series noise related to the surface anisotropy. Two other sites in France have overlapped areas that experience a similar multi-angular acquisition configuration.

The objective of this paper is to quantify the uncertainty of five selected high spatial resolution BRDF-adjustment techniques using the time series noise approach. The first four techniques rely on the VJB method and one on the MCD43 MODIS BRDF products. The first technique assumes no surface BRDF variation through time and space. The second technique corresponds to the average model [18] that considers the variation of the BRDF coefficient related to the NDVI variation with no variability in space. Then, we evaluate two methods based on MODIS-derived BRDF VJB parameters: one based on a NDVI-based disaggregation technique, one based on a Land-cover-based disaggregation technique [12]. The last technique, based on the MCD43 products, corresponds to the Look-Up Map [14].

2. The VJB BRDF Method

Looking for an improvement in the albedo temporal resolution that mitigated the assumption of a stable target, Vermote et al. [6] proposed the VJB method that assumes that the BRDF shape variations throughout a year are limited and linked to the Normalized Difference Vegetation Index (NDVI) variations. The method was developed based on the global MODIS data aggregated at 0.05° spatial resolution in order to remove the spatial variability observed at MODIS full spatial resolution [21]. Some studies have explored the potential and the consistency of this method to model the vegetated surfaces’ anisotropy. First, Breon and Vermote [18] and Franch et al. [22] have shown that VJB model leads to almost no difference when compared to MOD43 normalized surface reflectance and albedo products, respectively. Claverie et al. [23,24] have shown the reliability of the model to adjust the BRDF when compared to Landsat and Formosat-2 SR data. Finally, de Abelleyra and Verón [25] have shown the potential of VJB to reduce time series noise of daily MODIS 250 m resolution SR data over a cropped area.

The VJB method is based on the state-of-the-art BRDF model where the surface reflectance (ρ) is decomposed into three kernels (isotropic, volumetric, and geometric) represented by three coefficients (k₀, k₁, k₂):

$$\text{ρ}\left(\theta ,\lambda \right)=k_0\left(\lambda \right)+k_1\left(\lambda \right)F_1\left(\theta \right)+k_2\left(\lambda \right)F_2\left(\theta \right)$$

(1)

where λ is the spectral band and θ defines the geometric configuration of the acquisition described by three angles: the sun zenith angle, the view zenith angle and the relative azimuth angle (difference between view and sun azimuth angle). F₁ is the volume scattering kernel, based on the Ross-thick function derived by [26] and F₂ is the geometric kernel, based on the Li-sparse model [27] but considering the reciprocal form given by Lucht et al. [28]. Vermote et al. [6] introduced two proxy variables, V (as Volumetric) and R (as Roughness), to describe the BRDF “shape”. V and R are defined as k₁/k₀ and k₂/k₀, respectively. Equation (1) can thus be written as:

$$\rho \left(\theta ,\lambda \right)=k_0\left(1+V\left(\lambda \right)F_1\left(\theta \right)+R\left(\lambda \right)F_2\left(\theta \right)\right)$$

(2)

Using the MODIS CMG SR dataset, Vermote et al. [6] found that the V and R variations with time were well correlated to the NDVI variations over vegetated areas. Thus, they generated two linear regressions (Equations (4) and (5), including four coefficients: V₀, V₁, R₀ and R₁) to estimate V and R as a function of the NDVI.

$$V\left(day,\text{λ}\right)=V_0\left(\text{λ}\right)+V_1\left(\text{λ}\right)\times NDVI\left(day\right)$$

(3)

$$R\left(day,\text{λ}\right)=R_0\left(\text{λ}\right)+R_1\left(\text{λ}\right)\times NDVI\left(day\right)$$

(4)

A correction of the directional effect consists of putting the measurement in standard observation geometry, as follow:

$$\rho \left({\theta }_{\text{N}},\lambda \right)=\rho \left({\theta }_i,\lambda \right)\frac{1+(V_0\left(\text{λ}\right)+V_1\left(\text{λ}\right){NDVI}_i)\times F_1\left({\theta }_i\right)+(R_0\left(\text{λ}\right)+R_1\left(\text{λ}\right){NDVI}_i)\times F_2\left({\theta }_i\right)}{1+(V_0\left(\text{λ}\right)+V_1\left(\text{λ}\right){NDVI}_i)\times F_1\left({\theta }_{\text{N}}\right)+(R_0\left(\text{λ}\right)+R_1\left(\text{λ}\right){NDVI}_i)\times F_2\left({\theta }_{\text{N}}\right)}$$

(5)

where the normalized reflectance (ρ(θ_N)) is related to the actual reflectance (ρ(θ_i)). The standard observation geometry (θ_N) is defined for the Sun at 45° from zenith and the observation at nadir.

3. Data

This study is based on the SPOT4 (Take5) experiment accurately described in [20]. Among the 42 SPOT4 (Take5) sites, we focused on three sites that include multi-angular acquisitions. SPOT4 data has a spatial resolution of 20 m in 4 wavelengths (555 nm, 660 nm, 850 nm, and 1600 nm). We limited this study to the Red (660 nm) and NIR (850 nm) bands (Figure 1).

Figure 1. Maps of the study areas, including the SPOT4 satellite tracks (cyan lines, daytime only) during SPOT4 (Take5) experiment. Overlap areas for ProvLanguedoc and Sudmipy sites are shown with yellow polygons; the overlap area for Maricopa site corresponds to the entire footprint. Images are displayed using a false color composite (red = NIR, green = Red, blue = Green). Black boxes in A, B and C correspond to the spatial extend of the referred data frames.

Figure 1. Maps of the study areas, including the SPOT4 satellite tracks (cyan lines, daytime only) during SPOT4 (Take5) experiment. Overlap areas for ProvLanguedoc and Sudmipy sites are shown with yellow polygons; the overlap area for Maricopa site corresponds to the entire footprint. Images are displayed using a false color composite (red = NIR, green = Red, blue = Green). Black boxes in A, B and C correspond to the spatial extend of the referred data frames.

SPOT4 data acquired over the Maricopa (Arizona, USA) site follow a requested multi-angular acquisition plan. Data of the same footprint were indeed acquired from two successive satellite tracks, inducing two distinct view angles (Table 1 and Figure 2). The two acquisitions are delayed by one day and follow the five-day repeat cycle of the experiment. Figure 1C,D shows the study area footprint with the SPOT4 satellite tracks. The study area, defined by a 120 × 60 km² footprint SPOT4, is mostly covered by semi-desert area. About 20% of the area is covered by irrigated crops.

Table 1. View zenith and azimuth angles of the three sites from East and West satellite track.

Sites	View Angles From East Satellite Track		View Angles from West Satellite Track
	Zenith Range	Main Azimuth	Zenith Range	Main Azimuth
Maricopa	(1.4, 8.7)	−105.9	(21.9, 28.0)	73.9
ProvLanguedoc (overlap area)	(24.7, 25.3)	−109.6	(0.2, 0.9)	70.4
Sudmipy (overlap area)	(12.2,12.6)	−109.7	(12.8, 14.2)	70.0

Table 1. View zenith and azimuth angles of the three sites from East and West satellite track.

Sites	View Angles From East Satellite Track		View Angles from West Satellite Track
	Zenith Range	Main Azimuth	Zenith Range	Main Azimuth
Maricopa	(1.4, 8.7)	−105.9	(21.9, 28.0)	73.9
ProvLanguedoc (overlap area)	(24.7, 25.3)	−109.6	(0.2, 0.9)	70.4
Sudmipy (overlap area)	(12.2,12.6)	−109.7	(12.8, 14.2)	70.0

Figure 2. Sun and view angles of the SPOT4 (Take5) data.

Figure 2. Sun and view angles of the SPOT4 (Take5) data.

SPOT4 data acquired over the ProvLanguedoc and Sudmipy sites (both in France) are covered by two parallel footprints of 180 km and 280 km long, respectively, and acquired from two successive satellite tracks (Figure 1B,E,F). The two footprints overlap over an area 10–20 km wide, which was used to analyze the SR differences due to the view angle variation. Since one of the five BRDF techniques (LC-dis, see Section 4.3.2) is based on surrounding land cover and MODIS BRDF a-priori information, we defined an area covering more than 50 times the overlap area to account for a sufficient amount of MODIS data and land-cover types.

MODIS data, used for this study, are derived from the daily M[O/Y]D09GA products. It includes Red (620–670 nm) and NIR (841–876 nm) spectral bands gridded in a 250 m spatial resolution grid in sinusoidal projection. The 250 m resolution is nominal since it is valid only for nadir view. However, with a scan angle of 50°, and accounting for the point spread function, the resolution reaches 800 m for extreme pixels of the scan [29]. Bréon et al. [30] have shown that the multitemporal MODIS measurement at nominal spatial resolution (0.5 km) is affected by a change in the effective resolution for off-nadir observation and by inaccurate registration (about 1/5th of a pixel [31]). They demonstrate that the noise in the apparent reflectance decreases as the spatial resolution coarsens. Consequently, M[O/Y]D09GA products required preprocessing to derive aggregated products at a coarser spatial resolutions. Bréon et al. [30] concluded that the optimal aggregated spatial resolution for MODIS time series analysis is 2000 m. By simply extrapolating this finding to the 250 m nominal resolution, 1000 m resolution may be optimal. We defined a 1250 m resolution ensuring a relatively stable footprint of the measurement even in case of off-nadir observations. The MODIS preprocessing and BRDF derivation is described in Section 4.1.

For the implementation of the LUM technique (Section 4.3.5), we used the 2013 Cropland Data Layer (CDL) available from the US Department of Agriculture (USDA) National Agricultural Statistics Service (NASS) [32]. The CDL is an annual crop-type mask derived from the Landsat data. It accounts for major crops growth in US, including those cropped in the Maricopa area: alfalfa, cotton, barley, and corn. The original CDL, which is in the Albers Conical Equal Area projection at a 30-m resolution, was reprojected to the SPOT4 grid, i.e., 20 m UTM zone 13N.

4. Methodology

4.1. Coarse Resolution Processing and BRDF Derivation

The BRDF of the coarse resolution was derived using the VJB method and the same inversion approach as Vermote et al. [6] used to adjust the BRDF of MODIS CMG data. Instead of aggregating MODIS full resolution to a 0.05 degree grid (~5 km), we defined a finer grid at 1250 m. M[O/Y]D09GA products were averaged over the 1250 m grid, considering only the best quality pixels according to M[O/Y]D09GA quality assessment field. To be consistent with the SPOT4 (Take5) dataset, we used the same SPOT4 metric projections (UTM zone 13N for Maricopa site and Lambert 93 for both French sites).

The V and R BRDF coefficients inversion approach, presented by [6] and based on MODIS CMG data, was applied with the 1250 m-aggregated MODIS data. Although the SPOT4 (Take5) experiment lasted from February to June 2013, we considered the data acquired over the full 2013 year to run the inversion approach in order to cover the full phenology. It resulted in a spatially explicit 1250 m resolution map of slopes and intercepts of V and R coefficients as presented in Equations (4) and (5).

4.2. Land Cover Classification

Unsupervised land cover classifications were derived over the 3 sites using the non-BRDF-adjusted SPOT4 (Take5) data sets. The classifications were used to run the Land-cover-based disaggregation technique (LC-dis, Section 4.3.2) and to stratify the results. To avoid the cloud contamination, we created monthly NDVI composites using the average of all cloud-free and shadow-free pixels. We then gap-filled March to May monthly composites using a temporal linear interpolation technique. Remaining gaps (mostly for February and June data) were extrapolated using nearest temporal data. Finally we applied an unsupervised K-means clustering algorithm [33] on these monthly gap-filled NDVI data sets. A fixed number of 10 classes was defined.

4.3. High Resolution BRDF-Adjustment Techniques

Five high resolution BRDF-adjustment techniques were applied. A summary of the main characteristics is given in the Table 2 ([2,6,12,14,18,32]). Four of them rely on the VJB model but they differ by the use of the MODIS BRDF a-priori information. The average model and the Land-cover-based disaggregation technique were already published by Bréon et al. [18] and Franch et al. [12], respectively. The constant model and the NDVI-based disaggregation techniques were developed for this paper. The fifth method relies on the LUM approaches (introduced by Gao et al. [14]) that used cropland data layer and the MCD43 products as BRDF a-priori information.

Table 2. High resolution BRDF-adjustment techniques summary.

Techniques (short Name)	Reference	BRDF Approach	BRDF Dynamic	Use of Medium Resolution NDVI for Retrieving the BRDF	Use of Moderate Resolution BRDF	Use of Land-Cover Classification
Cst	-	VJB [6]	No BRDF dynamic	No	No	No
Av	[18]		Temporal variation with NDVI	Yes
VI-dis	-		Spatial and temporal variations		MODIS VJB coefficients (aggregated at 1250 m)
LC-dis	[12]			No		Yes (unsupervised classification)
LUM	[14]	MCD43 [2]			MODIS (MCD43 at 500 m)	Yes (CDL [32])

Table 2. High resolution BRDF-adjustment techniques summary.

Techniques (short Name)	Reference	BRDF Approach	BRDF Dynamic	Use of Medium Resolution NDVI for Retrieving the BRDF	Use of Moderate Resolution BRDF	Use of Land-Cover Classification
Cst	-	VJB [6]	No BRDF dynamic	No	No	No
Av	[18]		Temporal variation with NDVI	Yes
VI-dis	-		Spatial and temporal variations		MODIS VJB coefficients (aggregated at 1250 m)
LC-dis	[12]			No		Yes (unsupervised classification)
LUM	[14]	MCD43 [2]			MODIS (MCD43 at 500 m)	Yes (CDL [32])

4.3.1. The Average Technique (Av)

The average model was introduced by Bréon et al. [18]. They derived a single set of slopes and intercepts of V and R derived from simulations carried out on globally-distributed pixels. V and R are consequently derived from these parameters and the local NDVI, i.e., at high spatial resolution NDVI.

4.3.2. The Land-Cover-Based Disaggregation Technique (LC-dis)

The Land-cover-based disaggregation technique was developed by Franch et al. [12]. It consists of using the land cover unsupervised classification to disaggregate the MODIS BRDF derived with the VJB method. With this technique, they can estimate the BRDF parameters at high resolution for each class considered. Notice that they use the whole image to disaggregate the BRDF parameters at high spatial resolution, requiring enough variance to lead to a successful inversion. Therefore, following their recommendations, we filtered SPOT4 images with cloud-free pixel percentage higher than 70%.

4.3.3. The Constant Technique (Cst)

The constant model was developed by our colleagues from the Centre d’Etudes Spatiales de la BIOsphère using SPOT4 data acquired at the Maricopa site. It relies on the VJB model and assumes a single BRDF for all surfaces, described by a constant set of V and R coefficients, derived as followed: Let k be the ratio between SR of two consecutive east and west SPOT4 acquisitions of the Maricopa sites (Equation (6)). A unique value of k was retrieved for each pair of cloud-free data, using the mean SR of the entire footprint. Introducing V and R coefficients from Equation (2) and assuming no change of the isotropic kernel between two consecutive days (${{\text{ρ}}_0}_{{day}_i}={{\text{ρ}}_0}_{{day}_{i+1}}$), we can derive a linear relationship between the two coefficients (Equations (6) and (7)) for each pair. The nine relationships retrieved based on the data pairs acquired in Maricopa sites are reproduced in Figure 3. For both bands, there are recurrent crossing points (except one outlier) from which we defined the V and R constant coefficients values, as displayed in Figure 3.

$${\text{ρ}}_{{day}_i}=k{\text{ρ}}_{{day}_{i+1}}$$

(6)

$$R=\frac{{F_1}_{{day}_i}-k{F_1}_{{day}_{i+1}}}{k{F_2}_{{day}_{i+1}}-{F_2}_{{day}_i}}\times V+\frac{1-k}{k{F_2}_{{day}_{i+1}}-{F_2}_{{day}_i}}$$

(7)

Figure 3. V and R VJB (Vermote-Justice-Bréon) coefficients retrieval for the constant technique. Each line refers to the relationship as defined in Equation (7) for one specific day_i (color bar in the right). (a) red band; (b) NIR band.

Figure 3. V and R VJB (Vermote-Justice-Bréon) coefficients retrieval for the constant technique. Each line refers to the relationship as defined in Equation (7) for one specific day_i (color bar in the right). (a) red band; (b) NIR band.

4.3.4. The NDVI-Based Disaggregation Technique (VI-dis)

The NDVI-based model relies on the VJB coarse resolution BRDF coefficients (see Section 4.1) used as a-priori information of the surface. The model assumed a uniformity of the BRDF behavior across NDVI variations (i.e., the V and R relationship with NDVI does not vary) within the coarse resolution pixel. The coefficients V and R are consequently derived from the V/R slope/intercept of the corresponding coarse resolution pixel and the high resolution SPOT4 NDVI.

4.3.5. The Look-Up-Maps Technique (LUM)

The LUM technique, developed by Gao et al. [14], is based on the crop type-based BRDF information derived from MODIS products at 500 m (MCD43A1). The Crop Data Layer (CDL) is used to identify “pure” 500 m MODIS pixels in a 1 × 1 (latitude and longitude) cell. To generate the LUMs for each crop type, BRDF parameters were collected from the “pure” MODIS pixels on an eight-day basis while characterizing the data into five NDVI intervals (0–0.2, 0.2–0.4, 0.4–0.8, 0.8–1.0). This temporal and NDVI range characterization allows to account for phenological changes. The LUM is then applied for each SPOT4 pixel, accordingly to the CDL class, by using an actual-to-normalized ratio derived from the MODIS BRDF parameters of the “pure” MODIS pixels and the SPOT4 sun-view geometries.

4.4. The Evaluation Method

To quantify the impact of the BRDF-adjustment techniques, we defined an evaluation metric based on the estimate of the Surface Reflectance (SR) time series (TS) noise. This estimate assumes no variation of the surface between two consecutive days (corresponding to acquisitions from east and west satellite tracks). Noise is computed using the root-mean-square deviation (RMSD) between the two measured SR of day_i and day_i + 1 (Equation (8)), and considering one of the five BRDF-adjustment techniques SR TS (Cst, Av, VI-dis, LC-dis, LUM) and non-BRDF-adjusted SR TS (No-Adj). Except for No-Adj, we considered normalized SR, i.e., nadir view angle and 45° sun zenith angle.

$$Noise=\sqrt{{\displaystyle \sum }_{i=1}^n\frac{{\left(\rho \left(da{y}_i\right)-\rho \left(da{y}_i+1\right)\right)}^2}{n}}$$

(8)

The noise ratio reduction (NoiseRatio, Equation (9)) is defined to evaluate the Noise reduction induced by the BRDF-adjustment techniques. It corresponds to the ratio between the Noise of the normalized SR from each technique and the noise of the non-BRDF-adjusted SR (No-Adj). A value higher than 1 indicates a noise increase, i.e., the BRDF-adjusted SR TS is noisier than the non-BRDF-adjusted SR TS.

$$NoiseRatio\left(technique{s}_j\right)=\frac{Noise\left(technique{s}_j\right)}{Noise\left(No Adj.\right)}$$

(9)

Noise and NoiseRatio are consequently computed per pixel for the entire time series. Images, distributions and distribution statistics of Noise and NoiseRatio are used to evaluate the results.

5. Results

5.1. Maricopa Site

Figure 4 shows Surface Reflectance (SR) time series (TS) averaged per class (Section 4.2). Non-BRDF-adjusted (No-Adj) data correspond to the original SR TS while the five BRDF-adjustment techniques correspond to normalized SR TS (nadir view and 45° sun zenith angle). One can notice the high noise occurring with No-Adj SR TS for almost all surface types, enhancing the surface anisotropy. Considering distinct measurements for each of the two view geometries, the SR TS appears smooth indicating the quality of the absolute calibration and the atmospheric correction [34,35]. In most cases, the BRDF-adjustment techniques reduce the noise significantly by reducing the discrepancy between the east and west SR TS. In general, all techniques agreed in finding a similar range of normalized SR.

Even if LC-dis technique displays low level of noise, some high discrepancies occur with this technique. For class #1, corresponding to dark low vegetated shrubland, the LC-dis technique yielded significant differences mostly at the beginning of the experiment. This class, which includes shaded shrubland in mountainous areas, dark bare soil and shrubland surrounding agricultural fields, is weakly homogenous and spread all over the footprint. The retrieved BRDF was consequently not representative of all pixels belonging to this class. The introduction of additional classes would have been necessary to improve the BRDF retrieval of this class. Other issues with LC-dis techniques are detected for some specific days over class #7 (~alfalfa). For instance, the SR TS of class #7 displays a strong discontinuity on 31 March. On this day, bright surfaces were flagged as saturated on SPOT4 data and not used in the processing. This saturation did not concern class #7 pixels exclusively, but it affects most of the coarse resolution pixels (encompassing these class #7 pixels) that were not used for the BRDF inversion (see Section 4.3.2). It yielded to a lower confidence of the retrieved BRDF because it relies on a lower number of used coarse resolution pixels to do the inversion.

Figure 4. Red and NIR per-class Maricopa time series (TS) of the Surface Reflectance (SR) without BRDF adjustment (No-Adj) and time series of the normalized Surface Reflectance (SR, nadir view and 45° sun zenith angle) retrieved from one of the five BRDF-adjustment techniques: Constant model (Cst), average model (Av), NDVI-based disaggregation model (VI-dis), Land-cover-based disaggregation model (LC-dis), LUM technique (LUM). Two consecutive cloud-free measurements with one-day lag (i.e., acquired from two satellite tracks) are connected with a line to highlight the deviation related to the anisotropy effect. Interpreted land-cover class and computed Noises of the displayed TS are reported in the right panels.

Figure 4. Red and NIR per-class Maricopa time series (TS) of the Surface Reflectance (SR) without BRDF adjustment (No-Adj) and time series of the normalized Surface Reflectance (SR, nadir view and 45° sun zenith angle) retrieved from one of the five BRDF-adjustment techniques: Constant model (Cst), average model (Av), NDVI-based disaggregation model (VI-dis), Land-cover-based disaggregation model (LC-dis), LUM technique (LUM). Two consecutive cloud-free measurements with one-day lag (i.e., acquired from two satellite tracks) are connected with a line to highlight the deviation related to the anisotropy effect. Interpreted land-cover class and computed Noises of the displayed TS are reported in the right panels.

One case (class #5-NIR) displays lower noise of non-BRDF-adjusted SR TS than noise of normalized SR TS due to over adjustment of all techniques during vegetation peak (March–April). Compared to all other classes, class #5 (equivalent to barley) during vegetation peak displays an opposite anisotropic behavior: NIR measurement from west track (day_i) is higher than the one from east track (day_i + 1). This specific anisotropic behavior is captured by none of the five BRDF-adjustment techniques. The main hypothesis to explain this anomaly is high SR in the forward scattering configuration not accounted in any of the used BRDF models. Indeed the highest value (first of the two consecutive acquisitions, from west satellite track) corresponds to a 200° relative azimuth angle, a 42° sun zenith angle and a 30° view zenith angle, which makes the observation close to the principal plane in the forward scattering configuration (same view-sun zenith angles, 180° relative azimuth angle) and could be affected by the specular reflection. The forward scattering have been shown to be particularly strong with barley and wheat crop fields [36]. It does not affect class #5 Red band because the SR level is very low, nor the LC since they are poorly vegetated at this time of the year.

Figure 5. Scatterplots of Maricopa surface reflectance (ρ) between consecutive days: Day_i and day_i + 1 corresponding to the western and eastern acquisitions, respectively. The two subplot rows refer to the two spectral bands, the six subplot columns refer to the five BRDF-adjusted SR (see Figure 4 for naming) and the non-BRDF-adjusted SR. (a): No-Adj-Red; (b): Cst-Red; (c): AV-Red; (d): VI-dis-Red; (e): LC-dis-Red; (f): LUM-Red; (g) No-Adj-NIR; (h): Cst-NIR; (i): Av-NIR; (j): VI-dis-NIR; (k): LC-dis-NIR; (l): LUM-NIR.

Figure 5. Scatterplots of Maricopa surface reflectance (ρ) between consecutive days: Day_i and day_i + 1 corresponding to the western and eastern acquisitions, respectively. The two subplot rows refer to the two spectral bands, the six subplot columns refer to the five BRDF-adjusted SR (see Figure 4 for naming) and the non-BRDF-adjusted SR. (a): No-Adj-Red; (b): Cst-Red; (c): AV-Red; (d): VI-dis-Red; (e): LC-dis-Red; (f): LUM-Red; (g) No-Adj-NIR; (h): Cst-NIR; (i): Av-NIR; (j): VI-dis-NIR; (k): LC-dis-NIR; (l): LUM-NIR.

Figure 5 shows the density scatterplots between the SR between two consecutive days (i.e., with opposite view geometry), considering the complete SPOT4 (take-5) data set of Maricopa. Without BRDF-adjustment (subplots a and g), Red and NIR bands display a positive bias, meaning that the western acquisitions generate systematically lower values than the eastern ones. The five BRDF techniques yield to a lower deviation and in most cases (except “Av–Red”) to an over-adjustment of high SR values. The VIS-dis technique displays the lowest error and the best matching with almost no slope in the regression. Fewer data were used for the LC-dis technique. Indeed this technique relies on a per-scene inversion of the BRDF coefficient and we had to filter images with a cloud cover higher than 70% (similarly to [12]).

The Noise (Equation (8)) was calculated per pixel and for the whole time series, considering non-BRDF adjusted and BRDF adjusted TS. Figure 6 shows a spatial subset of the Noise map, centered on an agricultural area of the zone. There is a clear overall decrease of the Noise from No-Adj to any of the five techniques. Most of the noisy fields in the NIR are reduced, with a good performance of the Cst technique, which have been calibrated based on agricultural land cover zones. One can notice some fields circled in subplot (l) for which BRDF-adjusted Noises are higher than the non-BRDF-adjusted Noise. These fields belong to the class #5 (mostly barley) related to the specular reflection, previously mentioned.

Figure 6. Subsets of Maricopa surface reflectance (SR) time series noise. The two subplot rows refer to the two spectral bands, the six subplot columns refer to the non-BRDF-adjusted SR and the five BRDF-adjusted SR (see caption of Figure 4 for naming). Circle in subset-l identifies barley fields with issue in the BRDF-adjustments. (a): No-Adj-Red; (b): Cst-Red; (c): AV-Red; (d): VI-dis-Red; (e): LC-dis-Red; (f): LUM-Red; (g) No-Adj-NIR; (h): Cst-NIR; (i): Av-NIR; (j): VI-dis-NIR; (k): LC-dis-NIR; (l): LUM-NIR.

Figure 6. Subsets of Maricopa surface reflectance (SR) time series noise. The two subplot rows refer to the two spectral bands, the six subplot columns refer to the non-BRDF-adjusted SR and the five BRDF-adjusted SR (see caption of Figure 4 for naming). Circle in subset-l identifies barley fields with issue in the BRDF-adjustments. (a): No-Adj-Red; (b): Cst-Red; (c): AV-Red; (d): VI-dis-Red; (e): LC-dis-Red; (f): LUM-Red; (g) No-Adj-NIR; (h): Cst-NIR; (i): Av-NIR; (j): VI-dis-NIR; (k): LC-dis-NIR; (l): LUM-NIR.

Figure 7. Noise (first column) and NoiseRatio (second column) distributions retrieved from Maricopa normalized SR TS (using five BRDF-adjustment techniques in colored lines, see caption of Figure 4 for naming) and from non-BRDF-adjusted SR TS (No-adj, black line).

Figure 7. Noise (first column) and NoiseRatio (second column) distributions retrieved from Maricopa normalized SR TS (using five BRDF-adjustment techniques in colored lines, see caption of Figure 4 for naming) and from non-BRDF-adjusted SR TS (No-adj, black line).

Considering the entire SPOT4 footprint, Figure 7 shows the Distribution of the SR TS Noise and the NoiseRatio (Equation (9)). The two techniques relying on the VJB MODIS BRDF coefficients (VI-dis and LC-Dis) show the best results. The drop from No-Adj to any of the five techniques is related to the reduction of the deviation between two acquisitions with opposite geometry. The introduction of the NoiseRatio highlights that between 10% and 22% of the pixels that have a value higher than 1 are misadjusted, similarly to the circled field of Figure 6l.

5.2. ProvLanguedoc and SudMipy Sites

The analyses were reproduced over the overlap area of ProvLanguedoc (Figure 8) and SudMipy (Figure 9) sites. LUM technique was not applied since it relies on CDL data available only over the US. Similarly to Maricopa, VI-dis and LC-dis techniques display the lowest level of Noise, similarly to Maricopa site. The Cst technique, calibrated with Maricopa data, performed correctly over these two other sites. One can notice that Cst distribution is the highest Noise values among the five techniques. The Av techniques display narrow distribution, signifying a very good stability of the technique. Considering the two French sites, only 5% of the NoiseRatios are higher than 1, while other techniques display percentage from 6% to 13%. The good performances of these two simple techniques (Cst and Av), which do not include spatial variability of the BRDF shape, suggest that BRDF variability within the range of the view geometry (see Table 1) is limited.

Figure 8. Noise and NoiseRatio distributions retrieved from ProvLanguedoc site. See caption of Figure 7 for details.

Figure 8. Noise and NoiseRatio distributions retrieved from ProvLanguedoc site. See caption of Figure 7 for details.

Figure 9. Noise and NoiseRatio distributions retrieved from SudMipy site. See caption of Figure 7 for details.

Figure 9. Noise and NoiseRatio distributions retrieved from SudMipy site. See caption of Figure 7 for details.

5.3. Overall Results

We extracted the median values from all NoiseRatio distributions, considering distinction between sites and major land covers (LC). Results are presented in Figure 10. In most cases apart from few exceptions, the NoiseRatios are lower than 1, meaning a reduction of the Noise compared to non-BRDF-adjusted SR TS.

Figure 10. Median NoiseRatio values for the five BRDF-adjustment techniques, the two spectral bands, the three sites (Mar = Maricopa, SMP = Sudmipy, PL = ProvLanguedoc), and five majors Land-Cover (LC, Agriculture, Grassland, Shrubland, Forest, Bare soil). Some LC-site combinations were empty and removed. LUM, which was applied only on Maricopa site, was removed from “All sites”. LC were deduced from the LC classification (Section 4.2) based on visual interpretation. Cells are colored using the color bar located in the upper left and Median Noise values of each band-site-LC configuration, corresponding to half-a-row. Lower value of each configuration is written in bold with a star. Last column corresponds to the number of SPOT4 pixels before BRDF-adjustment (in millions).

Figure 10. Median NoiseRatio values for the five BRDF-adjustment techniques, the two spectral bands, the three sites (Mar = Maricopa, SMP = Sudmipy, PL = ProvLanguedoc), and five majors Land-Cover (LC, Agriculture, Grassland, Shrubland, Forest, Bare soil). Some LC-site combinations were empty and removed. LUM, which was applied only on Maricopa site, was removed from “All sites”. LC were deduced from the LC classification (Section 4.2) based on visual interpretation. Cells are colored using the color bar located in the upper left and Median Noise values of each band-site-LC configuration, corresponding to half-a-row. Lower value of each configuration is written in bold with a star. Last column corresponds to the number of SPOT4 pixels before BRDF-adjustment (in millions).

Overall, the VI-dis and LC-dis techniques yield the best performances in most of the LC-site configurations, with an overall NoiseRatio level of 0.5 and 0.6 for Red and NIR bands, respectively (last row of Figure 10). These correspond to a Noise reduction of 50% and 40%, respectively. The Cst and Av techniques perform well specially over grassland LC. Since these two techniques do not include any spatial variability of the BRDF shape, this finding suggests that very simple BRDF-adjustment techniques are good enough to adjust the BRDF over a limited variation of sun-view geometry (see Table 1). LC-dis displays varying performance caused by the fact that the technique is scene-dependent and relies on a matrix inversion that can fail for some classes (see for instance results of class #1 in Figure 4). Applied over Maricopa site only, LUM yields to a similar level of Noise as VI-dis over Agriculture and grassland LC, and slightly higher over shrubland LC. The VI-dis technique appears to be most stable across the three sites, the three selected LC and the two spectral bands.

6. Conclusions

This paper has presented an evaluation of high spatial resolution Bidirectional Reflectance Distribution Function (BRDF) adjustment techniques using multi-angular acquisitions of SPOT4 during the SPOT4 (take5) experiment (February to June 2013). The evaluation was based on the normalized (nadir view, 45° sun zenith angle) SR Time Series (TS) noise, defined as the average deviation between two consecutive observations (one-day lag) from opposite azimuth. Five BRDF-adjustment techniques were evaluated:

• The constant technique (Cst), which considers a uniform BRDF shape for all surfaces;
• The average technique (Av) which supposes the BRDF should shape vary uniformly for all surfaces with NDVI variations;
• The NDVI-based disaggregation technique (VI-dis) based on the disaggregation BRDF coefficient of the Vermote Justice Bréon (VJB) model using NDVI;
• The Land-cover-based disaggregation technique (LC-dis) based on the disaggregation of the MODIS BRDF coefficients of the VJB model using Land-cover; and
• The LUM technique (LUM) based on BRDF coefficients of the MCD43 products and using the US crop data layer.

The SR BRDF-adjustment using any of the five techniques yielded to significant decrease of the SR time series noise. Among the five techniques, the VI-dis one and the LC-dis one produced the lowest level of noise. It reduced the SR TS noise from directional SR to normalized SR by 40% and 50% in average, for red and NIR bands, respectively. However, the LC-dis technique yielded to some contrasted results. Since it is scene-dependent and relies on a matrix inversion, some misadjustment of the BRDF effects occurred on poorly represented land cover. The good performances of the Cst and the Av techniques suggested that a simple approach that does not include spatial variation of the BRDF (Av assumes the BRDF shape variation correlated to the NDVI variation) are good enough to adjust the BRDF in such geometric configurations (the view zenith angle varies from nadir to 25°). Applied only over the Maricopa site, the LUM techniques performed very well over grassland and agriculture land cover. The performances over shrubland were worse because this technique is based on the crop data layer, which is poorly documented outside the crop area. Finally, the VI-dis technique shows one of the best performances with a very good stability across the three sites, the five selected LC and the two spectral bands. This technique, which relies on BRDF a-priori information derived from MODIS, is easy to implement since it is not scene-dependent.

Forestland cover, related to important anisotropy effects, was weakly analyzed in this paper because forest does not cover a large part of the study areas. For this reason, we requested multi-angular acquisitions of a forest site in Finland (Sodankylä) during the SPOT5 (Take5) experiment occurring from April to August 2015. SPOT5 (Take5) data have the same temporal frequency as SPOT4 (take5) (five-day orbit). However, due to the northern latitude, the site is observed from four tracks. Moreover, the low sun zenith angle will accentuate the anisotropy effects making the BRDF-adjustment complex. The resulting observations will offer an ideal data set to evaluate BRDF-adjustment techniques in more complex configurations.