Agronomy Journal Grow Your Career With ASA
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published in Agron J 99:1278-1287 (2007)
DOI: 10.2134/agronj2006.0211
© 2007 American Society of Agronomy
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Shirazi, M. A.
Right arrow Articles by Reporter, M.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Shirazi, M. A.
Right arrow Articles by Reporter, M.
Agricola
Right arrow Articles by Shirazi, M. A.
Right arrow Articles by Reporter, M.
Related Collections
Right arrow Inverse Procedures/Parameter Estimation
Right arrow Remote Sensing
Right arrow Other Integrated Systems

Remote Sensing

A Diurnal Reflectance Model Using Grass

Surface-Substrate Interaction and Inverse Solution

Mostafa A. Shirazia,* and Minocher Reporterb

a Western Ecology Division, NHEERL, U.S. Environmental Protection Agency, 200 SW 35th St., Corvallis, OR, 97333
b Botany and Plant Pathology, Oregon State Univ., Corvallis, OR 97331-7306

* Corresponding author (Shirazi.Mostafa{at}epa.gov)


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 
The accuracy of using remote sensing data from earth orbiting radiometers can be improved by using a model that helps to separate the green-fraction in a canopy reflectance ({rho}) from thatch and soil background, accounts for their diurnal changes, and inverts to a solution of a biophysical plant property of interest. Previous studies addressed one or more of these needs separately. Because reflectance components are interdependent, difficulties remain in obtaining a combined inverse solution. We combined a conditional probability method with a novel experimental procedure to predict grass dry weight (dw). Using simple ratio (SR) = near-infrared reflectance (NIR)/red that varied with the normalized time T = (local time – sunset)/daylength, we predicted the mean differential grass dry weight, {Delta}dw = dw1 – dw2, of two grass patches. SR1 and dw1 defined the first patch, which included a background and the second, SR2 and dw2, a predefined background. The inverted solution for {Delta}dw was an ellipse with axes formed by the diurnal reflectance SR1 and SR2 coordinates. It described previously studied soil line and the zone of canopy x canopy x ground interactions. The standard error of predicting {Delta}dw was 17%. We separately tested for plant height SR = f(h) or fresh weight SR = f(fw) using SR, and for dw as a function of normalized difference vegetation index NDVI = f(dw). SR = f(dw) produced superior results. Potential applications include noninvasive prediction of other biophysical plant properties in a single or in hyperspectral bidirectional reflectance in agronomic and ecological remote sensing.

Abbreviations: B, biophysical parameter for Canopy 1 (B1) or Canopy 2 (B2) • c1, clipped once • c2, clipped twice • cg, canopy cut to ground level • cstdv, conditional standard deviation • dw, biomass dry weight • E(|) conditional expectation • fw, biomass fresh weight • h, canopy height • ir, index of relationship • LAI, leaf area index • MSS4–MSS7, four Landsat wavelength bands • n1, Natural Canopy Number 1 • NDVI, normalized difference vegetation index • NIR, near-infrared reflectance • RI, reflectance index • SR, simple ratio • stdv, unconditional (commonly used) standard deviation • T, normalized local time • {Delta}dw, differential biomass dry weight • {theta}s, solar zenith (or elevation) angle • {theta}v, viewing zenith angle • {rho}, canopy reflectance ({Gamma}grass/{Gamma}Labsphere) • {varphi}s, solar azimuth angle • {varphi}v, viewing azimuth angle • {Gamma}, canopy radiance (W m–2 sr–1) • {sum}pxq, variance-covariance matrix p rows by q columns

Received for publication July 18, 2006.

A Diurnal Reflectance Model Using Grass

Surface-Substrate Interaction and Inverse Solution

Mostafa A. Shirazia,* and Minocher Reporterb

a Western Ecology Division, NHEERL, U.S. Environmental Protection Agency, 200 SW 35th St., Corvallis, OR, 97333
b Botany and Plant Pathology, Oregon State Univ., Corvallis, OR 97331-7306

* Corresponding author (Shirazi.Mostafa{at}epa.gov)

Received for publication July 18, 2006.
The accuracy of using remote sensing data from earth orbiting radiometers can be improved by using a model that helps to separate the green-fraction in a canopy reflectance ({rho}) from thatch and soil background, accounts for their diurnal changes, and inverts to a solution of a biophysical plant property of interest. Previous studies addressed one or more of these needs separately. Because reflectance components are interdependent, difficulties remain in obtaining a combined inverse solution. We combined a conditional probability method with a novel experimental procedure to predict grass dry weight (dw). Using simple ratio (SR) = near-infrared reflectance (NIR)/red that varied with the normalized time T = (local time – sunset)/daylength, we predicted the mean differential grass dry weight, {Delta}dw = dw1 – dw2, of two grass patches. SR1 and dw1 defined the first patch, which included a background and the second, SR2 and dw2, a predefined background. The inverted solution for {Delta}dw was an ellipse with axes formed by the diurnal reflectance SR1 and SR2 coordinates. It described previously studied soil line and the zone of canopy x canopy x ground interactions. The standard error of predicting {Delta}dw was 17%. We separately tested for plant height SR = f(h) or fresh weight SR = f(fw) using SR, and for dw as a function of normalized difference vegetation index NDVI = f(dw). SR = f(dw) produced superior results. Potential applications include noninvasive prediction of other biophysical plant properties in a single or in hyperspectral bidirectional reflectance in agronomic and ecological remote sensing.

Abbreviations: B, biophysical parameter for Canopy 1 (B1) or Canopy 2 (B2) • c1, clipped once • c2, clipped twice • cg, canopy cut to ground level • cstdv, conditional standard deviation • dw, biomass dry weight • E(|) conditional expectation • fw, biomass fresh weight • h, canopy height • ir, index of relationship • LAI, leaf area index • MSS4–MSS7, four Landsat wavelength bands • n1, Natural Canopy Number 1 • NDVI, normalized difference vegetation index • NIR, near-infrared reflectance • RI, reflectance index • SR, simple ratio • stdv, unconditional (commonly used) standard deviation • T, normalized local time • {Delta}dw, differential biomass dry weight • {theta}s, solar zenith (or elevation) angle • {theta}v, viewing zenith angle • {rho}, canopy reflectance ({Gamma}grass/{Gamma}Labsphere) • {varphi}s, solar azimuth angle • {varphi}v, viewing azimuth angle • {Gamma}, canopy radiance (W m–2 sr–1) • {sum}pxq, variance-covariance matrix p rows by q columns


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 
REMOTE SENSING DATA are commonly used in net primary production and albedo estimations and in determining canopy structure characteristics (Irons et al., 1988, Roughgarden et al., 1991). These data are most extensively used in agricultural modeling of crop yield and development (Moulin et al., 1998). Ground-level radiometric observations of vegetative surface are needed to interpret remote sensing data (Walter-Shea and Beal, 1990). The variables used in modeling plant canopy interaction with solar radiation include a biophysical parameter of interest (B) and reflectance index (RI). The first describes plant biology (e.g., cover, height, biomass, and leaf area index, LAI). B is typically a defined component of a natural plant canopy and must be distinguished from other background components that include thatch and the soil. The second variable of the model, RI, is typically defined in the near infrared NIR and red wavelengths and the index describes the radiometric response of the canopy as a whole. The ratio of near infrared NIR to red was first defined by Jordan (1969), and because of its simplicity, Middleton (1991) used and named it simple ratio (SR = NIR/red). The NDVI = (NIR – red)/(NIR + red) originally defined by Deering (1978) is widely used, e.g., Shanahan et al. (2001). Wiegand et al. (1979) calculated several reflectance indices derived from the four Landsat wavelength bands MSS4 (504–597 nm), MSS5 (598–700 nm), MSS6 (702–798 nm), and MSS7 (808–1042 nm). Each emphasized one or more aspects of reflectance interaction with vegetation and soil in the red and NIR wavebands. The practical utility of reflectance indices have been documented, but the inaccuracy in deciphering the interdependency of RI and biophysical property remains.

Experiments using biophysical plant properties (B) and RI's are commonly performed to develop an empirical relationship of RI as a function of B, that is RI = f(B). There are several unresolved problems in developing and in using the relationship RI = f(B) (Wiegand et al., 1979; Middleton, 1991; Maas, 1997, 1998). The problems in the development of the function include the difficulty of isolating different sources of {rho} (e.g., the thatch and the soil response from the green part) and in dealing with diurnal response because of changing sun angles and asymmetry in canopy geometry (Sandmeier et al., 1999). The problem in applying the function RI = f(B) is the difficulty in inverting the function to solve for biophysical property, B = g(RI). The function g is a mathematical inverse of f, g = f–1. In other words, in the application of the relationship RI = f(B), we want to insert a new observed RI, obtained from a similar vegetation type, into B = g(RI) to estimate a new biophysical parameter. The accuracy of the inversion determines the accuracy of the estimate which can be recovered from the model. More research is needed in this area.

The inversion of a reflectance function is complicated by interdependence of reflectance and plant biophysical property at different viewing angles (Middleton, 1991, 1992; Sandmeier et al., 1999). This is further complicated by nonuniformity of plant biophysical property such as plant height and density and by the interaction of reflectance from the exposed bare ground and dead vegetation below the canopy. For example, the two-dimensional distribution of trajectories of vegetation and soil states in the red–NIR coordinate system was plotted by Jensen (2006). The plot displayed a triangular region whose base defined a soil line oriented along the diagonal in red–NIR coordinate. The wet soil was nearer to the origin and the dry soil at the opposite end of the soil line. The two sides of the region joined together along the NIR axis. The greater the amount of green vegetation, the greater the reflectance in NIR, and the smaller in the red waveband. The trajectories that formed the two sides of the triangular region merged at a point to define high canopy biomass (Jensen, 2006). In other words, a model must be capable to solve for the biophysical plant property in terms of red and NIR observations anywhere in this region. In addition, the model must include the consideration of diurnal change in reflectance and its interactions.

Middleton (1991) studied Tallgrass Prairie and plotted one-dimensional sun angle-red (or NIR) trajectories. Each trajectory defined a diurnal biophysical plant property for various treatments including natural, burned, and grazed vegetation states. After examining the trajectories, Middleton (1991) concluded that trajectory shapes "could not be ascertained a priori due to complexity of the illumination effects on the vegetation and substrate." In other words, the solution of the inverse problem for various treatments was needed but was difficult to separate because of diurnal {rho} interactions in red and NIR wavebands. A successful analysis of this problem can have a wide range of application in agronomic and ecological remote sensing.

Our objective in this study was to propose mathematical and experimental procedures that address the problem of inversion and {rho} interaction, together, for any vegetation type. Specifically, we used grass as an example of a vegetation canopy (Walter-Shea and Beal, 1990, Middleton, 1991), dw as an example of biophysical parameter (Middleton, 1991), the SR = NIR/red (Middleton, 1991), and in the MSS wavebands (Wiegand et al., 1979) defined SR = MSS7/MSS5 as an example of a RI. The conditional probability served as the model (Shirazi et al., 2001) for separating diurnal reflectance interactions and model inversion.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 
Bidirectional Reflectance
The reflectance from a non-Lambertian surface such as plant canopy or its parts is bidirectional. This means that the reflectance is dependent on the relative geometry of incident and viewing angles (Daughtry, 1990; Walter-Shea and Beal, 1990; Middleton, 1991; Starks et al., 1991). The reflectance from a perfect diffusing surface (i.e., Lambertian) is independent of these directions. Such surfaces may be called isotropic. This is in contrast with the anisotropic vegetation surface (Middleton, 1992; Deering et al., 1992) which is direction dependent. There are two spherical coordinates that describe each of these directions, and they are commonly designated by {theta} for the zenith (or elevation) angle and {varphi} for the azimuth angle (Irons et al., 1988). Each of these angles is subscripted to identify the pair with the solar direction ({theta}s) and solar azimuth angle ({varphi}s) or the viewing zenith angle {theta}v and viewing azimuth angle ({varphi}v). Because reflectance is bidirectional, the RI is also dependent on the angles, such as RI = f({theta}s, {varphi}s, {theta}v, {varphi}v, B). Although we observed RI distribution with changing diurnal sun angles, for the sake of convenience we simply wrote RI = f(B) without the angles but we developed the function as a diurnal relationship.

We used an Exotech Inc. (Pompano Beach, FL) Model 100 BX radiometer with Multi Spectral Scanner (MSS4, 504–597 nm, MSS5, 598–700 nm, MSS6, 702–798, nm, MSS7, 808–1042 nm) with a 15° field of view lens. The four waveband tubes were parallel but not coincident (Jackson et al., 1980). An Omnidata International, Inc., (Logan, UT) Model 516C-32-A Polycorder was used to collect radiance data. The radiometer mounted on a tripod 55–70 cm above ground level to view the center of each patch and thus limit the field of view to the patch being studied. The viewing angle was fixed at nadir ({theta}v, {varphi}v = 0). The tripod legs were occasionally repositioned to avoid shading the patch.

Each data line in the Polycorder involved patch identification, local time, and date of observation, along with four sensed voltages corresponding to four MSS wavebands. Voltage readings from the MSS bands were converted to radiance values in W m–2 sr–1 using factory provided calibration tables. Clear sky conditions were monitored visually. Any sign of cloud cover, haze, jet trails, and so forth were noted and the corresponding radiance records were not used in the analysis. Clear sky conditions were intermittent, making observation times and sample numbers highly irregular. Consequently, radiance records for a particular canopy were patchy and contained segments from more than 1 d.

The sun zenith and azimuth angles ({theta}s, {varphi}s) were obtained from the local time (24-h day) and geographic coordinates using a calculation procedure developed by Walraven (1987). Other calculations were the times of sunrise, noon, sunset, and the daylength (= sunset – sunrise) for the day of measurement. We calculated a defined normalized local time, T = (local time – sunrise)/(daylength), and used it as an independent temporal scalar in our data analysis. For example, on the 60th day of the year in Corvallis, OR, sunrise = 7.1988 h, sunset = 17.6018 h, and daylength = (17.6018 – 7.1988) = 10.4030 h. Therefore, on this day at local time = 12.3950 h, the normalized time T = 0.5, that is, the solar noon in Corvallis. As an example of seasonal change during the test, sunrise = 5.0485 h and daylength = 14.1349 h in mid-May in Corvallis. Because daylength changes seasonally, the normalized time T is both a seasonal and diurnal scalar at any location.

Radiance {Gamma} as a Function of Normalized Time {Gamma}(T)
Using the radiometer setup described above, we measured radiance {Gamma} over a calibrated Labsphere 46- by 46-cm BaSO4 plate diurnally and in different seasons (dates). We fitted diurnally and seasonally observed radiance over reference plate to a fourth-order polynomial using the normalized time T = (local time – sunrise)/(daylength). Jackson et al. (1992) fitted similar data with fourth-order polynomial using azimuth angle. Figure 1 displays an example of our approach, fitting radiance {Gamma}(T) for the MSS7 waveband to a fourth-order polynomial using normalized time T. The data were recorded in Corvallis, OR, diurnally and during different seasons. If plotted against the absolute local time, each daily radiance curve begins at a different sunrise time and spans a wider length of the time axis in long summer days than in winter. The midday radiance peaks for different days of the year do not coincide at a single point along the local time axis. When plotted against scale time T (Fig. 1), the normalized curves for varying local sunrise and daylength collapsed together within an error band.


Figure 1
View larger version (25K):
[in this window]
[in a new window]

  		Fig. 1. Seasonal and diurnal radiance in Landsat MSS7 waveband recorded nadir over BaSO4 plate in Corvallis, OR, and fitted to a fourth-order polynomial in the normalized local time, T. The normalization, which adjusted local time for varying sunrise and daylength, helped to collapse seasonal data to within an error band shown.

 
Similar fits were obtained for all the four MSS wavebands. These are shown without the data in Fig. 2 . Notice that nested curves with different peaks at T = 0.5 are observed for different wavebands in Fig. 2, but for a fixed waveband (Fig. 1), the data collapse together within an error margin. We assumed from the small standard errors for these fits that seasonal and diurnal scaling by T was reasonably successful. In other words, the radiance calculated at any time T can be used to normalize a radiance observation over grass canopy at equal time T and equal waveband. The advantage of not having to make simultaneous observations over grass canopy and reference plate is obvious.


Figure 2
View larger version (24K):
[in this window]
[in a new window]

  		Fig. 2. Seasonal and diurnal radiance in four Landsat wavebands each fitted to observed radiance (data points not shown) over BaSO4 plate in Corvallis Oregon using a fourth-order polynomial in normalized local time T.

 
Grass Canopy Treatments Used in This Study
The grass plot was a 9- by 28-m experimental low-maintenance grass field at the Environmental Protection Agency's Environmental Research Laboratory, Corvallis, OR. The field was a mix of bluegrass (Poa trivialis L.), clover (Trifolium repens L.), and western weeds which included thistle [Cirsium arvense (L.) Scop.], nightshade (Solanum nigrum L.), yarrow (Achillea millefolium L.), milkvetch (Astragalus Spp.), and pineappleweed (Matricaria discoidea DC.). We did not measure the percentage composition of the species. The ground-level composition was either thatch or bare soil. The grass was not watered, grew naturally, and was mowed once after a peak growth in spring. Nearby buildings were located north of the plot and did not shade any part of the grass that we studied.

We collected radiance over grass patches during May 1994 and March 1995. Distinct 35- by 35-cm test patches were studied by recording the canopy height (h) and then the canopy radiance ({Gamma}) at 5-min intervals from sunrise to sunset. Some grass patches were studied in consecutive days. For these, the intact (undisturbed) canopy was studied during the first day. A grass layer of constant thickness was clipped from the top of the patch and the clipped canopy was studied the next day. A second clipping from the top of the same patch was continued to carry out the study of a twice-clipped patch in subsequent days. Finally, the ground-level radiance record of the exposed thatch and the soil was obtained after a full harvest. Grass layers were clipped to uniform heights with a manual sheep shear, taking care to collect all loose grass blades. The weight of the freshly clipped vegetation was recorded, oven-dried at 85°C, and its dry biomass was used in the statistical model.

Patches have one of the four canopy descriptions: a natural canopy with the top growth intact designated by a small scale letter (n) and a number (e.g., n1, n7, to distinguish it from a separate natural grass patch), a canopy after the top grass has been clipped once (c1, where the letter c stands for clipped), the same canopy after a second clipping (c2), and the bare ground after a total harvest of the grass (cg, canopy cut to ground level). The data for the ground-level radiance were gathered over multiple bare grounds as various patches in the plot were harvested. For example, in 1995 the bare-ground radiance was recorded from 7 March through 3 April over four neighboring bare-ground patches in the plot.

Summarized radiance observations in this study are presented in Table 1. The first column in the table describes the patches, their codes, and the date of the record for the patch. The next three columns are the observed canopy properties for the patch. The fourth column is the number of data lines of clear sky diurnal radiance observations recovered from the Polycorder. The next two columns are the times of first and the last recovered lines for the same patch. The last two columns in Table 1 list the standard errors of fitting grass radiance observations over the patch in MSS5 and MSS7 wavebands to separate cubic polynomials in the normalized time T. Figure 3 shows an example of radiance fitted to four separate cubic polynomials in MSS wavebands using the first patch in Table 1.


View this table:
[in this window]
[in a new window]

  		Table 1. Grass patches are identified by a letter and number codes as seven natural patches (n1, n2, ..., n7), patches where the top is clipped once (c1), clipped twice (c2), and (cg) harvested to ground level, and by the date of the test.{dagger}

 

Figure 3
View larger version (28K):
[in this window]
[in a new window]

  		Fig. 3. An example of diurnal radiance in four Landsat wavebands recorded nadir over a patch of grass in March 1995 with naturally grown canopy code n1 and dry wt. of 229.3 gm–2. The radiance in each waveband is fitted to a cubic polynomial in normalized local time T (Table 1, Line 1).

 
Reflectance {rho} as a Function of Normalized Time {rho}(T)
The conversion of {Gamma} to reflectance was made by evaluating the fourth-order polynomial (Labsphere data) and the cubic polynomial (canopy data) for the same MSS waveband at equal normalized local time T to produce the ratio: reflectance = {Gamma}/(Labsphere radiance). Examples of reflectance plotted for four MSS wavebands and four different canopy descriptions (n1, c1, c2, and cg, Table 1) are shown in Fig. 4 . Reflectance in MSS4 and MSS5 wavebands differ very little. But only MSS5 waveband was used in our study. Although we collected data beginning each day at sunrise and ending at sunset, clear sky conditions did not prevail during early and late hours of the day. After examining our records, we settled on a start (Tstart = 0.3) and end time (Tend = 0.8) that was common to all records. This explains why in Fig. 4 the reflectance was evaluated for 0.3 ≤ T ≤ 0.8. The evaluations in this figure (and in the remaining diurnal calculations) were performed at 15 equal increments of normalized time T between these limits (i.e., {Delta}T = 0.0333).


Figure 4
View larger version (26K):
[in this window]
[in a new window]

  		Fig. 4. Diurnal reflectance in four Landsat wavebands showing examples of naturally grown grass canopy with top intact (code = n1), the top once clipped (c1), twice clipped (c2), and canopy fully harvested to the bare ground (cg). Canopy codes, dates and dry weights match the listings in Table 1.

 
Curves similar to those shown in Fig. 4 were reported by Deering et al. (1992) over prairie grassland. They obtained their results using a constant zenith angle, but changing diurnal view angles. By contrast, we fixed the view angle and recorded the changing solar angle. The two procedures are equivalent because of Helmholtz Reciprocity Theorem where the locations of the incident source and the viewing instrument can be interchanged without affecting the distribution of bidirectional reflectance (Silva, 1978). Experimentally, there is a clear advantage in this alternative approach which we used.

Linking Biophysical Parameter and Reflectance Index
To separate reflectance interactions between the canopy and the background, we developed two distinct empirical diurnal records of reflectance, RI1 and RI2, where the subscripts refer to two separate canopies. Then, the difference between RI1 and RI2 was modeled to obtain a statistical estimate of the difference between biophysical properties (B1 – B2). Therefore, the first record was developed over a canopy expected to include multiple sources of reflectance interactions. The second was restricted, in that B1 > B2. When B2 = 0, diurnal canopy interaction was with the bare ground. When B2 > 0, the statistical prediction included difference between the two geometries over the diurnal sun angles. In other words, by our choice of a predefined background record, RI2 = f(B2), we intended to use the model to separate the background reflectance interaction from RI1 = f(B1).

More specifically, a diurnal record was developed of the SR1 over a patch of grass with a standing biomass dw1 to obtain the first function SR1 = f(dw1). This was repeated over a second patch of grass with a standing biomass dw2, where dw1 > dw2 (or dw1 >> dw2) to obtain the second function SR2 = f(dw2). The first patch could be intact, that is, a naturally grown grass canopy or it could be grazed because we intentionally clipped a thin layer from its top. The second patch could be natural, clipped, or totally harvested, leaving the roots intact. The dw2 in the harvested situation was assumed zero.

The functions SR1 and SR2 were combined in a conditional probability model (Shirazi et al., 2001) that described the mean difference of the standing dry biomass weights dw1 and dw2, that is, {Delta}dw = dw1 – dw2. The model calculated the mean difference {Delta}dw from the statistical difference of the diurnal functions SR1 and SR2. This difference was defined in a precise way by splitting the observational data as commonly done in a conditional probability calculation. The formal statement for the model that calculated the mean or expectation E of {Delta}dw was E({Delta}dw | SR1, SR2). The symbol | indicates that the mean value {Delta}dw is dependent on the joint occurrence of the diurnal reflectance SR1, SR2. A geometric explanation of the approach is gained by the fact that a constant probability relationship between the three variables {Delta}dw, SR1, and SR2 is defined on the surface of a solid ellipsoid. The projection of that surface is a conditional probability defined by an ellipse whose axes are SR1 and SR2. Thus, {Delta}dw may be defined by a polynomial of the form {Delta}dw = a + bSR1 + cSR2 + cSR1SR2 + eSR12 + fSR22, where a, b, c, d, e, and f are constants. Therefore, the conditional probability serves precisely the desired inverse function.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 APPENDIX
 REFERENCES
 
Model Development
Using the first five patches in Table 1, we combined two patches at a time (each having biomass dw1 > dw2, {Delta}dw = dw1 – dw2). Each coupled patch defined a treatment. Treatments were used to train the model. There are only experimental and practical limits to the number of treatments that may be used to train a model. For example, any number of coupled patches with different {Delta}dw values may be used in training the model. We previously showed that {Delta}dw may be expressed in terms of a binomial number of combinations of SR1 and SR2 (see above). We thus assumed that the number of treatments for each diurnal record must be no less than six, which is the number of unknowns that defined {Delta}dw in the previously defined model.

We coupled a patch in Table 1 repeatedly with a second patch to form treatments. For example, in Lines 1 through 4 (Table 2), the same patch was coupled: with an intact neighboring patch (Line 1), with itself after clipping one layer from its top (Line 2), with itself again after a second clipping (Line 3), and with a bare ground (line 4). We formed a total of 10 treatments (coupled patches) shown listed in the top part of Table 2. In all cases the dw of the first patch was larger than the second. It mattered that we trained the model using a range of treatments. Thus, one coupled pair was natural (n2n1), two naturals coupled with ground (n2cg, n1cg), four naturals coupled with clipped (n2c1, n2c2, n1c1, n1c2), one with both patches clipped (c1c2), and two clipped patches coupled with ground (c1cg, c2cg). The independent dataset which was used to test the model (Table 2, bottom part) had seven natural canopies coupled each with a natural patch, four naturals coupled with ground, and two naturals coupled with clipped patches (n5c1, n4c1). In other words, more than one-half of test dataset was paired natural patches compared with 10% of the model dataset.


View this table:
[in this window]
[in a new window]

  		Table 2. By combining two grass patches with unequal dry weights, dw1 > dw2, from Table 1, 10 combinations are produced as shown in the top part of the table to produce a model using either Eq. [A1] or Eq. [A2] and 13 separate combinations listed in the lower part of the table are produced in testing the models. The differential dry weight ({Delta}dw, gm–2) listed in the last column is marked with patch codes from Table 1. For example, coupling the two natural patches listed in the first line produced the rounded differential weight 97 gm–2 which is identified by the double code n2n1.

 
Recall that the SR was defined by dividing reflectance in the MSS7 waveband by the reflectance in the MSS5 waveband. Thus, for each coupled set of patches and a constant {Delta}dw, we developed a diurnal record of reflectance (SR1) from the first patch and a separate record (SR2) from the second patch (Table 2). The record from the first patch presumably contained reflectance interactions from multiple sources. The record from the second patch was assumed to contain a predefined background which must be separated from the first. The diurnal interactions of SR1 and SR2 for the 10 pairs of grass patches are plotted in Fig. 5 . Each produced a distinct path which we named a trajectory. The trajectories are numbered and coded with line numbers and codes in Table 2. For example, the top rightmost trajectory represents Line 1 in the table and is identified as 97/n2n1. The trajectories changed direction as the sun angle changed diurnally. The sun angles were not explicitly plotted but their values have been represented at 16 points along each trajectory (0.3 ≤ T ≤ 0.8). Trajectory numbers are placed near the start of the morning lobe at T = 0.3 while a triangular marker identifies the 16th point at the end of the afternoon lobe at T = 0.8.


Figure 5
View larger version (17K):
[in this window]
[in a new window]

  		Fig. 5. Trajectories of diurnal canopy reflectance interactions in SR1 and SR2 obtained by using paired treatments of grass patches: two natural patches (n2n1), natural with clipped (n2c1, n1c2), two clipped (c1c2), and canopies in combination with bare ground (c2cg, n1cg, and n2cg) listing the differential weight ({Delta}dw) in round numbers (Table 2). The triangular zone depicts increasing canopy x ground (left-to-right) interactions in the four bottom row trajectories, increasing canopy x canopy (bottom-to-top) interactions in the four trajectories on right side of the zone, and mixed interactions in the interior and diagonal side.

 
The diurnal change in the direction of a trajectory described the combined interaction of diurnal reflectance from paired patches. For example, when SR1 from the first patch started to increase with increasing sun angle, the direction of trajectory-change in the morning lobe began from left-to-right along SR1. A series of complex directional changes may take place along each trajectory because of the combined changes in SR1 and SR2 with increasing (morning lobe) and decreasing (afternoon lobe) sun angles. Figure 5 displays a wide range of directional changes along each trajectory and among the 10 trajectories that define a triangular zone of reflectance interactions in SR1 and SR2 coordinates. Trajectories, 1 to 4 along the SR2 axis define the right side of the triangle; 1,5,8 and 10 define the diagonal; and 10, 9, 7, and 4 along the SR1 axis define the lower side.

Focusing first on the zone as a whole, one observes that trajectories in the base of the triangle describe canopy interactions with the bare ground, and the effect of increasing canopy biomass from the left to the right. The trajectories on the right hand side describe diurnal interactions in decreasing biomass difference ({Delta}dw) between patches. The trajectories in the interior and along the diagonal side of the triangular zone represent mixed canopy x canopy x ground interactions. Although not identical, the zone described here recalls Jensen (2006), who plotted a triangular zone of interaction of the soil and vegetation states in the Red-NIR coordinates. Because we combined Red-NIR axes into a single reflectance ratio (SR) for each patch, the triangular zone in Fig. 5 is a differential biophysical property. Therefore, the base of this triangle represents the extension of Jensen's soil line and trajectories along parallel lines with the base representing different degrees of canopy x canopy interactions.

Next, focus on individual trajectory interactions within the triangular zone. Starting with the shape of the first trajectory in the base of the triangle (112/c2cg), we observe it is nearly flat and the individual trajectory points are evenly spaced. This indicates a uniform dominance of diurnal reflectance from Patch 1 over Patch 2. As the canopy biomass increases in 174/c1cg, the trajectory points in the morning lobe become more crowded than the points in the afternoon lobe. The crowding advances toward the afternoon lobe with increasing biomass difference in 229/n1cg until the trajectory in 326/n2cg completely reverses direction. This trend remains generally unchanged in all canopy x canopy interactions in four trajectories on the right side of the triangle. Then after, along the interior and diagonal side of the triangle, trajectories virtually flip half-way over (55/n1c1 and 118/n1c2), aligning with the SR2 axis. Direction change continues in 62/c1c2 until it completely reverses again in 112/c2cg. We believe that Jensen (2006) described a similar interaction in the Red-NIR coordinate and named it "migration of a pixel." What we just described was the mechanism of interaction that quantifies this migration in our system.

All the relationships described relative to Fig. 5 are in consideration using a conditional probability to estimate a specific mean {Delta}dw from correlating diurnal trajectories. Although the calculation of conditional probability is well documented in the literature (Johnson and Wichern, 1982), we thought it was helpful to interpret a few steps in these calculations in terms of the specific problem at hand. For example, the top rightmost trajectory, 97/n2n1, is a statement of mean differential dw {Delta}dw = 97 gm–2 constrained by the joint reflectance values 31.3 < SR1 < 33.3 and 24.7 < SR2 < 25.3. In other words, the conditional expectation of {Delta}dw is E({Delta}dw = 97 gm–2| 31.3 < SR1 < 33.3 and 24.7 < SR2 < 25.3), where in the model the values of the functions in SR1 and SR2 are known as conditions. It is true that the observed {Delta}dw is a constant, but {Delta}dw has a distribution along each diurnal trajectory when using the model. The average for a particular trajectory is determined by the model from considering relationship with all trajectories in the triangular zone. The manner that this calculation is accomplished with model dataset is described below (and more formally in the Appendix).

The model dataset provided the mean, standard deviation (stdv), and the cross-correlation among {Delta}dw, SR1, and SR2 (Table 3). Notice that there is nonzero cross-correlation among all variables, indicating first-order (linear) statistical interactions. The stdv's calculated in this table are one-dimensional (nonconditional). One can reduce variability for any variance in the table by partitioning it from the rest. For example, the four elements of the matrix in the lower right hand corner of the table exclude the differential weight ({Delta}dw), and partitioning {Delta}dw this way (see Eq. A1), obtains the conditional variance, var({Delta}dw | SR1, SR2). The magnitude of the conditional stdv obtained by partitioning is cstdv({Delta}dw) = 15.2 gm–2 against the one-dimensional stdv({Delta}dw) = 80.3 gm–2. This reduction in variance (and increased accuracy) was obtained precisely because of the nonzero correlations in Table 3. Using similar calculations, Shirazi et al. (2001) defined an index of relationship, ir = 100 (1 – cstdv/stdv), to describe the degree of explanation of experimental variance by a conditional probability model. This index was 81% for the model dataset.


View this table:
[in this window]
[in a new window]

  		Table 3. Statistical summary of model variables listing the mean, standard deviation (stdv), cross correlation (symmetric) matrix of the differential dry weight ({Delta}dw), and simple ratios SR1 and SR2.

 
The calculation of conditional mean (by partitioning) follows a similar path but requires inserting specific diurnal trajectory values of SR1 or SR2 into each calculation (e.g., similar to the ellipse equation at the end of MATERIALS AND METHODS). Unlike the variance model (Eq. [A1]), which describes the entire dataset (i.e., the triangular zone in Fig. 5) the conditional mean model (Eq. [A2], Appendix) describes instantaneous positions along the trajectory shown for each record in Fig. 5. There were 16 data points that described each trajectory. A separate differential weight, {Delta}dw, was predicted by Eq. [A2] (Appendix) for each data point along the trajectory. This distribution of the differential weight must be integrated (or simply averaged) over the entire trajectory to completely describe the diurnal response. That the temporal averaging, over the diurnal trajectories, is indeed reflective of a trajectory pattern is discussed below.

The results of averaging over 10 diurnal trajectories are presented in Fig. 6 . The points in this plot were numbered to link them with line numbers in the upper part of Table 2. The symbols represent the mean trajectory values and the bars through the symbols are the stdv's. Focusing first on the whole data set, the overall standard error of prediction for the model was 15 g m–2 and is represented by the solid line in Fig. 6. Next, we focus on individual trajectories and use the model for interpolation or self-prediction within the triangular zone of its own data. The smallest mean and stdv errors of prediction in Number 7 associate with collapsing trajectory points (Fig. 5). Five trajectory points in the morning lobe of Trajectory 7 and eight points in the afternoon lobe collapse atop each other, indicating nearly equal interaction of SR1 and SR2 at these points. Number 4, with two trajectory points collapsed from the morning lobe atop six from the afternoon lobe, has a small mean error and a relatively small stdv. Trajectory 1 is similar but the mean response is apparently dominated by SR1, placing it in the same class as numbers 10 and 9, which generally parallel that axis. Trajectories 5 and 6 are generally dominated by SR2 axis and trajectories 2 and 3 are dominated by SR1 in their morning lobes and by SR2 in their afternoon lobes. Finally, trajectory 8 is evenly but jointly dominated by SR1 and SR2, leading to small mean and stdv errors.


Figure 6
View larger version (14K):
[in this window]
[in a new window]

  		Fig. 6. The error of recovering modeled dry weight (i.e., self-prediction of dry weight difference) for the whole triangular region in Fig. 5 (solid regression line) and the mean for each trajectory (symbol) and standard deviation relative to the mean (vertical lines drawn through the symbol). Data points are linked with line numbers in the top part of Table 2. Trajectories nearest the solid line are associated with strong and joint diurnal reflectance of paired patches.

 
In summary, the patterns that are constant functions of SR1 and SR2 (7 and 4, Fig. 5) or proportional with SR1 and SR2 (8) are associated with the smallest errors in predicting {Delta}dw. Other trajectory patterns generally produce larger errors because of weaker diurnal relationships with both SR1 and SR2. These include trajectories that vary with SR1 in the morning lobe and with SR2 in the afternoon lobe (2, 3) and those that vary mainly either with SR1 (1, 9, 10) or with SR2 (5, 6). In other words, the diurnal prediction of {Delta}dw is two-dimensional except when the canopy alone or the background alone dominate the response.

Model Testing
We tested the above-described model, with no further modification, by using a separate dataset obtained independently of the model dataset. When discussing Fig. 5, we explained the ways that diurnal difference in sun illumination and biophysical difference in plant properties combine to form a trajectory for the model dataset. We also explained that the test dataset (bottom parts of Tables 1 and 2) was seasonally and structurally different from the model dataset. Some related quantitative differences were as follows: daylength was 3 to 4 h longer, and the grass was drier in May 1994 when the test dataset collected compared with March 1995 for the model dataset. The mean fw–dw ratio was 4.6 (SD = 1.2) in May 1994 compared with 6.8 (SD = 0.9) in March 1995. The test dataset produced the 13 trajectories shown in Fig. 7 which are substantially different, if not more uniform, than trajectories for the model. We do not know the exact ways that the above differences between the two datasets combine to produce such dissimilar trajectories. More research is needed to better explain trajectory forms relative to known additional biophysical parameters such as the dominance of natural x natural canopy structure (treatments) in the test dataset compared with the model dataset.


Figure 7
View larger version (21K):
[in this window]
[in a new window]

  		Fig. 7. Trajectories of diurnal canopy reflectance interactions obtained mainly by pairing natural canopies (7 of total 13 test dataset, bottom part of Table 2), four natural canopies paired with bare ground, and two with clipped canopies (n5c1, n4c1). The contrasting trajectories here and in Fig. 5 represent seasonal differences and dominance of natural-natural canopy treatments.

 
To obtain numerical validation, we inserted trajectory coordinates from Fig. 7 into Eq. [A2] to predict a differential mean dw, {Delta}dw, for the test dataset and compared it with observed {Delta}dw. The results are plotted in Fig. 8 . The points in this plot were numbered to link them with the line numbers in Table 2 (bottom part). The standard error of prediction for verifying the model performance against an independent, albeit slightly different, dataset was 26 gm–2. The relative error was calculated by dividing the overall standard error of the test by the overall average differential weight to obtain 100(26/145) = 17%. This compares with 100(15/154) = 9.7% error in the model itself (Fig. 6). Recall that seven patches of test dataset had paired natural patches compared with only one in the model. In other words, the model was used to predict greater diversity of treatments compared with its training dataset. Thus, a 17% error from applying the model in an extrapolation mode appeared reasonable to us.


Figure 8
View larger version (16K):
[in this window]
[in a new window]

  		Fig. 8. The error of predicting dry weight of test dataset by using all trajectories in Fig. 7 (solid regression line) and the mean for each trajectory (symbol) and standard deviation relative to the mean (vertical lines drawn through the symbol). Data points are linked with line numbers in the bottom part of Table 2.

 
Choices of B and RI in a Model
It is useful to have an index that calculates the effectiveness of different choices when pairing biophysical parameter and reflectance (B and RI). Such an index can be used to compare various mixes of variables on similar basis. Recall that we calculated the ir for the paired model dw and SR. Because ir described the strength of association between pairs, it offered a uniform but robust measure of comparison among pairing alternatives. Several examples are listed in Table 4 to describe this application using our model dataset (top parts of Tables 1 and 2). For each example, we calculated the one-dimensional and conditional stdv's to obtain the index ir. Although the choice, dw and SR, was superior to other combinations listed in the table, the combinations dw and NDVI and h and SR ranked among the lowest, with fw and SR ranking second after dw and SR.


View this table:
[in this window]
[in a new window]

  		Table 4. Index of relationship ir = 100(1 – cstdv/stdv) calculated for three vegetation indices: dry weight (dw), fresh weight (fw), and canopy height (h), and two reflectance selections: simple ratio (SR = MSS7/MSS5) and normalized difference vegetation index (NDVI) = (MSS7 – MSS5)/(MSS7 + MSS5).{dagger}

 

   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
APPENDIX
 REFERENCES
 
More research is needed to separate diurnal reflectance from the background and the green fraction in remote sensing studies. In this paper, a nondestructive field method is developed to predict biomass (and other canopy biophysical properties). Diurnal RI trajectories, from the nadir view angle with changing sun angles, were acquired at intervals between sunrise and sunset. A conditional probability model was trained to distinguish varying biomass difference, for example, between two paired plots on the basis of their trajectory differences. The goal of the proposed approach was to estimate biomass in a nondestructive manner across other plot pairs for a similar vegetation type (e.g., switchgrass, Panicum virgatum L.).

Using {Delta}dw as an example of plant biophysical property, model results showed that {Delta}dw can be predicted after separating background diurnal reflectance. However, when the background or the other diurnal reflectance dominates, the dominant reflectance may be used alone to predict the dw with a reduced number of variables.

Nondimensional time scalar collapsed diurnal and seasonal reflectance from reference plate by a polynomial fit for any waveband. When used at a fixed location in a radiometric study of plant canopy, such data eliminated the need for simultaneous radiance observations over the canopy and a reference plate. Further, it was operationally more convenient to fix the view angle in this study and collect diurnal radiance with changing sun angle than the reverse.

Although conditional probability was used for only three interdependent variables (i.e., SR1 and SR2 with only one biophysical parameter {Delta}dw or {Delta}fw or {Delta}h), more than one biophysical parameter could be used in the model to predict the joint occurrence of the mean parameters together. For example, two separate dw's may be used in place of differential dw or differential dw and differential h may be used together in this study and two representations of LAI (total and illuminated green, Middleton, 1991) could be used along with a RI. Likewise, additional reflectance indices may be used in the model, for example, to study hyperspectral bidirectional reflectance distribution over grass (Sandmeier et al., 1999). It is particularly helpful to expand our results to include LAI in these tests.


   APPENDIX
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
REFERENCES
 
We began with a dataset which we obtained from observing diurnal reflectance SR over separate grass patches having {Delta}dw = dw1 – dw2. We wanted to develop an inverse relationship which was solved for {Delta}dw. In other words, we observed the relationship (SR1, SR2) = f({Delta}dw), and we wanted to end up with the inverse relationship {Delta}dw = g(SR1, SR2). The first step was to calculate the variance–covariance matrix of the three interrelated variables {Delta}dw, SR1, and SR2 listed in Table A1. The left hand side of the table was a symmetric three-by-three matrix, which we designated by the symbol {sum}3x3. The diagonal elements were the one-dimensional variances of the three variables {Delta}dw, SR1, and SR2. The off-diagonal elements were a mixed covariance of two variables. On the right side of Table A1, we showed a four-way split of the original {sum}3x3 matrix. The first diagonal element, {sum}1x1, had a single member and was the one-dimensional variance in {Delta}dw. The second diagonal element, {sum}2x2, had four members, which totally excluded {Delta}dw. The off-diagonal elements were a one-by-two row-matrix, {sum}1x2, and a two-by-one column-matrix, {sum}2x1, each with mixed two-dimensional members. This manner of partitioning allowed calculation of conditional statistics {Delta}dw given that SR1, and SR2 have been diurnally sensed. This is formally written as E({Delta}dw|SR1, SR2) for conditional mean. The conditional variance was calculated from the difference between the one-dimensional variance, {sum}1x1, and any nonzero correlation present between {Delta}dw and SR1 or SR2 as follows (Johnson and Wichern, 1982):

Formula 1[A1]
The second term, {sum}1x2{sum}–12x2{sum}2x1, on the right hand side of Eq. [A1] defined the correction to be applied to the one-dimensional variance, {sum}1x1. Note that the variance–covariance matrix, {sum}2x2, which we previously described, appears inverted in the correction term in Eq. [A1]. This inverted matrix, {sum}–12x2, defined the solid ellipsoid model in the correction. This result was precisely what we sought out in the inversion g = f–1. Similar to Eq. [A1], the conditional probability model of the mean {Delta}dw corrected the one-dimensional mean differential weight µ1x1 by the magnitude of the cross-correlations previously described (Johnson and Wichern, 1982). Thus,

Formula 2[A2]
Note that this model required inserting specific diurnal values of SR1 or SR2.


View this table:
[in this window]
[in a new window]

  		Table A1. The variance–covariance (symmetric) matrix of model variables showing the partitioning strategy for calculating conditional variance ({Delta}dw |SR1, SR2).

 
Multidimensional biophysical parameters (e.g., dw, fw, cover, h, and LAI) and RI (e.g., hyperspectral diurnal reflectance) may be modeled with this procedure. Typically, the p by p variance–covariance matrix is first calculated then partitioned as described here in any combination. The goal is to produce a model of one biophysical parameter or a model of joint probability of several biophysical parameters predicted by reflectance and other known biophysical parameters.


   ACKNOWLEDGMENTS
 
Dr. Craig S. T. Daughtry USDA-ARS Hydrology and Remote Sensing Laboratory read the first draft of this manuscript and made valuable comments. The Associate Editor and the anonymous reviewers of Agronomy Journal made valuable comments and helped streamline the presentation. The U.S. Environmental Protection Agency through its Office of Research and Development funded the research described here. It has been subjected to Agency review and approved for publication.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX
 REFERENCES
 





This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Shirazi, M. A.
Right arrow Articles by Reporter, M.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Shirazi, M. A.
Right arrow Articles by Reporter, M.
Agricola
Right arrow Articles by Shirazi, M. A.
Right arrow Articles by Reporter, M.
Related Collections
Right arrow Inverse Procedures/Parameter Estimation
Right arrow Remote Sensing
Right arrow Other Integrated Systems

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Crop Science Vadose Zone Journal
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome