at Optical and Microwave wavelengths applied to vegetation canopies: Part 2
 
 

P. Lewis & P. Saich, RSU, Dept. Geography, University College London, 26 Bedford Way, London WC1H 0AP, UK.
 
 

Email: plewis@geog.ucl.ac.uk, psaich@geog.ucl.ac.uk
 
 

3. The Radiative Transfer Equation

3.1 The Radiative Transfer Equation

Radiative transfer models have been used extensively since the 1960s to model scattering from canopies at optical wavelengths (Ross, 1981). This approach first exploited in the microwave scattering context during the 1980s.
 
 

The models take as a starting point consideration of energy balance across an elemental volume. This links energy into the volume (either energy incident in the propagation direction, or energy that is scattered from other directions) and energy losses from the volume (either scattering out of the propagation direction, or absorption losses). Whilst optical modelling generally exploits a scalar radiative transfer equation, in microwave scattering, we deal usually with a vector of intensities (typically using the modified Stokes vector defined in section 1.4). The reason for this is that the propagating waves can have well-defined polarisation, and orthogonal polarisations are coupled by depolarising processes - therefore, we cannot consider radiative transfer equations for polarised waves separately from one another. Note that the intensities we are using are not themselves vectors - the introduction of a vector of intensities is only a mathematical convenience.

3.1.1 The Scalar Radiative Transfer Equation

The (one-dimensional) scalar radiative transfer (SRT) equation for a plane parallel medium (air) embedded with a low density of small scatterers defines the change in specific Intensity (Radiance) I(z,W) at depth z in direction W at any given wavelength with respect to z through:
 
 

(3.1)
 
 
 
 

(3.2)
 
 

3.1.2 The Vector Radiative Transfer Equation

The vector form of the vector radiative transfer (VRT) equation expresses the transfer of radiation through an elemental volume in the same way as the scalar form. Vectors are used to express polarisation coupling, as noted above. Depolarisation of incident horizontal and vertically polarised waves is an important part of the remote sensing signal. For a linearly polarised wave, the cross-polarised response from a large conducting sphere is zero (Ulaby and Elachi, 1990; p. 34). In atmospheric LiDAR sensing, for instance, a measure of the degree of depolarisation allows liquid and solid phases of water in the atmosphere to be distinguished. Spherical particles, such as wet haze, fog, cloud droplets, and small raindrops do not significantly depolarise, whereas non-spherical particles such as ice crystals, snow flakes or dust particles have a much higher degree of depolarisation.
 
 

The vector form of the radiative transfer equation is most conveniently expressed using Stokes vectors. In one dimension (after Ulaby and Elachi, 1990; p. 136):
 
 
 
 

(3.3)
 
 
 
 

(3.4)
 
 

In order to fully specify our scattering problem, we therefore need to define the forms of the extinction and phase matrices that appear in the VRT equation and solve the coupled integro-differential equations. The extinction and phase matrices are formed for an ensemble of canopy elements by a summation of extinction and Mueller matrices for different canopy element types (potentially of different size, shape and orientation). This can be otherwise expressed as the product of the number density and a proportion-weighted average of the terms in the extinction and Mueller matrices.
 
 

This is where the solution to Maxwell's equations reappears - in defining the scattering properties of individual particles for inclusion within the 'harness' of the radiative transfer formulation. As we have seen earlier, it is usual to encounter particular approximations in the course of this. Examples of this include far-field scattering, or small, thin or planar scatterers. The phase matrix is closely related to the scattering cross-sections of the 'average particles'.
 
 

3.1.3 Intrinsic Canopy Properties

In modelling canopy scattering, we typically wish to state the scattered quantity as an intrinsic property of the canopy, rather than stating a scattered intensity as a function of incident intensity. This allows us, for instance, to compare measurements made under differing illumination intensities.
 
 

The fundamental intrinsic scattering quantity at optical wavelengths is known as the Bidirectional Reflectance Distribution Function (BRDF) (sr^-1):
 
 

(3.5)
 
 

where: represents polarisation of the receive/transmit wave (x=r or i, p_x=h,v); F_i is the irradiance (Wm^-2) on the surface; and I_r is the radiance (Wm^-2sr^-1) (Tomiyasu, 1988). The BRDF of an ideal diffuse (Lambertian) surface is  (for an unpolarised reflector) and is independent of viewing and illumination angles (=for a polarised detector). As this is defined for an infinitessimal soild angle, it is more usual to use the Bidirectional Reflectance Factor (BRF)  (with implicit wavelength etc. depoendence). This can be defined as the ratio of radiance leaving the surface around direction _r, , due to irradiance to the radiance on a flat totally reflective Lambertian surface under the same illumination conditions, i.e.:
 
 

(3.6)
 
 

The fundamental intrinsic property for the microwave case is the differential scattering coefficient  (Tomiyasu, 1988), a dimensionless term. This is generally used as , i.e. the radar backscatter coefficient, often denoted  for p=p_r, q=p_i. This is typically specfied in dB (log power ratio) and is defined as for v or h transmit/receive as (Ulaby and Elachi, 1990; p. 138):
 
 

(3.7)
 
 

for an area A at distance r (in the far field). It is more conveniently calculated from the first two components of the Stokes vector as:
 
 

(3.8)
 
 

where is the angle between the surface normal and the transmit/scattering direction.
 
 

Tomiyasu (1988) notes a formal relationship between the differential scattering coefficient and BRDF for a diffuse surface as:
 
 

3.2 Extinction Coefficient and Beer?s Law

The volume extinction coefficient,  (the ?total interaction cross section?, ?extinction loss? or ?number of interactions? per unit length) is a measure of attenuation of radiation in a canopy (or other medium). For a scalar radiation Intensity I (Radiance or Brightness, in Wm^-2sr^-1) travelling in a homogeneous medium of randomly located scatterers, the Intensity is exponentially attenuated over a distance l:
 
 

(3.9)
 
 

where I₀ is the Intensity at l=0. Equation 3.9 is known as Beer?s Law (also the Beer-Lambert Law). From it, we can see that:
 
 

(3.10)
 
 

We can see equation 3.2 as a no-source statement of the SRT equation - zero order scattering solution. can be defined as a function of travel of the direction of the radiation, W (Fung, 1994):
 
 

(3.11)
 
 

can also be expressed in terms of volume absorption and scattering coefficients, and respectively through:
 
 

(3.12)
 
 

where dependence on W is implicit and on polarisation allowed, for all terms. These two terms represent loss due to absorption by the particles (leaves) and scattering by the particles away from the direction of propagation (Fung, 1994; p11). They are related to number density through particle absorption and scattering cross sections similarly to above.
 
 

We also note the related term ?optical thickness?,  (Fung, 1994; p.16):
 
 

(3.13)
 
 

which is often used in radiative transfer formulations.
 
 

A further term, the single scattering albedo of a particle, , can be defined as the probability of radiation being scattered rather than ?extinguished?:
 
 

(3.14)

3.2.1 Optical Extinction Coefficient for Oriented Leaves

For optical wavelengths comprising oriented objects (leaves) which are large compared to the wavelength of the radiation, we can define an effective ?particle? extinction cross section Q_e(W ) in terms of leaf area as:
 
 

(3.15)
 
 

ignoring any dependence of canopy properties on z, for constant leaf area A_l.  is the foliage area orientation function (the ?G-function?), a dimensionless geometry factor equal to the projection of a unit area of foliage on a plane perpendicular to the direction , averaged over elements of all orientations:
 
 

(3.16)
 
 

For a spherical leaf angle distribution (), this is simply 0.5. It can be similarly shown that for an azimuthally independent distribution, for the special case of a completely horizontal distribution (,  in an elliptical description), . For a completely vertical distribution,  (Ross, 1981).
 
 

Myneni et al. (1989; p. 29) show G-functions for a variety of measured canopies. After Ross (1981), they comment that:
 
 

• the range of G-functions for the measured canopies is relatively small (0.3-0.8) and is considerably smoother than the measured leaf inclination distributions;
• for near planophile canopies, the G-function is high (>0.5) for low and low (<0.5) for high 

• the G-function is always close to 0.5 between 50 and 60 degrees, being essentially invariant at 0.5 over different leaf angle distributions at 

Equation 3.16 expresses the G-function as a simple geometric attenuation factor for ?blocking? of unpolarised radiation. Note that it does not vary with wavelength, but that it will generally vary with . From equations 2.1, 3.11 and 3.15, we can write:
 
 

(3.17)
 
 

which is the more usual way of expressing the extinction coefficient for canopies at optical wavelengths (Ross, 1981). Inserting this into equation 3.9, we can see that the attenuation of radiation is exponential and controlled by path length l, leaf area density, and the normalised leaf projection function. In a plane parallel canopy, l can be expressed as:
 
 

(3.18)
 
 

where  (figure 3.1). Note that z is defined from 0 at the top of the canopy to -H at the soil layer in these notes (figure 3.1). The optical depth at the bottom of the canopy (z=-H) is then:
 
 

(3.19)
 
 

i.e. the radiation at the bottom of the canopy is for a spherical leaf angle distribution, for a horizontal distribution.
 
 

Note that for the optical case, we typically define the single scattering albedo of a particle (leaf), w(l), which varies as a function of wavelength, through reflectance, r_l(l) and transmittance, t_l(l), so
 
 

The optical case is ?straightforward? and relatively intuitive, as we can use geometric optics principles to consider blocking (interception) of the radiation in defining the extinction coefficient or optical depth. For a more general case, we turn to the example of a so-called Rayleigh medium within which the particles are assumed spherical, of radius a. It is further assumed that the particle radius is small relative to the wavelength in the dielectric material of the particle. This condition can be stated as:
 
 

where k is the wavenumber (), with in the same linear units as a (Ulaby and Elachi, 1990; p.141). This condition holds for scattering by gasses in the atmosphere (Fröhlich and Shaw, 1980) and for long wavelength microwave (P-band) extinction by vegetation canopies. Where this is not the case, we must apply modifications to the extinction and scattering coefficients, as shown for discs and needles in section 1.3 (equations 1.7 and 1.8). We do not attempt to derive the extinction coefficient here, rather we state it and examine its properties. A Rayleigh scattering layer can be characterised by two parameters: the single scattering albedo  (equation 3.14); and the optical depth (equation 3.13) for a path length d for constant .
 
 

For a canopy in air, the scattering coefficient for a random medium of spherical Rayleigh scatterers is (Fung, 1994; p. 122-124; Ulaby and Elachi, 1990; p. 142):
 
 

(3.21)
 
 

where f is the volume fraction of scatterers (), (=) is the dielectric constant of the sphere material (leaves) and  is the wavenumber in air. The absorption coefficient is:
 
 

(3.22)
 
 

From this we see that:
 
 

(3.23)
 
 

We can see that the single scattering albedo increases for increasing leaf radius (figure 3.3) and decreases with increasing wavelength. Equation 3.23 shows that it is effectively a function of the ratio of leaf linear dimension to wavelength, modulated by leaf moisture (dielectric constant). As the leaf dielectric constant is a function wavelength, the relationship of single scattering albedo with normalised leaf radius (leaf radius over wavelength) is not obvious, but figure 3.4 (with dielectric constants for 70% GMC) demonstrates that the impact of the variation in dielectric constant with wavelength is typically small compared to the overall behaviour.
 
 

Note that when the leaf radius is more than about 0.04 times the wavelength, the single scattering albedo is no longer small (?small? being conidered as 0.2 or less here). The condition for validity of the model given above implies a/l << 0.16, so a threshold of 0.04 would seem generally reasonable. We can therefore state that for leaf dimensions and wavelengths within the region of validity of the Rayleigh scatterer assumption, the single scattering albedo is low.
 
 

The optical depth is given by:
 
 

(3.24)
 
 

at the bottom of the canopy (z=H), where m is the cosine of the zenith angle (figure 3.1). We see that the optical thickness of the canopy depends on wavenumber, the size of the scatterers, the volume of scatterers per unit area (fH) and the real and imaginary parts of the leaf dielectric constant.

Figure 3.5 shows the dependence of normalised optical depth (tm/(fH)) on leaf radius. The general trend is an increase in optical depth with increasing leaf radius.

Figure 3.6 shows the same information, on a linear scale, for leaf sizes limited to within the assumed range of validity of the model (w<=0.2). Note that the rate of increase in optical depth increases significantly with decreasing wavelength, and that it is essentially flat for L- and P-band over the range of validity of the model.
 
 

For spherical scatterers, the extinction matrix  is equal to a 4x4 identity matrix multiplied by . Since the extinction term is effectively scalar, we can use Beer?s Law to express extinction within the Rayleigh medium (see equation 3.9):
 
 

(3.25)
 
 
 
 

 
 

Figure 3.6 Normalised optical depth for 70% leaf GMC as a function of leaf radius for range of validity of model (threshold at w=0.2)

The phase matrix , of a Rayleigh medium is defined as a function of incident and scattering polar angles () with subscripts i and s respectively is (Ulaby and Elachi, 1990):
 
 

(3.26)
 
 
 
 

where

(3.27)
 
 

The defined phase matrix is dependent on the incident and scattering geometry (but not on the orientation of the scatterers as they are spherical) and on . This depends in turn on  (thus essentially on leaf water content), on k⁴ (so l^-4 - scattering decreases as wavelength increases), and on the total volume of scatterers (product of fractional volume and a³).
 
 

3.2.3 A Solution to the Vector Radiative Transfer Equation

For a medium of low scattering (single scattering albedo, w, is low) we can use the iterative method to solve the RT equation (vector or scalar). At visible wavelengths, w is generally low, but it increases dramatically to close to 1.0 in the near infrared (see figure 2.6). In the microwave, w (see section 3.2.2 for Rayleigh scatterers) can often be small, particularly for low leaf moisture, small leaves, or long wavelengths.
 
 

Under these conditions, a solution to the VRT (or SRT) equation can be obtained by first computing a zero-order solution, then plugging the resultant terms in as source terms to compute first order scattering (Ulaby and Elachi, 1990; pp. 146-151). Second+ order scattering terms are then calculated by inserting the results of the previous order as source terms. Under the condition that w is low, the terms rapidly converge to zero for increasing scattering order. Generally, only up to second order terms are used with this method.
 
 

Considering a two-stream model concerning the propagation of upward- and downward-going intensities, I⁺ and I^- respectively, we can write the VRT equation as:
 
 

(3.28)
 
 

with source functions:
 
 

(3.29)
 
 

subject to boundary conditions (figure 3.9):
 
 

(3.30)
 
 

where is a specular reflectivity matrix for a ?flat? soil, describing scattering at the interface between medium 2 (air) and medium 3 (soil), and d(x) is a delta function which is unity for x=0 and zero elsewhere. It is obtained by exchanging subscripts 23 for 12 in equation 2.6.
 
 

Following Ulaby and Elachi (1990; p.148) we rephrase the differential equations in 3.28 as integral equations:
 
 

(3.31)
 
 

where T(-x)=e^-x is an attenuation (extinction) term due to Beer?s Law.

In the upper equation we see (figure 3.7) the upward intensity defined by an attenuation of the upward intensity at the canopy bottom over distance z-(-H), along with an integral (sum) from the canopy base to height z of the upward source term, attenuated over vertical distance z - z?.
 
 

In the lower equation we see the downward intensity at the canopy top is attenuated over the distance z, with the addition of the integral of a source term from canopy top to height z attenuated over a distance z-z? (figure 3.8).
 
 

Inserting the boundary conditions (figure 3.9):
 
 
 
 

(3.32)
 
 

we see the upward intensity now to consist of:
 
 

1. a double attenuation of the incident intensity (over distance z-(-H) at angle m _s and distance -H at angle m ₀);
2. an integral of the downward intensity (attenuated over distance -z-z? within the integral) over the canopy height, reflected at the canopy base an attenuated over distance z-(-H);
3. an integral from the canopy base to height z of the upward intensity attenuated over distance z-z?.

3.2.3.1. Zero-Order Solution

The zero order solution is arrived at by setting the source terms in equation 3.28 to zero (Ulaby and Elachi, 1990; p.148):
 
 

(3.33)
 
 

The downward intensity is the incident intensity attenuated by a Beer?s Law term over height z. The upward intensity comprises the incident intensity modulated by a downward attenuation over height H, reflectance at the ground interface, and attenuated over height z-(-H).
 
 

3.2.3.2. First-Order Solution

Using the iterative method, we obtain the solution to first-order scattering by inserting the zero-order solution as the source terms in equation 3.29.
 
 
 
 
 
 

(3.34)
 
 

The upward intensity for first order scattering (Ulaby and Elachi, 1990: p. 149) is:
 
 

where k₁= k_e(1/m₀+1/m_s); k₂= k_e(1/m₀-1/m_s). The solution for upward intensity at z=0, for m_s=m₀, f_s=p+f₀ is:
 
 
 
 

(3.35)
 
 

Noting f or small , so for . We can recognise the various terms in this result as:
 
 

1. 

 

 
 

direct intensity reflected at the canopy base, doubly attenuated through the canopy;
 
 



Figure 3.10 First-order microwave scattering mechanisms




2. 

 

 
 

A ?double bounce? term involving a ground interaction, (downward) volumetric scattering by the canopy, and an additional ground interaction. The term includes a double attenuation (on the way to the ground the first time & back from the ground at the end) and attenuation as part of the volume scattering. This term is generally very small and is often ignored.
 
 

3. 

 

 
 

A downward volumetric scattering term, followed by a soil interaction, including a double attenuation on the upward and downward paths.
 
 

4. 

 

 
 

A ground interaction followed by volumetric scattering by the canopy in the upward direction, again including double attenuation.
 
 
 
 

5. 

The phase function for spherical scatterers is given in equations 3.26 and 3.27. We noted that the extinction coefficient for spherical scatterers was essentially scalar, allowing us to treat it simply in the formulation above. In the general case both the extinction terms and the phase matrix will also be functions of scatterer orientation, and may include cross-polarised components.
 
 

The backscatter coefficient is given by Ulaby and Elachi (1990; p. 150):
 
 

(3.36)
 
 

The second order solution is similarly obtained by setting the first-order solutions as the source terms. Note that whilst there are no cross-polarisation terms for spherical (Rayleigh) scatterers, they do occur in second+ order scattering, i.e. cross-polarisation for spherical scatterers is a result of multiple scattering. The second order solution is given in Ulaby and Elachi (1990; pp. 150-151).
 
 
 
 

3.2.4 A Solution to the Scalar Radiative Transfer Equation
 
 

It is interesting to note at this point the term ?first-order interactions? used in the microwave literature to describe these five terms. Essentially, it is used to express terms which have a single interaction with the vegetation canopy. In the optical canopy modelling literature, only terms 1 and 5 would be included in a first-order scattering model (the other terms having multiple interactions with soil and vegetation).
 
 

Other than this, we can proceed for the scalar approach for optical wavelengths in much the same way as section 3.2.3. One other issue we must be aware of is that we are not (generally) interested in a backscattering term - we need to keep a dependence on incident and scattered intensities in the fomulation.
 
 

We have already dealt with the extinction coefficient for the optical case (section 3.2.1), where we noted  (equation 3.17). We will therefore generally also need to allow this extinction term to vary with the angle of the incident radiation.
 
 

The phase function at optical wavelengths is often expressed as:
 
 

(3.37)
 
 

where u_l is the leaf area density, m? the cosine of the incident illumination zenith angle, and G the area scattering phase function.
 
 

(3.38)
 
 

where r_l and t_l are the leaf directional reflectance and transmittance factors respectively (Ross, 1981). We can see this as a double projection of the leaf angle distribution, modulated by reflectance for the upper hemisphere and transmittance for the lower hemisphere, in much the same way as the G-function expresses a projection in the direction of the radiation for the extinction coefficient.
 
 

A typical assumption used in building a canopy reflectance model is to let the leaves be biLambertian scatterers. In this case, the angular dependence is removed from the reflectance and transmittance, which can be taken outside the integral. This is sometimes modified by the addition of a ?rough surface? specular component from the leaves (Nilson and Kuusk, 1989). If the reflectance is assumed equal to transmittance (or a linear function thereof) we can more simply express the spectral dependence in terms of the single scattering albedo and a weighting of the upper and lower hemisphere integrals.
 
 

To formulate for first order scattering in the optical (scalar) case, we consider first-order interactions (as in only one interaction with canopy or soil elements) in a scalar expression of equation 3.32:
 
 

(3.39)
 
 

(3.40)

(3.41)
 
 

Performing the integral to give the intensity at z=0:
 
 

 
 

(3.42)
 
 

(3.43)
 
 

(3.44)
 
 

1. 

 

 
 

radiation travelling through the canopy to the soil, being reflected, and travelling out of the canopy in the scattering direction. The radiation is subject to a double attenuation due to the two paths through the canopy. The (Beer?s Law) attenuation depends on the (zenith) angles of viewing and illumination and on the leaf projection functions. The exponential term includes a dependency on LAI - one-sided leaf area per unit ground area, so we can see the soil contribution as small for high LAI (or high zenith angles of high G-functions), but small for low LAI.
 
 

2. 

We note that for the special case of a spherical leaf angle distribution for biLamberian scatterers:
 
 

(3.45)
 
 

Further, if we assume reflectance to be linearly related to transmittance ():
 
 

(3.45)
 
 

(3.46)
 
 

Note that the single scattering albedo for leaves is not always low. This is particularly true in the near infrared. In this case, we still need to account for further orders of scattering. The iterative method used above for the microwave case is not directly appropriate, but a range of methods exist to solve for multiple scattering. In general terms, we can state that the multiple-scattered component will be high for high LAI and high single scattering albedo. As a multiple-scattered term, it is less dependent on the specific leaf, viewing, and illumination geometries, and is typically an upward shallow ?bowl? shape when plotted as a function of viewing zenith angle.
 
 

3.2.4 Modifications to a Simple radiative Transfer Approach
 
 

3.2.4.1 Optical Canopy Hotspot

Under the far field approximation, where scatterers are distant from one another, the optical gap probabilities of incoming and outgoing radiation can be assumed independent and can be described by Beer?s Law. Thus, the joint gap probability, Q is:
 
 

 
 

However, when the viewing vector is leaving in the same direction as the incident solar radiation, using this equations we obtain:
 
 

 
 

 
 

 
 

in the retro reflection direction. Away from the retro reflection direction, we require:
 
 

From a mathematical standpoint then, we have the basis of a functional form for this correction factor. The enhanced joint gap probability in the retro reflection direction gives rise to a localized peak in reflectance. This feature is known as the hot spot. The angular width of the joint gap probability (of the hot spot feature) can be shown to be a function of the ratio of the average leaf size to canopy height - a dimensionless 'roughness' parameter.
 
 

In studying soils or other rough surfaces, the angular width is again found to be directly related to a dimensionless measure of roughness.
 
 

There have been many formulations of a hot spot correction factor, both from an empirical standpoint and from physically-based modelling. The model of Pinty et al. (1989) is physically-based, being derived from considerations of the common volume between two cylindrical voids in the canopy. The depth at which a cylinder of radius r from the viewing direction will overlap with one from the illumination direction is denoted L_i. A more physically-consistent hot spot model, provided by Nilson and Kuusk (1989), is based on considerations of random overlapping of horizontal circular disc leaves.
 
 

Other hot spot formulations exist: the vast majority are simply variants on the themes presented above, but many do not consider reciprocity of the joint gap probability well, preferring to opt for empirical 'correction' terms.
 
 
 
 

References

Jacquemoud, S., Ustin, S.L., Verdebout, J., Schmuck, G., Andreoli, G., and Hosgood, B., 1996, "Estimating leaf biochemsitry using the PROSPECT leaf optical properties model", Remote Sensing of Environment, 56, 194-202.

Jacquemoud, S., and Baret, F., 1990, "PROSPECT: a model of leaf optical properties spectra", Remote Sensing of Environment, 34, 75-91.

El-Rayes, M.A., and Ulaby, F.T., 1987, "Microwave dielectric spectrum of vegetation - Part 1: Experimental observations", IEEE Trans. Geoscience and Remote Sensing, GE-25, 541-549.

Chuah, H.T., Lee, K.Y., and Lau, T.W., 1995, "Dielectric constants of rubber and oil palm leaf samples at X-band", IEEE Trans. Geoscience and Remote Sensing, GE-33, 221-223.