Lumped Parameter Modelling - DOC

Document Sample
Lumped Parameter Modelling - DOC Powered By Docstoc
					Vegetation Science – MSc Remote Sensing UCL                                                        Lewis

Lumped Parameter Modelling
P. Lewis RSU, Dept. Geography, University College London, 26 Bedford Way, London WC1H 0AP,


1. Introduction

The aim of these notes is to introduce concepts, model, and applications of ‘simple’ lumped parameter
models of canopy reflectance and scattering.

In the previous sections of this course, we introduced the radiative transfer (RT) equation as a framework
for the calculation of (optical) reflectance as a function of canopy and soil biophysical variables (leaf
dimensions and density or LAI, leaf and soil moisture content, leaf biochemistry etc). Formulating the RT
equation and solving for a given canopy description allows us to investigate (through analytical or
numerical means) the relationship between these ‘fundamental’ descriptions of canopy state and the
remote sensing signal. It allows us to investigate issues such as the sensitivity of a form of signal (e.g.
Landsat TM waveband reflectance) to canopy properties to allow us to make decisions on the type of data
we might be able to use to ‘solve the remote sensing problem’ – i.e. to derive estimates of biophysical
quantities from a remote sensing signal.

In this section, we will investigate how more generalised/approximate forms of canopy
reflectance/scattering models can be developed and exploited for a range of practical tasks.

2. Linear Models
An important concept in many of the models we will be dealing with is that of the ‘linear model’. We can
define a linear model as follows: For some set of independent variables x ( = {x0, x1, x2, … , xn}), we have
a model of a dependent variable y which can be expressed as a linear combination of the independent

The following are examples of linear models:

y  a 0  a1 x1                                                                            (2.1)

y  a 0  a1 x1  a 2 x 2                                                                  (2.2)

y      a x
        i 0
               i   i                                                                       (2.3)

Equation 2.1 is the form of a ‘standard’ linear regression you will have come across (‘y = mx + c’).
Equation 2.2 is a ‘multi-linear’ form. Equation 2.3 is the general form for n+1 terms. As far as we are
concerned here, one of the major features of a linear model is that we can use matrix inversion to solve
for the model parameters (a ( = {a0, a1, a2, … , an}). We will see how to do this later, but for the moment
can note that:
Vegetation Science – MSc Remote Sensing UCL                                                                           Lewis

y  ax                                                                                                      (2.4)

Many other equations can be expressed in this form, given a suitable transformation, so that the following
are also linear models:

                       i n
y  a 0   a i sin  x i   bi cos x i                                                                 (2.5)
                       i 1

                       i n
y  a 0   a i sin x i  bi                                                                             (2.6)
                       i 1

y        a x                     a 0  a1 x 0  a 2 x 0  ...  a n x 0
                              i                          2               n
                   i     0                                                                                   (2.7)
          i 0

y  a 0 e  a1 x                                                                                             (2.8)

Equation 2.5 is a Fourier series – it has 2 terms per ‘i’, but is otherwise directly in the form of equation
2.4. Clearly, equation 2.6, a magnitude-phase representation of a Fourier series, can be expressed in an
equivalent form to 2.5. Equation 2.7 sometimes confuses people: a polynomial, such as a cubic, is a linear
function. Equation 2.7 requires transformation to a linear form:

ln  y   ln a 0   a1 x                                                                                  (2.9)

since ln(ab) = ln(a)+ln(b).

2.1 Linear Mixture Modelling
As an example, we can consider the ‘linear mixture model’ (or ‘spectral mixture modelling’) often used in
remote sensing1,2. We can consider a remote sensing measure, e.g. of spectral reflectance, r, of a pixel as
a summation or integration of a set of n ‘component’ measures (reflectances), i weigthed by their relative
proportions (‘fractional covers’), Fi, where i 0 Fi  1 . These components will vary according to the
                                               i  n 1

application and the area being studied, but may, for instance, be different mineral soils, or soil and
vegetation. Then, we can write:

r  i 0  i Fi
        i  n 1

we can express this as:

r  F                                                                                                      (2.11)

where  is a vector of n reflectance terms (e.g. in a given waveband), and F is a vector of n fractional
cover terms.

  G. I. Metternicht and A. Fermont, 1998, “Estimating Erosion Surface Features by Linear Mixture Modeling”, Remote
Sensing of Environment, Volume 64, Issue 3, June 1998, Pages 254-265
  Settle, J.J. and Drake, N.A., 1993., “Linear mixing and the estimation of ground cover proportions”. International Journal of
Remote Sensing 14, pp. 1159–1177
Vegetation Science – MSc Remote Sensing UCL                                                          Lewis

This model assumes that only first-order interactions are important in defining r. We know that multiple
scattering can be dominant in some cases – e.g. vegetation canopy reflectance in the near infrared, so in
this case, the component reflectance i would be specified as an equivalent ‘canopy’ reflectance, rather
than using a leaf-level measure. Even so, equation 2.10 misses out multiple interactions between
components – e.g. scattering interactions between the canopy and soil, which are assumed small. We shall
return to further develop this model later.

One of the main uses of equation 2.10 is to attempt to derive estimates of fractional cover from a
multispectral measurement of r, assuming that a set of ‘end-member’ reflectance spectra i are

r     F                                                                                (2.12)

The remote sensing problem in this case is to derive an estimate of F from a set of measurements of .
We can express equation 2.12 for all wavebands considered in vector-matrix form:

r  F                                                                                      (2.13)

where r = {r, r, … rm, 1.0} is a vector of length m+1 spectral reflectance measurements (with a 1.0 at
the end) over m wavebands, F is as above, the proportions (fractional cover) vector of length n, and  is
an nx(m+1) matrix, the columns of which are the end-member reflectance spectra i with 1.0 on the
end of each vector.

 r 0      0 0         0 1        0 2             0 n 1   P0 
                                                                               
 r1       1 0         1 0        1 2             1 n 1   P1 
                                                                        P       (2.14)
                                                                              2 
 rm 1     m 1 0     m 1 1     m 1 21      m 1 n 1      
 1.0      1.0                                                                P 
                            1.0             1.0              1.0             n 1 

                                                                                                 i  n 1
The role of the 1.0 terms on the end of each vector is (one way of) expressing the constraint    i 0
                                                                                                            Fi  1 .
If n=m+1,  is a square matrix, and we can solve for F from:

F  r                                                                                      (2.15)

where  is the inverse of  . So, if we have a measurement of spectral reflectance, e.g. in 2 wavebands
(m=2), and a set of 3 (n=3) ‘end-member’ spectra in these wavebands, we can determine the proportion of
those spectra in the measurement.
Vegetation Science – MSc Remote Sensing UCL                                                       Lewis

                                   Figure 2.1 Linear Mixture Model

Figure 2.1 shows an example mixture model, demonstrating the method graphically for m=2, n=3. End-
members 1-3 are shown. The measured reflectance r is seen to be a combination of its proportionate
distance from 1 to 3 (around 0.5 here) and its proportionate distance from 1 to 2 (around 0.2 here).
So, r = 1+0.5 (3-1) + 0.2 (2-1) = 0.31 + 0.2 2 +0.5 3. The figure also demonstrates well one feature
of the linear mixture model – that the end-member spectra must define vertices of the convex hull of all

Such methods can be used for practical applications, but one must be aware of the following issues:

1. The method, as described, is not robust to error in measurement or end-member spectra;
2. Proportions must be constrained to lie in the interval (0,1) – this is effectively the convex hull
   constraint described above;
3. Only m+1 end-member spectra can be considered in the mixing;
4. The method is dependent on the prior definition of end-member spectra;
5. The method cannot directly take into account any variation in component reflectances (e.g. due to
   topographic effects).

Of these limitations, number 1 is often the most serious, particularly in the presence of noisy data (e.g.
noise introduced by atmospheric correction, on top of any sensor/quantisation noise).

2.2 Linear Mixture Modelling in the presence of Noise
The linear mixture model developed above allowed an estimation of m+1 parameters (proportions) from
m observations (m wavebands). This is because we have m+1 linear simultaneous equations from these m
observations. If there are fewer than m independent observations, the model is under-determined for m+1
parameters, and we cannot determine the parameters uniquely. If, there are more than m independent
observations, we have an over-determined problem (i.e redundancy), which allows us to use statistical
methods to reduce uncertainty in the derived parameters.
Vegetation Science – MSc Remote Sensing UCL                                                                                                          Lewis

The simplest way to do this is using the Method of Least Squares (MLS). To employ this for a linear
model, we consider that our model may have an error associated with it:

r  F  e                                                                                                                                  (2.16)

where eis a vector of residuals (discrepency between model and measurement). The MLS attempts to
minimise the sum of the squares of the error e, i.e. e  e . This is achieved by setting the partial derivatives
of r   F   r   F   l 0
                                              l  m 1
                                                         r       F       2
                                                                                     e  e with respect to each of the model parameters Fi, to
zero. Thus:

 l 0                 r    F 
                                                                                                     0
              l  m 1              2

                                       2 l m1  r    F     F
                                           l 0   
                                                                Fi                                 
                                                                                                    

If the model is linear,     F Fi    i , so:    
                  0                2l 0
                                              l  m 1
                                                         r                  
                                                                     F   i       
       l  m 1
       l 0
                  r   i               
                                                 l  m 1
                                                 l 0
                                                                
                                                                        F   i   
Equation 2.17, for i = 0, 1, …, n-1 gives a set of n simultaneous linear equations with n unknowns (F)
which can be used to solve for F.

We can write this in matrix form, O  M P , as following the pattern:

            rl  l 0                 l 0  l 0                     l 1  l 0              l n1  l 0  F0 
                                                                                                                                   
l  m 1
            rl  l 1      l  m 1
                                          l 0  l 1                     l 1  l 1              l n1  l 1  F1 
                           
                                                                                                                                     (2.18)
    l 0
                            l 0                                                                                                 
            r    
            l       l n 1 
                                              
                                             l 0     l n 1                    l 1  l n1      l n1  l n1  Fn1 
                                                                                                                                    

where O is an ‘observations’ vector, M is the ‘model’ matrix (contains only terms assosciated with the
end-member spectra – the model), and P is the ‘parameter’ vector – the model parameters we wish to
solve for.

This can be solved by matrix inversion as above. Note the ‘patterns’ in the matrix and vectors –
remembering this will allow you to formulate for any similar linear case. The vector on the LHS
comprises a summation over all observations (wavebands here) of the product of the observation and the
model term (end-member spectrum here) associated with that location in the vector. The matrix is solely a
function of model terms (end-member spectra) with a clear pattern. Note that we have not dealt with
applying any constraints here – for the mixture modelling case, we must constrain34 the proportions to lie
in the correct limit (0 to 1) and also constrain their sum to 1 as above.
    This can be achieved using Lagrange Multipliers.
 Lewis, P. (1995), “The utility of kernel-driven BRDF models in global BRDF and albedo studies”, Geoscience and Remote
Sensing Symposium, 1995. IGARSS '95. 'Quantitative Remote Sensing for Science and Applications', International , Volume:
2 , 10-14 Jul 1995 Page(s): 1186 -1188 vol.2
Vegetation Science – MSc Remote Sensing UCL                                                           Lewis

The solution to this is a considerable improvement over equation 2.13 in providing a robust estimate, as it
provides a ‘best fit’ (in the least squares sense) of the model parameters for an over-constrained case.

2.3 Best-fit Line
An example you may have come across before that uses this approach is in providing the ‘best fit’ of a
line to a set of data points. In this case, the model is:

y  c  mx                                                                                   (2.19)

We can use the work we did above to go straight to the solution for this. Following the pattern of equation
2.18, O  M P :

       yl  l n1  1
l  n 1                                    xl  c 
   y x    x
                                              
                                            xl2  m 
 l 0  l l   l 0  l                          

which we may more conveniently write:

 y  1                          x  c 
                                  
 yx   x                         2 
                                 x  m 
                                      

where the bar represents a mean over the observation set, by dividing both sides by n. The inverse matrix
M can be written as:

     1         1  x2  x
                           
M                                                                                           (2.21)
                xx   x 1 
                           

where  xx is the variance in x,  xx  x 2  x . For larger matrices, we will use alternative methods to solve
         2                           2

for the inverse. It is left as an exercise to show that:

           y   2
                 xx    x xy
                                    2
                    xx
                                      2

where  xy is the covariance between x and y.

Using this example, we will examine a number of issues relevant to more general linear modelling. First,
recall from above that the sum of the squares of the errors (residuals) is:

           l  n 1
e2           y  c  mx 
                       i               i
            l 0

We can clearly calculate this once the model parameters (m and c) have been worked out (it can also be
calculated more directly). This is an important term in seeing how appropriate the model (straight line
here) is to the dataset. Sometimes a similar term is used as an estimate of noise in the data, although such
Vegetation Science – MSc Remote Sensing UCL                                                                  Lewis

a measure tends to make the noise estimate too highly variable. Rather than e2, the sum of the squared
errors, we tend to use the normalised measure Root Mean Squared Error (RMSE):

RMSE                                                                                               (2.22)

where n is the number of observations and m here is the number of model parameters (2 here) – (n-m) is
the number of degrees of freedom (DOF) of the system of equations.


                                               x2          x                   x1         x

                                          Figure 2.2 Linear Regression

You may be used to placing ‘error bars’ on graphs for which you have fit a straight line to (figure 2.2). A
more general concept is the the idea of the ‘weight of determination’ which allows you to predict
uncertainty in any of the parameters of a linear model, or a linear combination thereof. This is discussed
in more detail by Lucht and Lewis (2000)5.

We can write the weight of determination, 1 w at some location

                1  c 
yx   Q  P      
                 x   m                                                                          (2.23)
                   


1     1
  Q M Q

where T denotes the transpose operation. Plugging 2.23 into 2.24:

 Lucht, W. and P. Lewis (2000) “Theoretical noise sensitivity of BRDF and albedo retrieval from the EOS-MODIS and
MISR sensors with respect to angular sampling”. International Journal of Remote Sensing 21(1) 81-89.
Vegetation Science – MSc Remote Sensing UCL                                                          Lewis

   1
       xx     

w        xx

The uncertainty (‘mean error to be feared’) associated with a prediction of y(x),  is:

 e                                                                                        (2.26)

where e is the uncertainty associated with the measurements. We typically use the RMSE to approximate
e, as noted above. Lucht and Lewis (2000)5 describe the weight of determination as a ‘noise inflation
factor’, as it relates noise in the data to uncertainty in model prediction (equation 2.26).

We can see from equation 2.25 that the minimum error in the prediction of y from a straight line fit is
obtained at x = x - the error (‘uncertainty’) will increase as we depart from this value as a quadratic in x.
So uncertainty in a model prediction is determined by the noise in the data and the sampling over x – for
example – the sampling of x in figure 2.2 will allow a better prediction of y(x2) (effectively an
interpolation) than it will y(x1) (effectively an extrapolation). We can see from equation 2.26 that the
higher the variance of x, the lower the uncertainty.

3 Lumped Canopy Reflectance/Scattering Models
In the previous lecture, we developed a description of the reflectance or scattering from a canopy which
was driven by what we might term ‘fundamental’ biophysical parameters – these are the terms that we
generally consider to drive physically-based models. There are many applications in Earth surface remote
sensing when we do not need to know these terms – instead a much simpler parameterisation may often
suffice, particularly if an alternative parameterisation might provide a more robust measure than a ‘full’
description of the surface.

Examples of this are: (i) when considering shortwave energy budgets – when generally a measure of
albedo will suffice; (ii) for atmospheric correction; (iii) in tracking changes in the remote sensing signal
to detect changes, a parameterisation of the dynamics of spectral reflectance may suffice; (iv) in erosion
modelling, a measure of canopy cover is a much more important requirement than any details of leaf
angle distribution or other ‘fundamental’ terms; (v) when we have sufficient ground observations to
consider using a calibrated model. Examples (i) to (iii) can be, as we shall see, somewhat related – they
require a method of interpolating and extrapolating an observed signal. Example (iv) is a slightly different
form of problem, requiring some generalisation of canopy description. Example (v) is used widely in
remote sensing to estimate some quantity by from a signal which may have some variation due to
system/satellite effects (e.g. varying viewing/illumination angles over a scene) which need to be
accounted for in an otherwise essentially empirical relationship.

There are several ways of arriving at such models – the two main approaches we might take are:

1. Empirical modelling – using a simple analytical function (e.g. polynomial) to describe a set of
   observations. If the model is appropriate, this can be used to interpolate and extrapolate from a set of
   observations. However, great care should be taken if attempting to predict measures near the sampling
   limits of the dataset (point x1 in figure 2.2) or if predicting any form of integral of the the signal
   requiring extrapolation.
2. Semi-empirical modelling – in this approach, rather than starting from a purely empirical function,
   the point of departure of model development is some physically-based solution. A semi-empirical
   model is developed by simplifying terms or making approximations in the physically-based model,
Vegetation Science – MSc Remote Sensing UCL                                                                      Lewis

    perhaps applying empirical ‘linkages’ between components of the model (e.g. an empirical weighting
    of scattering by different mechanisms).

In this section, we shall develop the following semi-empirical models: (i) a linear approximation to an
optical canopy reflectance model; (ii) a lumped parameter model for estimating canopy cover from
optical measurements; (iii) a lumped parameter microwave model.

3.1 Linear Kernel-driven Modelling of Canopy Reflectance

This form of model was first developed for remote sensing applications by Roujean et al. 6 The motivation
was to provide a model to correct for (‘normalise’) angular variations in observed reflectance from
moderate spatial resolution, high temporal frequency, sensors.

                           Satellite, Day 1                    Satellite, Day 2


                                          Figure 3.1 Wide FOV sampling

The main way in which such sampling is achieved is from (near) sun-synchronous orbits, using sensors
with a wide field of view (FOV). Sensors in this category, with FOVs of around 110o are NOAA
EVISAT has a slightly narrower FOV, but operates in much the same way. Thus, even though the repeat
time for the sub-satellite point may only be of the order of two weeks, daily sampling (or better at higher
latitudes) can be achieved as a very wide swath is sampled by the wide FOV. The problem, as
demonstrated in figure 3.1, is that the viewing vector may be very different for subsequent ‘looks’ at the
surface. The local illumination angle may also change considerably – the local time at the sub-satellite
point is kept approximately constant for a Sun-synchronous orbit, but the local time will vary across a
wide swath.

Since there are generally considerable variations in the surface reflectance due to varying viewing and
illumination geometries (‘bidirectional reflectance effects’ or BRDF effects), this means that the observed
reflectance can vary significantly, even if there are no changes in the surface state. This effect is not
generally compensated for in any spectral transformations (e.g. NDVI).

 J-L Roujean, M. Leroy, and P.Y. Deschamps, 1992, “A bidirectional reflectance model of the Earth’s surface for the
correction of remote sensing data”, Journal of Geophysical Research, vol.97, pp. 20455-20468.
Vegetation Science – MSc Remote Sensing UCL                                                                                                                                     Lewis


                            NDVI   0.25
                                                                                                Julian Day

                                                                  original NDVI                 MVC               BRDF normalised NDVI

                                     Figure 3.2 AVHRR NDVI over Hapex-Sahel, 1992

Figure 3.2 shows a temporal sequence of NDVI data. The data have been screened for cloudy pixels
(using thermal channel information and brightness thresholds), and are plotted as a function of day of
year. The ‘raw’ data (points) show a large degree of day-to-day variation (noise), so that it is not at all
obvious when there is a change in the vegetation state (it is around day 218). A typical method for
compositing AVHRR (or other moderate resolution data) scenes has been the Maximum Value
Composite (MVC). Using this method, the maximum NDVI that occurs within a temporal window is
assumed to be clear of cloud or other contaminations (cloud and cloud shadow reduce NDVI). However,
the resultant signal is ‘stepped’ (as a maximum value is used), and known to be sensitive to various other
contaminations e.g. due to directional reflectance effects.

Roujean et al. developed a semi-empirical model to account for BRDF effects and allow a normalisation
of the signal to a common reference frame (constant viewing and illumination geometry).

The model developed can be stated:

   , ,    f iso    f vol  k vol ,    f geo  k geo ,                                                                                        (3.1)

where   , ,  is the spectral bidirectional reflectance factor, a function of wavelength  , and viewing
and illumination vectors, ,  respectively. f iso   , f vol   , and f geo   are functions of wavelength
and describe the magnitude of the three components of scene scattering. k vol ,   and k geo ,   are
the so-called ‘kernels’ which are functions of viewing and illuminations angles only, and are derived from
abstractions of physically-based models.

It is assumed that the directional reflectance of a scene can be modelled by a weighted combination of
‘volumetric’ scattering (RT-type scattering) and geometric-optics (’geo’) type shadowing/shadow-hiding
effects. Accepting that these are essentially single-scattering mechanisms, the model assumes that
multiple-scattering will be isotropic (i.e. have no directional dependency).

The volumetric kernel is developed from the first-order solution to the radiative transfer equation for
infinitesimal scatterers. Roujean et al. and later Wanner et al. (1995) developed solutions based on
assumptions of low LAI and high LAI.
Vegetation Science – MSc Remote Sensing UCL                                                        Lewis

The solution for first order scattering, assuming equal leaf reflectance and transmittance, for a canopy
with a spherical leaf angle distribution and a Lambertian soil/understorey reflectance of  s is:

                                   
             2 sin   cos  2   
 1 ,   l 
                           
                                                     
                                         1  eX  seX                                   (3.2)

where  l is the leaf single scattering albedo,  is the phase angle between viewing and illumination
vectors,  and   are the cosines of the viewing and illumination zenith angles respectively, and the
‘lumped’ term X is:

      L   
X                                                                                         (3.3)
      2  

where L is the canopy LAI. If L is small, we can write e  X  1  X , so equation 3.2 becomes:

                                   
             2 sin   cos  2    
                                      L            L    
 1 ,   l                          2      s 1  2   
                                                                                       (3.4)
             3                                                

We can further approximate the term involving soil reflectance transmitted through the canopy as zero,
since L is low:
                                    
              2 sin   cos  2    
                                       L     
  1 ,   l                          2      s
                                                                                     (3.5)
              3                               

We can now write this approximation to first order ‘volumetric scattering’ for an optically thin canopy as:

 thin  , ,   a0    a1  k vol ,                                         (3.6)


                                   
                sin   cos  2   
                                      
k vol ,                                                                            (3.7)
                                        2

This ‘kernel’ is defined so that it equals zero at viewing and illumination angles of zero
(   0;      1 ). The ‘lumped’ model parameters are:

            L l
a0            s                                                                      (3.8)

            L l
a1                                                                                     (3.9)
Vegetation Science – MSc Remote Sensing UCL                                                       Lewis

Since the kernel is zero for nadir viewing and illumination zeniths, a0 can be thought of as the nadir
reflectance – a normalised reflectance term. The lumped parameters tell us about the behaviour of single-
scattered reflectance for an optically thin canopy: the reflectance is essentially controlled by the product
of LAI and leaf single scattering albedo, and varies with angle according to the kernel defined in equation
3.7. This kernel is sometime referred to as RossThin, after Ross who first formulated the expression for
first order canopy scattering.

The RossThick kernel is derived in a similar way – for an optically thick canopy, we assume the LAI term
in the exponential to dominate over angular variations, and approximate the exponential:

    L    
    2     exp LB
exp         
             

where B is an average of the angular dependence. In this case, we obtain a similar linearised formula. In
fact, we find that there is not a great deal of difference between RossThick and RossThin, in the sense that
they are approximately linearly related over a range of angles. For operational algorithms, we tend to
chose RossThick as the kernel which fits observed data best most of the time.

We consider geometric-optics scattering (shadowing and protrusion effects) in a similar way – deriving
an approximate formula from the GO models of Li and Strahler (figure 3.3). These models abstract a
canopy as a distribution of spheroids on trunks (‘spheroid on stick’ models). They are parameterised by
the density of the objects, the brightness of ground and crown protrusions, and the geometry of the
spheroids (base-height ratio and tree height-crown height ratios). To obtain a linearised form, we must fix
the spheroid geometries at ‘typical’ values, and assume the crown and ground brightness to be constant
(C). Again, we have a choice of approximations, and typically assume a ‘sparse’ (no overlapping shadows
or projections) and ‘dense’ (overlapping terms) case. The resulting kernels are known as LiSparse and
LiDense. We tend to chose the LiSparse kernel as a generally-applicable model.


                                                                               shadowed crown
           r                                    Sunlit crown
       b                            Projection

                                                                          shadowed ground

                          Figure 3.3 Geometric Optics Model considerations
Vegetation Science – MSc Remote Sensing UCL                                                                                 Lewis

                                                          RossThick        LiSparse



                               -75   -60   -45   -30      -15          0     15       30     45       60       75
                      kernel value




                    Retro reflection (‘hot spot’)
                                                          view angle / degrees

                                     Figure 3.4 Volumetric (RossThick) and Geometric
                                     (LiSparse) kernels for viewing angle of 45 degrees

Figure 3.4 shows a plot of the RossThick and LiSparse kernels in the solar principal plane. Note that
RossThick is generally an ‘upward bowl’ shape, with LiSparse generally downward. As RossThick is
formulated for infinitesimal scatterers, it has no ‘hot spot’. The Li kernels, however, do implicitly contain
a hotspot as they are derived from shadowing considerations. The angular width of the Li kernel hotspot
is fixed, as the protrusion geometry is fixed in the linearisation.

We suppose the scene reflectance to be composed of some proportion (  ) of scattering from geometric-
optics effects (thus 1-  from volumetric effects), and build a scene reflectance model as:

  , ,   1   a0vol    a0 geo     mult    1   a1vol  k vol ,   a1geo  k geo , 

which we see to be the form required in equation 3.1.  is the isotropic multiple scattering term. The
compound isotropic term is the total canopy reflectance at nadir viewing and illumination.

To summarise, a linerarised BRDF model is developed by considering a linear mixture of two main
scattering mechanisms – volumetric and geometric terms, along with an assumed isotropic multiple
scattering term. Physically-based models for the two mechanisms are abstracted and linearised to give a
linear formula. The model parameters are compound (‘lumped’) terms – e.g. functions of LAI and leaf
single scattering albedo. They cannot be used directly to determine canopy biophysical parameters (LAI),
but they do suggest general couplings in canopy scattering.

Similar models have been developed to describe variations in emitted (thermal) radiation as a function of
viewing and illumination geometries.

3.2 Using Linear BRDF Models for angular normalisation
In this section, we will consider some of the uses of the linear models developed above. The original uses
such models were put to were angular normalisation of reflectance from wide FOV sensors, and the
calculation of angular integrals of reflectance – terms related to albedo. These models are used for this
purpose in the operational MODIS BRDF/albedo product.
Vegetation Science – MSc Remote Sensing UCL                                                        Lewis

                 Figure 3.5 MODIS TERRA sampling for a latitude of 20oS in March

As noted above, wide FOV sensors produce high temporal frequency data with inherent variations due to
BRDF effects. For some applications, such as the calculation of vegetation indices, we would like to
normalise for these effects.

Figure 3.6 demonstrates these effects – the upper figure shows MODIS 500m reflectance data for two
pixels over vegetation in Southern Africa (one burned – filled circle, the other unburned – unfilled circle),
plotted as a function of day of year. For any particular sample, we see a large variation (of the order of
0.06) in reflectance between the burned and unburned vegetation. However, the day to day variations are
even larger than this – of the order of 0.1. The lower figure shows the same reflectance data plotted as a
function of view zenith angle – wide field of view sensors tend to sample close to a single ‘plane’ in the
viewing hemisphere (figure 3.5). There is a small amount of solar zenith angle variation, but the plot
becomes smooth plotted in this way – we see that the dominant day to day variations are due to BRDF
effects and the difference between burned and unburned vegetation is clear.
Vegetation Science – MSc Remote Sensing UCL                                                                     Lewis

                             Figure 3.5 MODIS 500m channel 5 reflectance data

The ‘downward’ bowl shape of the observed BRDF is indicative of geometric optics scattering effects –
exactly what we would expect for the sparse vegetation in the area considered.

We obtain a set of model parameters by ‘fitting’ the linear model to the observations. For sensors such as
MODIS, we must gather samples over some time period as we only typically have one sample per day.
This makes the assumption that the surface state does not change over the time ‘window’ of the inversion
– practically, this smooths the signal somewhat – the larger the moving window, the more the smoothing.
We saw in section 2.3 that the error in a prediction increases as we move away from where the
observations were taken. We could use the linear model parameter f iso   as an angular normalisation
term, but the error in this may be large. Instead, we tend to use a normalised measure of reflectance at the
mean local solar zenith angle. This could, for instance, be the reflectance at nadir illuminated at the mean
solar zenith angle7. If we write equation 3.1 as:

      , ,   P  K                                                                              (3.10)


 X. Zhang, M.A. Friedl, C.B. Schaaf, A.H. Strahler, J.C.F. Hodges, F. Gao, B.CV. Reed and A. Huete, (2003), “Monitoring
vegetation phenology using MODIS”, Remote Sensing of Environment, 84:471-475.
Vegetation Science – MSc Remote Sensing UCL                                                          Lewis

      f iso                                 1      
                                                      
 P   f vol                     K   k vol ,  
      f                              k , 
      geo                               geo           

then we can calculate this normalised reflectance by setting the terms in K to those of the required
geometry. The uncertainty in the normalised reflectance is given through equation 2.24, with Q set to K.

3.3 Using Linear BRDF Models for albedo

Directional-hemispherical reflectance,   ,   can be phrased as an integral of BRF for a given
illumination angle over all illumination angles. It is a measure of total reflectance due to a directional
illumination source (e.g. the Sun) and is sometimes called ‘black sky albedo’. Radiation absorbed by the
surface is simply 1-   ,   .

                   1
 ,                , , d                                                    (3.11)
                   0

Given a linear model of the form of equation 3.1, we can write:

 ,   
                   1

                   0
                          f    f  k ,   f  k , d
                           iso         vol     vol           geo   geo

Since only the kernels are functions of viewing angles, we obtain:

  ,   P  K                                                                          (3.12)


                          
               1          
      1                 
K    k vol , d                                                                 (3.13)
        0
      1                 
      0
      k geo ,  d 

Similarly, the bi-hemispherical reflectance    , is a measure of total reflectance over all angles due to
an isotropic (diffuse) illumination source (e.g. the sky). This is sometimes known as the ‘white sky
albedo’. It can be written:

              1  1
                    , ,  dd                                             (3.14)
           0   0

similarly to above, we can write:

    P  K                                                                              (3.15)

where K is the bi-hemispherical integral of the kernels. The spectral albedo,    – the total energy
reflected by a surface in a particular waveband, an important surface parameter in energy budget (e.g
climate) modelling, can be given as:
Vegetation Science – MSc Remote Sensing UCL                                                      Lewis

                                        
             D  P   K  (1  D ) P   K  
                                          
        P   D K  (1  D ) K    P   K 

where K  is a proportionate mixture of the direct and diffuse kernels, weighted by the proportion of
diffuse illumination D   . Angular integrals of the kernels can be pre-computed, so that the spectral
albedo can be rapidly calculated from the linear model parameters P.

In fact, we often find the directional hemispherical integral of reflectance a more useful measure for the
calculation of vegetation indices, so we can use equation 3.12 to calculate a normalised reflectance term
by inserting the angular integral of the kernels for the mean solar zenith angle into K .

Angular normalisation and the calculations of albedo are seen to be examples where we are not
necessarily interested in the estimation of surface biophysical parameters – the emphasis is on providing a
reliable description of the reflectance field from a limited sample set of BRDF measurements.

3.4 Using Linear BRDF Models to track change
Another example where we do not necessarily require a ‘full’ parameterisation of the surface is in
tracking change. The linear BRDF models derived above are being used to develop an operational
algorithm for detecting burn scars from MODIS data. The motivation for using a BRDF model is
provided by the observed variations in reflectance due to BRDF effects which may mask the change we
are trying to detect (figure 3.5). We first examine which MODIS channels are sensitive to the change we
want to detect. MODIS can detect active fires (using thermal information), but this does not pick up all
burning events as a fire may not be burning when observed. Thus we look for evidence of change in the
surface reflectance.

     Figure 3.6 Mean and +/- 1 SD for burned (filled circles) and unburned (open circles) pixels
Vegetation Science – MSc Remote Sensing UCL                                         Lewis

     a)                                                b)
Figure 3.7. MODIS Channel 5 reflectance for days 275 (a) and 277 (b) over Namibia

Figure 3.7. MODIS Channel 5 reflectance for days 275 (a) and 277 (b) over Namibia

  Figure 3.8 MODIS Channel 5 predicted reflectance (277) and Z-score pseudocolour
Vegetation Science – MSc Remote Sensing UCL                                                       Lewis

Figure 3.6 shows summary statistics for a sample set of more than 15000 pixel samples before (unburned)
and after (burned) a fire event was picked up by the active fire detection algorithm. Note that the samples
contain some BRDF effects. We can see form this that channels 2, 5, and 6 (all NIR) provide the greatest
separability between burned and unburned pixels – the signal is not strong in the other channels. We
choose to develop an algorithm based mainly on channel 5 reflectance, as this has the strongest

Figure 3.7a shows an observation in channel 5 over Namibia for day 275 of 2001 – recent burning is
evidenced by the dark non-linear features. Figure 3.7b shows the same area two days later. Some areas of
burning are evident between the two dates, but BRDF variations mean that they cannot all be reliably

To use a BRDF model to detect change, we invert the linear model against observations leading up to an
observation, and make a prediction of the reflectance based on that evidence, at the viewing and
illumination angles of the new observation. We also predict the uncertainty in the prediction from
equation 2.26. We then define a measure of disassociation of the observed and predicted measures – a 
measure (this was originally done with a Z-score):

                     predicted                        predicted 
                                     
          observed                              observed

                e2   2                                       1
                                                    e 1

