Robust estimation of BRDF model parameters

Remote Sensing of Environment 89 (2004) 63 – 71 www.elsevier.com/locate/rse

Robust estimation of BRDF model parameters Junichi Susaki a,*, Keitarou Hara a, Koji Kajiwara b, Yoshiaki Honda b a

b

Tokyo University of Information Sciences, 1200-2, Yatoh-cho, Wakaba, Chiba 265-8501, Japan Center for Environmental Remote Sensing (CEReS), Chiba University, 1-33 Yayoi, Inage, Chiba 260-8522, Japan Received 3 February 2003; received in revised form 23 September 2003; accepted 5 October 2003

Abstract The effect of the Bidirectional Reflectance Distribution Function (BRDF) is one of the most important factors in correcting the reflectance obtained from remotely sensed data. Estimation of BRDF model parameters can be deteriorated by various factors; contamination of the observations by undetected subresolution clouds or snow patches, inconsistent atmospheric correction in multiangular time series due to uncertainties in the atmospheric parameters, slight variations of the surface condition during a period of observation, for example due to soil moisture changes, diurnal effects on vegetation structure, and geolocation errors [Lucht and Roujean, 2000]. In the present paper, parameter estimation robustness is examined using Bidirectional Reflectance Factor (BRF) data measured for paddy fields in Japan. We compare both the M-estimator and the least median of squares (LMedS) methods for robust parameter estimation to the ordinary least squares method (LSM). In experiments, simulated data that were produced by adding noises to the data measured on the ground surface were used. Experimental results demonstrate that if a robust estimation is sought, the LMedS method can be adopted for the robust estimation of a BRDF model parameter. D 2003 Elsevier Inc. All rights reserved. Keywords: BRDF; Parameter estimation; Robust estimation; LMedS (least median of squares); M-estimators

1. Introduction Spaceborne and airborne remotely sensed data are observed from varying viewing zenith angles, which causes differences in the observed reflectance even though the sensor focuses on the identical target surface. The effects of this variation are corrected by applying the Bidirectional Reflectance Distribution Function (BRDF), which is a function that is based on viewing and illumination geometries. The BRDF effects can be seen in satellite data having a wide swath, e.g. MODIS data have an approximately 2300km swath, which yields a maximum viewing zenith angle collected using a swath of approximately 55j (Strahler, Muller, & MODIS Science Team Members, 1999). The BRDF can be used to correct reflectance to the reflectance at nadir (Kimes et al., 1985; Roujean, Leroy, & Deschamps, 1992). Therefore, the BRDF is one of the most important factors to consider when utilizing satellite data and deriving

* Corresponding author. Tel.: +81-43-236-4633; fax: +81-43-2364633. E-mail address: [email protected] (J. Susaki). 0034-4257/$ - see front matter D 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2003.10.004

biophysical information from the ground surface (Chopping, 2000; Wanner & Strahler, 1995). Various studies related to BRDF have been conducted, including Bidirectional Reflectance Factor (BRF) data measurement (Chopping, 2000; Kimes et al., 1985), retrieval of information from BRDF (Lucht, Shaaf, & Strahler, 2000; Pinty & Verstraete, 1991; Privette, Eck, & Deering, 1997), and BRDF model development (Hu, Wanner, Li, & Strahler, 1997; Lacaze & Roujean, 2001; Roujean et al., 1992; Wanner & Strahler, 1995). Kernel-driven BRDF models are regarded as robust semiempirical BRDF models and are applicable to any type of land cover (Hu et al., 1997). The kernel is a function determined by viewing and illumination geometries. A kernel-driven BRDF model is formulated as follows (Lucht et al., 2000): Rðh; #; /; KÞ ¼ f0 ðKÞ þ f1 ðKÞK1 ðh; #; /Þ þ f2 ðKÞK2 ðh; #; /Þ

ð1Þ

where h and #; are the solar and viewing zenith angle, respectively, / is the view-solar relative zenith angle, K is a waveband of width Dk (k is wavelength), R(h, #, /, K) is

64

J. Susaki et al. / Remote Sensing of Environment 89 (2004) 63–71

the reflectance factor in waveband, fk(K) (k = 0, 1, 2) is the BRDF kernel model parameter in waveband K, f0, f1 and f2 represent the weights for isotropic, volume and geometric scattering, respectively, and Kk(h, #, /) (k = 1, 2) is the BRDF model kernel. Robustness of BRDF model parameter estimation (Gao, Li, Strahler, & Schaaf, 2000; Gao, Schaaf, Strahler, & Lucht, 2001; Li, Gao, & Strahler, 2001; Lucht & Roujean, 2000) is one of the most important aspects when correcting remotely sensed data in order to obtain meaningful physical parameters. Estimation of BRDF model parameters can be deteriorated by various factors; contamination of the observations by undetected subresolution clouds or snow patches, inconsistent atmospheric correction in multiangular time series due to uncertainties in the atmospheric parameters, slight variations of the surface condition during a period of observation, for example due to soil moisture changes, diurnal effects on vegetation structure, and geolocation errors (Lucht & Roujean, 2000). Li et al. (2001) have discussed the importance of a priori knowledge accumulation and its application to linear BRDF model inversion. Linear kernel-driven BRDF models were designed to ease the difficulties involved in inverting nonlinear physical models at the expense of the approximation of the original physics. Normally, unknown parameters of kernel-driven BRDF models are solved by minimizing the sum of the squared fitting error. After a successful inversion of the three parameters from kerneldriven BRDF models, the results are used to obtain the albedos of the pixel, i.e. the white-sky albedo and the black-sky albedo. However, noise and poor sampling reduce the accuracy of the BRDF model inversion, and a failed inversion may result from the combination of either poor sampling or poor directional range, and noisy sample(s), or both. A priori knowledge can indicate when retrieved kernel weights or albedos are outside the expected bounds. This approach is based on Bayes inference theory, which is considered to be the best way to make use of a priori knowledge. A lack of robustness exists when the least squares criterion is implemented due to a small number of large errors in a data set. Parameters of kernel-driven models are generally estimated based on least squares criterion. The least squares error function to be minimized is as follows (Lucht et al., 2000): e2 ðKÞ ¼

1 X ðqi ðhi ; #i ; /i ; KÞ  Ri ðhi ; #i ; /i ; KÞÞ2 d i xi ðhi ; #i ; /i ; KÞ

where qi is a measured reflectance factor, xi (hi, #i, /i, K) is a weight that is assigned to each respective observation (e.g. xi(hi, #i, /i, K) = 1 or xi(hi, #i, /i, K) = qi(hi, #i, /i, K)), and d are the number of degrees of freedom (the number of observations minus the number of parameters fk). Gao et al. (2001) have reported a multikernel least-variance approach

to retrieve and evaluate the albedos from limited bidirectional measurements in order to avoid a biased estimation. In operational application, MODIS BRDF products are produced in the following procedure. If at least seven cloud-free observations of the surface are available during a 16-day period, a full model inversion is attempted. The available data are first evaluated in order to discard outliers, and additional checks are performed in order to assure that the kernel weights are positive. If the data pass these evaluations, a full inversion, or a normal inversion, is performed in order to establish the BRDF parameter weights that provide the best RMSE fit. For those cases with insufficient or poor sampling, or a poor fit, a magnitude inversion is performed, which exploits a priori knowledge (Schaaf et al., 2003). In the present research, the authors investigate the validity of applying the least squares criterion to BRDF model parameter estimation. We focused on the M-estimator and the LMedS method as robust parameter estimators, and the formulation is described herein. Experimental results for the BRDF model parameter estimation obtained by applying the M-estimator, the LMedS method and a conventional least squares method (LSM), respectively, are compared. After a discussion of the experimental results, the conclusions of the present study are presented.

2. Robust estimation Robust estimation techniques have been improved not only in statistics but also in various application fields. Recently, computer vision is one of the major fields that have focused on robust estimation because techniques in computer vision, e.g. orientation and mosaicing of several images, deal, in most cases, with images contaminated by noises (Meer, Mintz, Rosenfeld, & Kim, 1991). If a robust estimation of BRDF model parameters can be achieved, the benefits will be brought to the estimation of BRDF-effectcorrected reflectance and albedo, and to the application of various biosphere models. The T-statistics screening procedure was introduced by Anderson (1958), and was improved for the detection of outliers. Robert (1969) introduced the idea of ‘‘cost’’, and proposed an approach to cost-benefit comparison. Stefansky (1971) introduced the maximum normed residual, defined to be the largest among the absolute values of the normed residuals. Similarly, Tietjen, Moore, and Beckman (1973) proposed a procedure of standardizing the residuals by dividing them by their estimated standard deviations. In the 1970s, M-estimators were proposed, and since then several functions for M-estimators have been tested. M-estimators have been applied for practical data processing, such as image processing (Black, Sapiro, Marimont, & Heeger, 1998; Can, Stewart, Roysam, & Tanenbaum, 2002; Nestares & Heeger, 2000). As similar estimators, W-esti-

J. Susaki et al. / Remote Sensing of Environment 89 (2004) 63–71

mators, R-estimators and L-estimators have been proposed and analyzed. However, they are not so robust in terms of the breakdown point. The breakdown point can be defined as the smallest proportion of the data that can affect an estimate by an arbitrary amount. If the number of data sets is n, the breakdown point of the mean and the median are equal to 1/n and 0.5, respectively. On the other hand, the LMedS method has a high breakdown point, i.e. 0.5. Whereas the LMedS method has a deficiency in that convergence is too slow, computational efficiency can be achieved by applying the MonteCarlo technique. The efficiency and advantages of the LMedS method has been reported in image processing research (Meer et al., 1991; Meer, Stewart, & Tyler, 2000; Rosin, Hervas, & Barredo, 2000). In this chapter, Mestimators and the LMedS method are introduced.

Then, the scale estimator can be estimated from the residuals as rˆ ¼ C  MAD ¼ Cmed fAri  medfri gAg

i

where w(r) is a symmetric, positive-definite function having a unique minimum at zero, and is referred to as the ‘‘robust loss function’’. Eq. (2) is replaced by Eq. (3), in which iteratively reweighted least squares are expressed using a weight function x(r), defined as x(r)=(1/r)((dw(r))/dr). X min xðri Þr 2i ð3Þ i

The weight x(ri) should be recomputed after iteration. The first iteration starts with x(r) = 1, which is equivalent to the ordinary LSM. The LSM, not a robust estimator, can be regarded as one of the M-estimators that has w(r) = r 2 and x(r) = 1. The solution of Eq. (3) is not equivalent with respect to scale. Therefore, the residuals should be standardized by ˆ means of some estimate of the standard deviation r. Then, the residual ri is replaced by the scale-normalized residual ˆ ui = ri/r for the solution. One possibility for the scale estimator is to use the median absolute deviation (MAD) scale estimator, MAD ¼ medfAri  medfri gAg

ð4Þ

ð5Þ

where C = 1.4826 if Gaussian noise is assumed for inliers. Several functions w(ri) have been proposed for M-estimators. Tukey’s biweight (Huber, 1981) is one of the most robust functions and has an advantage in that its influence gradually decreases to zero as a residual becomes large. Therefore, the probability distribution of x(u) in the Tukey’s biweight has trimmed tails. 8 2 c > 2 2 3 > < ½1  ð1  u =c Þ if AuA<¼ c 6 wðuÞ ¼ 2 > > :c if AuA> c 6

2.1. M-estimators M-estimators are robust estimators for parameter estimation and generalize in a straightforward manner for multiparameter problems, even though M-estimators are not automatically scale-invariant and must be supplemented for practical applications by an auxiliary estimate of scale of the typical inlier errors (Huber, 1981). The least squares best-fit method minimizes Sir i2 where ri is the residual of the ith datum, i.e. ri = qi  Ri. The Mestimators reduce the effect of outliers by replacing the quadratic function of the residuals r i2 by another function, yielding X wðri Þ ð2Þ min

65

xðuÞ ¼

8 < ½1  u2 =c2 2

if AuA<¼ c

:

if AuA > c

0

ð6Þ

ð7Þ

The parameter c in Eqs. (6) and (7) is important for the estimation because the parameter value affects the estimation results. However, the optimal value for the parameter c is different depending on the application field (Meer et al., 2000), and there is no agreement with respect to the determination of the parameter c. Compared to M-estimators, in LSM, x(u) = 1 in all u. This means that weights in LSM are equal to all data while weights in M-estimators decrease to zero as residuals become large. However, M-estimators have a low breakdown point, which leads to M-estimators not being robust for the estimation of parameters when more contaminated data are used. 2.2. Least median of squares (LMedS) method Instead of minimizing the sum of squared residuals, Rousseeuw and Leroy (1987) proposed the least median of squares method minimizing their median as follows min med r i2

ð8Þ

The LMedS method cannot be reduced to a least-squares based technique, unlike the M-estimators. The solution of the LMedS must be obtained by a search in the parameter space generated from the data (Meer et al., 1991). Since the space is too large, a randomly chosen subset of data can be analyzed. A Monte Carlo technique is used to select m random subsamples of p different data. In the present research, p is selected as p = 3 because the number of the BRDF parameters to be estimated is three. For each subsample, med r2i is calculated. The subsample that minimizes med r2i can be determined as the solution of Eq. (8).

66

J. Susaki et al. / Remote Sensing of Environment 89 (2004) 63–71

weighted least squares procedure that has high Gaussian efficiency. The robust standard deviation estimate qffiffiffiffi 5 gmed ri2 ð11Þ rˆ ¼ 1:4826f1 þ np can be immediately obtained. The factor 1.4826 is for consistent estimation in the presence of Gaussian noise, and the term 5/(n  p) is recommended as a finite sample ˆ correction (Rousseeuw & Leroy, 1987). Based on r, a weight for each data can be determined, 8 Ari A > > <¼ 2:5 < 1 if rˆ ð12Þ xðri Þ ¼ > > : 0 if Ari A > 2:5 rˆ Fig. 1. Measurement of BRF data at Yohkaichiba in Chiba, Japan. This image was taken on June 7, 2002.

The number of random subsamples m can be determined as follows. Let 1  P be a probability of error and by requiring 1  P1, an equation, P ¼ 1  ½1  ð1  hÞp m

ð9Þ

is obtained where h is the fraction of data contaminated by outliers. Eq. (9) can be rewritten as, m¼

logð1  PÞ log½1  ð1  hÞp

ð10Þ

Then, the parameters can be estimated by solving the weighted least squares problem, X min xðri Þri2 ð13Þ i

Therefore, the estimation of BRDF model parameters by LMedS can be replaced by solving Eq. (13) for the subsamples of data selected by a Monte Carlo technique. This procedure improves the computational efficiency of the LMedS method.

3. Experiments

For example, if p = 3, P = 0.99, and h = 0.3, then m = 11 for any size of data set. Even if p = 3, P = 0.99, and h = 0.5, then m = 35. When Gaussian noise is present in inliers in addition to outliers, the relative efficiency of the LMedS method is low. Rousseeuw proposed combining the LMedS method with a

BRF data measurements were conducted for paddy fields at Yohkaichiba in Chiba, Japan, August 9, 2002. The measurement site is located at 140j33V57E, 35j43V45N, and is a flat area having a homogeneity of more than 2 km2. The measurement area is shown in Fig. 1. For BRF data measurement, we used a car to which a controllable arm was

Fig. 2. Geometry at BRF data measurement. Measurement points are expressed as crosses. Measurements at nadir were conducted four times, that is, every plane contains a nadir measurement. (Solar zenith angle was approximately 20j.)

Fig. 3. Estimated BRDFs for near-infrared band (841 – 876 nm) original data in the principal plane (relative azimuth angle = 0j). In the horizontal axis, positive values indicate backward sensor zenith angles. The data used for the estimation were measured during the period from 12:17 to 12:35 a.m. JST, August 9, 2002. The results of the LSM and the M-estimator in ‘‘(i) added noise =  0.5’’ are indistinguishable in the plots.

J. Susaki et al. / Remote Sensing of Environment 89 (2004) 63–71 Table 1 Reflectances estimated at nadir, + 60j and  60j LSM M-estimator LMedS

all subsamples 4 subsamples 9 subsamples 21 subsamples 35 subsamples

Nadir

+ 60

 60

0.429 0.429 0.430 0.422 0.442 0.434 0.432

0.548 0.547 0.608 0.585 0.523 0.605 0.616

0.431 0.430 0.442 0.442 0.404 0.446 0.442

The BRDF model parameters were estimated using the least squares method (LSM), an M-estimator and the least median squares (LMedS) method. The data used for the estimation were the data in the principal plane. The results in all subsamples of the LMedS method represent the results obtained by searching all possible subsamples from the original data set. Other results for the LMedS method represent the results obtained by searching a limited number of subsamples selected by Eq. (10). ‘‘4 subsamples’’ was determined for the condition that p = 3, P = 0.99, and h = 0.0833. Correspondingly, ‘‘9 subsamples’’, ‘‘21 subsamples’’ and ‘‘35 subsamples’’ were determined for the condition that h = 0.250, h = 0.417 and h = 0.500, respectively. h = 0.0833, h = 0.250 and h = 0.417 are equivalent to cases (II), (III) and (IV) in Fig. 2, respectively.

attached. A spectroradiometer was attached to the rotatable platform at the edge of the arm. The height of the spectroradiometer was approximately 4.0 m from the ground, and the instantaneous field of view (IFOV) of the spectroradiometer was approximately 21j. The viewing azimuth and zenith angles of the platform were controlled via software run on a PC from the ground. The spectroradiometer can measure both sample and reference simultaneously. Thus, both data were used to calculate the reflectance of the paddy fields. The paddy fields to both the left and right sides of the car were measured, and the collected data were combined as a

67

set of BRF data for the analysis. Among several data sets, this experiment used two data sets measured during the period from 12:17 to 12:35 a.m. Japanese Standard Time (JST). In Fig. 2, measurement points are expressed as crosses. In four different relative azimuth angle planes, 0j, 45j, 90j and 135j, sensor zenith angle was set as  60j to 60j at intervals of 10j. Positive sensor zenith angles indicate backward sensor zenith angles. Therefore, reflectances at nadir were measured four times, and the mean reflectance was used as the reflectance at nadir. Reflectances were calculated using both sample digital count and reference digital count. Because the solar zenith angle was 20.4j to 21.4j during the period from 12:17 to 12:35 a.m. JST, data at points having 0j relative azimuth and 20j viewing zenith angles were contaminated by the shadow of the arm. Therefore, the data was excluded in the present experiments. The height of the top of the paddy was approximately 90 to 95 cm. The leaf area index (LAI) was determined by cutting the paddy leaves within a 1 m2 area and scanning the leaves using a scanner. Paddy leaves were collected from three locations. The sample LAIs were 5.03, 4.37 and 3.38, and the mean LAI was 4.26. The kernel-driven model reported by Lucht et al. (2000) was adopted for BRDF/albedo standard product, i.e. MOD43, of the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor onboard the Terra satellite, launched on December 1999 (Strahler et al., 1999). This model adopts RossThick and LiSparse kernels as volume and geometric scattering kernels, respectively (Schaaf et al., 2003). The parameters required in the model were set as h/ b = 2.0 and b/r = 1.0, based on MOD43 products. Nearinfrared (841 – 876 nm) wavelengths were selected in accordance with the MODIS band wavelength. The results

Table 2 Errors of nadir-adjusted reflectances Added LSM noise Reflectance (I) No noise –  0.05 (II) One datum contaminated (1/12 = 8.33%)  0.1  0.2 (III) Three data  0.05 contaminated (3/12 = 25.0%)  0.1  0.2 (IV) Five data  0.05 contaminated (5/12 = 41.7%)  0.1  0.2

0.429 0.425 F 4.24  10 3

M-estimator Error

Reflectance

LMedS Error

– 0.429 –  0.971 F 0.989% 0.425 F 4.71  10  3  0.857 F 1.10%

Reflectance 0.430 0.422 F 1.4110  2

Error –  1.73 F 3.28%

0.421 F 8.48  10  3 0.412 F 1.70  10  2 0.417 F 6.66  10  3

 1.94 F 1.98%  3.88 F 3.95%  2.91 F 1.55%

0.427 F 4.09  10  3  0.387 F 0.954% 0.427 F 7.44  10  3  0.638 F 1.73% 0.431 F 7.19  10  3 0.447 F 1.68% 0.429 F 6.89  10  3  0.274 F 1.60% 0.418 F 8.80  10  3  2.46 F 2.05% 0.418 F 1.64  10  2  2.75 F 3.82%

0.404 F 1.33  10  2 0.379 F 2.66  10  2 0.408 F 7.57  10  3

 5.83 F 3.10%  11.7 F 6.20%  4.86 F 1.76%

0.410 F 1.63  10  2 0.407 F 3.08  10  2 0.409 F 9.11  10  3

 4.41 F 3.81%  5.13 F 7.19%  4.67 F 2.12%

0.416 F 2.17  10 2 0.424 F 1.11  10  2 0.409 F 2.19  10  2

 3.17 F 5.04%  1.29 F 2.57%  4.87 F 5.09%

0.387 F 1.51  10  2 0.346 F 3.03  10  2

 9.71 F 3.53%  19.4 F 7.05%

0.389 F 1.83  10  2 0.350 F 3.75  10  2

 9.22 F 4.27%  18.3 F 8.75%

0.392 F 4.28  10  2 0.401 F 5.60  10  2

 8.78 F 9.95%  6.78 F 13.0%

The BRDF model parameters were estimated using the least squares method (LSM), an M-estimator and the least median squares (LMedS) method. Noise (value =  0.05,  0.1 and  0.2) was added to the original data randomly. For example, in (III), three data were selected among 12 data, and then 12C3 = 220 cases were implemented. The upper and lower values of reflectances are the mean and the standard deviation of the reflectances for all implemented cases. The errors were calculated as the difference between the reflectance of each case and the reflectance of case (I) divided by the reflectance of case (I).

68

J. Susaki et al. / Remote Sensing of Environment 89 (2004) 63–71

estimated from near-infrared band data are reported in the present paper. For M-estimator use, the optimal value for the parameter c in Eqs. (6) and (7) differs according to the application field, e.g. 2.5 or 4.0 (Meer et al., 1991). In the present research, the parameter c in Eq. (7) was set to 4.0. As a result, it was determined that the estimation results were almost the same in the present experiments, even if c was changed from 2.5 to 4.0. The M-estimator was initialized by LSM estimation, and the M-estimator was then determined based on an iteratively reweighted least squares criterion. The parameter C in Eq. (5) was set as C = 1.4826, which is the most appropriate value for Gaussian noise for inliers. The parameter estimation results are shown in Figs. 3 – 6, which plot the measured BRF data, as well as the LSM, Mestimator and LMedS method estimation results. In the horizontal axis, positive values indicate backward sensor zenith angles. The data used to create Fig. 3 were estimated using all of the data measured in the principal plane, excluding the data contaminated by the shadow of the arm, 0j relative azimuth and 20j viewing zenith angles. First, the estimation results for the original data set were compared quantitatively in Fig. 3 and Table 1. The reflectances adjusted at nadir, + 60j and  60j view zenith angles were computed based on the estimated BRDFs. Then, the LMedS method was examined in terms of the robust selection of subsamples from the original data set. As shown in Eq. (10), if p, P, and h are given, then the number of random subsamples m can be determined. In Table 1, the results in all subsamples of the LMedS method represent the results obtained by searching all possible subsamples from the original data set. The other results in the LMedS method represent the results obtained by searching a limited number of subsamples selected by Eq. (10). In each implementation, a different number of subsamples was randomly selected, and the optimal parameters were determined. ‘‘4 subsamples’’ in Table 1 was determined for the condition p = 3, P = 0.99, and p = 0.0833. Correspondingly, ‘‘9 subsamples’’, ‘‘21 subsamples’’ and ‘‘35 subsamples’’ were determined for the condition h = 0.250, h = 0.417 and h = 0.500, respectively. h = 0.0833, h = 0.250 and h = 0.417 are equivalent to cases (II), (III) and (IV) in Fig. 2, respectively. As a result, h = 0.500 can be a stable value, in which the estimated parameters are closest to the parameters obtained by searching all possible subsamples. Therefore, h = 0.500 was adopted for the application of the LMedS method in the following experiments. Finally, we compared the LSM, the M-estimator and the LMedS method by changing the proportion of noise contamination and the noise absolute value. The results are compared in Table 2. The BRDF model parameters were estimated by the LSM, an M-estimator and the LMedS method. Noises were randomly added to the original data, and the proportion of noise addition was varied. We simulated one-data contamination (1/12 = 8.33% contamination),

three-data contamination (3/12 = 25.0% contamination), five-data contamination (5/12 = 41.7% contamination). For example, in the case that three data were selected and

Fig. 4. Examples of estimated BRDFs for near-infrared band (841 – 876 nm) data in the principal plane (relative azimuth angle = 0j). Noise (value =  0.05,  0.1 and  0.2) was randomly added to one original datum. In the horizontal axis, positive values indicate backward sensor zenith angles. The data used for the estimation were measured during the period from 12:17 to 12:35 a.m. JST, August 9, 2002. The results of the LSM and the M-estimator in ‘‘(i) added noise =  0.5’’ are indistinguishable in the plots.

J. Susaki et al. / Remote Sensing of Environment 89 (2004) 63–71

contaminated among 12 data, 12C3 = 220 cases were implemented. The ratio of noise contamination h in Eq. (10) was determined as the contamination ratio, and the number of subsamples was changed for the one-data, three-data and five-data contamination cases. The added noise value was set as  0.05,  0.1 and  0.2 in Table 2, respectively. In each case, the reflectances adjusted to nadir were calculated from the BRDF model parameters estimated by the LSM, an Mestimator and the LMedS method. The errors in the tables were calculated as the difference between the adjusted reflectance and the reflectance estimated from parameters using the original data, divided by the estimated reflectance. The mean and standard deviation of the reflectances were calculated for all of the implemented cases. Figs. 4 –6 show one of the results of BRDF model parameter estimation. For application of the LMedS method, first, Eq. (13) was implemented for the subsamples of the data selected via the Monte Carlo technique to replace med r2i in Eq. (8). However, the reflectance adjusted to + 60j became more than 0.70, which is an overestimated value considering Table 1 and Fig. 3. Therefore, we determined that the implementation of Eq. (13) was not appropriate for the present experiments, and as such med r2i in Eq. (8) was exploited in order to determine the optimal subsample and parameters.

69

the LMedS method are more robust than the nadiradjusted reflectances estimated by either the LSM or M-estimator.

4. Discussion Fig. 3 and Table 1 show the BRDFs estimated for the original data set. The figure and the table show that although the BRDFs estimated by LSM and by M-estimator are quite similar, those estimated by the LMedS method are a little different. Parameter estimation was affected by the noise that was added randomly to the reflectance in the experiments, and the difference between the results for the LSM, an Mestimator and the LMedS method are shown in Figs. 4 – 6. In Table 2, the results for the LMedS method are shown to be affected to a lesser degree by the decrease in the contaminated reflectance compared to the results for the LSM and the M-estimator. Figs. 4 and 5 show that the results for the M-estimator and the LSM are similar in most cases and that the results for the LMedS method are least affected by added noised among three estimators. Measured data are often contaminated by errors caused by imperfect cloud mask and atmospheric correction, and thus become outliers. Theoretically, the LMedS method can mitigate the effects of outlier existence because this method produces estimates based on the median value, which should be a robust statistics for outliers. On the other hand, an M-estimator can be applied for parameter estimation as a robust estimator because the weights for outliers approach zero. The experimental results reveal that while the LMedS method is appropriate for delineating the overall characteristics of the BRDF, an M-estimator is not so robust for outliers. In most of the cases shown in Table 2, nadir-adjusted reflectances estimated by

Fig. 5. Examples of estimated BRDFs for near-infrared band (841 – 876 nm) data in the principal plane (relative azimuth angle = 0j). Noises (value =  0.05,  0.1 and  0.2) was randomly added to three original data. In the horizontal axis, positive values indicate backward sensor zenith angles. The data used for the estimation were measured during the period from 12:17 to 12:35 a.m. JST, August 9, 2002. The results of the LSM and the M-estimator in ‘‘(i) added noise =  0.5’’ are indistinguishable in the plots.

70

J. Susaki et al. / Remote Sensing of Environment 89 (2004) 63–71

(Lacaze & Roujean, 2001). There is a possibility that hot spot data are dealt with as ‘‘contaminated data’’ for the estimation of the parameters. Applying the combination of these kernels sometimes returns negative parameters. Theoretically, the parameters should be non-negative values, but negative values were detected in the present experiments. The present research focuses on BRDF products derived from MODIS data, the kernels of which were exploited. However, the deficiency of BRDF models may lead to the misdiagnosis of discrimination between inliers and outliers. Convergence is also one of the most important aspects of such parameter estimations. In M-estimators, the convergence is not always guaranteed. Initialization of the parameters to be estimated may affect the estimation results and the convergence. In the experiment, parameters were initialized by the LSM. The initialization was based on the assumption that the optimal parameters may be around the parameters obtained by the LSM. In some cases, the assumption does not hold. However, the initialization was considered to be reasonable because there was no a priori information about the optimal parameters. In contrast, the LMedS method guarantees the convergence of the estimation. Indeed, the LMedS method requires much more estimation time than M-estimators and the LSM, because the iteration in the LMedS method was determined by Eq. (10). In addition, subsampling from the original data set is one of important aspects when applying the LMedS method because the subsampling affects the estimation results. If the fraction of data contaminated by outliers is determined as 0.5, then the estimation time will increase accordingly. However, the estimation results reveal that the LMedS method is a much more robust estimator than either an Mestimator or the LSM. Moreover, the LMedS method has high a breakdown point, i.e. 0.5, which is a powerful advantage for robust estimation. Therefore, we recommend the LMedS method as the most appropriate robust estimator of BRDF model parameters. Selection of parameter estimation method is quite important, not only because the estimation results can be affected, but also because the albedos derived from the results can be affected.

5. Conclusions Fig. 6. Examples of estimated BRDFs for near-infrared band (841 – 876 nm) data in the principal plane (relative azimuth angle = 0j). Noises (value =  0.05,  0.1 and  0.2) was randomly added to five original data. In the horizontal axis, positive values indicate backward sensor zenith angles. The data used for the estimation were measured during the period from 12:17 to 12:35 a.m. JST, August 9, 2002. The results of the LSM and the Mestimator in ‘‘(i) added noise =  0.5’’ are indistinguishable in the plots.

We would like to discuss parameter estimation in terms of the BRDF model itself. The RossThick-LiSparse kernels used in the experiment cannot represent real scattering completely. In particular, a hot spot phenomenon is reported to be difficult to express for the combination of these kernels

When BRDF model parameters are estimated, it is reasonable to consider that the data to be used for the estimation may contain noise, e.g. derived from imperfect atmospheric correction. In order to overcome this noise contamination and obtain stable results for parameter estimation, a robust estimation method should be adopted. In the present paper, the experimental results obtained in deriving the BRDF model parameters are described, and the robustness of the estimation was found to be a necessary condition. An Mestimator and the LMedS method were investigated as robust parameter estimators, and the LSM, an M-estimator and the LMedS method were compared to the estimation results for

J. Susaki et al. / Remote Sensing of Environment 89 (2004) 63–71

BRF data to which noises were randomly added. The experimental results demonstrate that the LMedS method has robustness for the estimation of the BRDF model parameter.

References Anderson, T.W. (1958). An introduction to multivariate statistical analysis. New York: Wiley. Black, M. J., Sapiro, G., Marimont, D. H., & Heeger, D. (1998). Robust anisotropic diffusion. IEEE Transactions on Image Processing, 7, 421 – 432. Can, A., Stewart, C. V., Roysam, B., & Tanenbaum, H. L. (2002). A feature-based, robust, hierarchical algorithm for registering pairs of images of the curved human retina. IEEE Transactions on Pattern Analysis and Machine Intelligence, 24, 347 – 364. Chopping, M. J. (2000). Testing a LiSK BRDF model with in situ bidirectional reflectance factor measurements over semiarid grasslands. Remote Sensing of Environment, 74, 287 – 312. Gao, F., Li, X., Strahler, A., & Schaaf, C. (2000). Acquiring a priori knowledge from ground and spaceborne BRDF measurements. IEEE Int. Geosci. Remote Sensing Symp. Proc., vol. 2, pp. 718 – 720. Gao, F., Schaaf, C. B., Strahler, A. H., & Lucht, W. (2001). Using a multikernel least-variance approach to retrieve and evaluate albedo from limited bidirectional measurements. Remote Sensing of Environment, 76, 57 – 66. Hu, B., Wanner, W., Li, X., & Strahler, A. H. (1997). Validation of kerneldriven semiempirical models for global modeling of bidirectional reflectance. Remote Sensing of Environment, 62, 201 – 214. Huber, P. J. (1981). Robust statistics. New York: Wiley. Kimes, D. S., Newcombe, W. W., Tucker, C. J., Zonnefeld, I. W., van Vijingaarden, W., de Leeuw, J., & Epema, G. F. (1985). Directional reflectance factor distributions for cover types of Northern Africa. Remote Sensing of Environment, 18, 1 – 19. Lacaze, R., & Roujean, J. L. (2001). G-function and HOt SpoT (GHOST) reflectance model application to multi-scale airborne POLDER measurements. Remote Sensing of Environment, 76, 67 – 80. Li, X., Gao, F., & Strahler, A. H. (2001). A priori knowledge accumulation and its application to linear BRDF model inversion. Journal of Geophysical Research, 106, 11925 – 11935. Lucht, W., & Roujean, J.L. (2000). Considerations in the parametric modeling of BRDF and albedo from multiangular satellite sensor observations. Remote Sensing Reviews, vol. 18, 343 – 379.

71

Lucht, W., Schaaf, C.B., & Strahler, A.H. (2000). An algorithm for the retrieval of albedo from space using semiempirical BRDF models. IEEE Transactions on Geoscience and Remote Sensing, 38, 977 – 998. Lucht, W., Schaaf, C. B., & Strahler, A. H. (2000). An algorithm for the retrieval of albedo from space using semiempirical BRDF models. IEEE Transactions on Geoscience and Remote Sensing, 38, 977 – 998. Meer, P., Mintz, D., Rosenfeld, A., & Kim, D. Y. (1991). Robust regression methods for computer vision: a review. International Journal of Computer Vision, 6, 59 – 70. Meer, P., Stewart, C. V., & Tyler, D. E. (2000). Robust computer vision: An interdisciplinary challenge. Computer Vision and Image Understanding, 78, 1 – 7. Nestares, O., & Heeger, D. J. (2000). Robust multiresolution alignment of MRI brain volumes. Magnetic Resonance in Medicine, 43, 705 – 715. Pinty, B., & Verstraete, M. M. (1991). Extracting information on surface properties from bidirectional reflectance measurements. Journal of Geophysical Research, 96, 2865 – 2879. Privette, J. L., Eck, T. F., & Deering, D. W. (1997). Estimating spectral albedo and nadir reflectance through inversion of simple BRDF models with AVHRR/MODIS-like data. Journal of Geophysical Research, 102, 29529 – 29542. Rosin, P. L., Hervas, J., & Barredo, J. I. (2000). Remote sensing image thresholding for landslide motion detection. In: 1st Int. Workshop Patt. Recognition Tech. in Remote ( pp. 10 – 17). Roujean, J. L., Leroy, M., & Deschamps, P. Y. (1992). A bidirectional reflectance model of the Earth’s surface for the correction of remote sensing data. Journal of Geophysical Research, 97, 20455 – 20468. Rousseeuw, P. W., & Leroy, A. M. (1987). Robust regression and outlier detection. New York: Wiley. Schaaf, C. B., Gao, F., Strahler, A. H., Lucht, W., Li, X., Tsang, T., Strungnell, N. C., Zhang, X., Jin, Y., Muller, J. P., Lewis, P., Barnsley, M., Hobson, P., Disney, M., Roberts, G., Dunderdale, M., Doll, C., d’Entremont, R. P., Hu, B., Liang, S., Privette, J. L., & Roy, D. (2003). Remote Sensing of Environment, 83, 135 – 148. Stefansky, W. (1971). Rejecting outliers by maximum normed residual. Annals of Mathematical Statistics, 42, 35 – 45. Strahler, A. H., Muller, J. P., & MODIS Science Team Members (1999). ‘‘MODIS BRDF/Albedo product: Algorithm theoretical basis document’’, NASA EOS MODIS Doc., version 5.0 (available at http:// modis.gsfc.nasa.gov/data/atbd/atbdmod09.pdf). Tietjen, G. L., Moore, R. H., & Beckman, R. J. (1973). Testing for a single outlier in simple linear regression. Technometrics, 15, 717 – 721. Wanner, W., & Strahler, A. H. (1995). On the derivation of kernels for kernel-driven models of bidirectional reflectance. Journal of Geophysical Research, 100, 21077 – 21089.