which will be a high negative value if there is a negative change in reflectance. The uncertainty measures
account for noise in the data in this comparison and uncertainty in the prediction. It is effectively a
measure of the overlap of the observed and predicted Normal distributions.

Figure 3.9 shows the tracking of a single pixel over time – the algorithm follows the observed reflectance
variations well until day 275 (when the burn occurred). The discrepancy measure is seen to have a large
negative value at this point, which then stays negative for around a week. We detect this change as being
of long duration, and infer that it is a burn scar. Short-term unexpected changes are attributed to residual
clouds and cloud shadow.
Vegetation Science – MSc Remote Sensing UCL                                                                    Lewis

      Figure 3.9 BRDF tracking showing predicted (filled triangle) and observed (open triangle)
         reflectance for a single pixel in MODIS channel 5, near the Angola-Namibia border

The algorithm produces ‘maps’ of the day of burning. Figure 3.10 shows an example for the area
considered. A clear progression of the fire fronts is shown.

                                    Figure 3.10 – rainbow scale of day of burning

The importance of the algorithm over previous efforts at burn scar detection is its use of a BRDF model to
predict reflectance – we are not primarily interested in the surface biophysical parameters, but we do want
the algorithm to work rapidly. A linear lumped parameter BRDF model is very suitable for such a task.

3.5 Other Lumped Parameter Optical Models
The kernel-based linear BRDF models developed and used above are not the only lumped parameter
formulations available. Another example which can be used in much the same way as above is the
Modified RPV (MRPV) model. This consists of three ‘scattering mechanism’ terms, grouped together as
a product. When a logarithm of the model is taken, the model is linearised and can be inverted

Gilabert et al8 develop a simplified canopy reflectance model. Canopy reflectance is written:

          s       exp  CL 

where L is the canopy LAI,     is an optically-thick canopy reflectance,  s   is the soil reflectance,
and C is a lumped parameter depending on clumping and leaf orientation. In one sense this isa simple
linear mixture model:

    (1  f )      f s  

 M.A. Gilabert, F.J. García-Haro, and J. Meliá, (2000), “A Mixture Modeling Approach to Estimate Vegetation Parameters for
Heterogeneous Canopies in Remote Sensing”, Remote Sensing of Environment, 72:328-345
Vegetation Science – MSc Remote Sensing UCL                                                                       Lewis

where f is the fractional soil cover equal to exp(-CL).     is defined in terms of leaf single scattering
albedo  l   :

     Al  B

where A and B are parameters related to canopy structure, with A determining first-order scattering and B
describing multiple scattering effects.

The model is in many ways similar to the volumetric scattering kernel developed above. It has three
‘lumped’ parameters (A,B,C), which will be a function of canopy structure (LAI, leaf angle distribution
etc.) and maintains a connection to leaf and soil optical properties. The lumped parameters can be
calibrated for a particular canopy (from measurements or a more complex model) and the simple formula
inverted from observations. It is not, however, a linear model (beyond being a linear mixture model).
Neither does it contain any explicit description of viewing and illumination angle effects. It does,
however, provide a linkage between a linear mixture model and canopy reflectance models.

3.6 A Lumped Parameter Microwave Model
The most commonly-used lumped parameter model for microwave applications is the ‘water cloud’
model, originally due to Attema and Ulaby9. A generalised form of the model is given by Champion et
al.10. It derives from consideration of the radiative transfer equation. The backscattered power, P, can be
given as:

                 2 h     2 h 
P                    S exp   
          1  exp                 
     2                          

Here, S is the radar scattering from the soil,  represents scattering and  attenuation properties of the
canopy, with canopy height h. Note that this form is very similar to what we obtain in the optical case. If
we lump terms together, we can write:

                        2 h
a                b
     2                   

and insert a dependency on LAI, then:

   10 log a1  e bL   Se bL 

Champion et al. argue for a slightly more generalised form:

   10 log aLe 1  e bL   Se bL 

and note that S can be approximated as:

S  C  Dm v

 E.P.W. Attema and F.T. Ulaby, (1978), “Vegetation modelled as a water cloud”, Radio Science, 13:357-364.
  I. Champion, L. Prevot, and G. Guyot, (2000), “Generalised semi-empirical modelling of wheat radar response”,
International Journal of Remote Sensing, 21(9):1945-1951.
Vegetation Science – MSc Remote Sensing UCL                                                                         Lewis

where m v is the volumetric soil moisture. The resulting model is able to mimic variations in observed
backscatter dependencies on soil moisture and LAI. However, the model parameters (a, b, C, D, e) vary
for different canopies, as canopy backscatter depends on more terms than just LAI and soil backscatter on
more than moisture. Effectively the model uses ‘calibration’ of the lumped parameter terms to hide the
fact that biophysical parameters will be correlated (e.g. LAI and leaf size, number density etc.). It also
need re-calibration for different soil textures and roughness measures (parameters C and D). All terms
need to be re-calibrated for different wavelengths and polarisations.

Model calibration is generally carried out using field observations, although, as with the optical lumped
parameter model given above, this could be achieved using more complex microwave scattering models.

The model is useful for localised applications (e.g. for the same crop and soil conditions as were used in
the calibration) and for examining and explaining11 the relative contributions of soil and canopy scattering
terms, and for describing when saturation (lack of dependency on model parameters) occurs. Note that
even when calibrated, the model has two unknowns: LAI and soil moisture, so that it cannot be inverted
from a single channel observation (e.g. ERS SAR) unless one of the parameters is assumed known. If the
cosine term is not included in the lumped definition of parameter b, it can be used with multi-angular
measurements to invert the parameters.

It is likely that such models will find new applications in parameter estimation from the new generation of
multi-channel SARs such as ASAR on ENVISAT.

4. Conclusions
This section of the course has developed a range of abstracted treatments of radiative transfer (and
associated) theories to develop simple (often linear) models of canopy reflectance and scattering. These
semi-empirical models are clearly not as accurate as ‘full’ solutions, but they contain fewer parameters,
and are more tractable to invert against limited remote sensing observations. It is seen that for certain
tasks, such as the calculation of albedo, angular normalisation, and the detection of change, simple
lumped parameterisations which describe expectations in the scattered signal (e.g. as a function of
viewing and illumination angles) are sufficient models to achieve practical tasks in remote sensing.

  Y. Inoue, T. Kurosa, H. Maeno, S. Uratsuka, T. Kozu, K. Dabrowska-Zielinska, and J. Qi, (2002), “Season-long daily
measurements of multifrequency (Ka, Ju, X, C, and L) and full-polarisation backscatter signatures over paddy rice field and
their relationship with biological variables”, Remote Sensing of Environment, 81:194-204.

Shared By: