Non-linear Canopy Reflectance Model Inversion

MSc. Remote Sensing

P. Lewis & M. Disney

Remote Sensing Unit,

Dept. Geography

University College London,

26 Bedford way,

London WC1H 0AP

+44-71-387-7050 x 5557

e-mail:

January 14th 2006

1.INTRODUCTION

So far, we have been considering descriptions of a relationship between canopy reflectance (or a VI) and some parameterization of the surface (e.g. LAI, chlorophyll concentration). Such a relationship is generally termed a forward model.

If we wish to attempt to obtain information on the surface parameterization from observed RS data, we need to invert such a model. Model inversion is therefore of key importance to the remote sensing problem. In this lecture, we will discuss the issues involved in attempting an inversion, and have an introduction to the range of methods available for the task. Many of the issues and options raised are generic in nature and applicable to a wide range of scientific problems.

1.2 ISSUES IN NON-LINEAR INVERSION

1.2.1. Error Surfaces

Earlier in the course, we considered linear model inversion as part of our study of lumped parameter modelling. We will see that by phrasing the problem as a minimization of some error function (e.g. RMSE) we can solve for the inverse of a linear model by matrix inversion. In this case, the RMSE ‘error surface’ (in N dimensions for N model parameters) is a quadratic function. It has only a single (unconstrained) minimum. We can apply (linear) constraints through simple methods such as Lagrange multipliers.

1.2.2 Parameter transformations and bounding

More generally, especially when considering inverting physically-based canopy reflectance models, we do not have an analytical solution of this type, often because there is no analytical solution to the (non-linear) simultaneous equations obtained from the partial differentiation. Also, more generally, the error surface is more complex than a quadratic function – this is easy to visualize if one considers non-linear parameters which are known to threshold (i.e. the remote sensing signal, thence RMSE, is not sensitive to changes in the parameter beyond a certain point) such as LAI. For high LAI, the error surface becomes ‘flat’.

One way top get around this problem is by approximately linearising the parameter space. The following table, taken from Combal et al[1] (2002), shows a set of transformations for a range of optical canopy reflectance model variables. By transforming the physical variables to these terms, and setting robust limits to the terms, we can attempt a more reliable inversion, obtain a ‘simpler’ (closer to quadratic) error surface, and (partially) avoid saturation problems.

Table 1. Biophysical parameter transformations

Figure 3 shows a typical error surface (the model parameters are transformed LAI and soil brightness) for a non-linear model – the error surface for non-transformed LAI would show saturation at high values – in the transformed space, the global; minimum is more clearly accentuated.

In fact, provided there is sufficient information in the remote sensing signal (i.e. sensitivity to the model parameters), we find that in many cases for optical remote sensing, the error surface is relatively ‘well-behaved’. We will demonstrate the importance of this point later. Note that the same is often not the case for microwave signals – this is partially a sensitivity to a large number of model parameters, and partially a lack of information in the data received from satellite data.

With certain numerical (error) minimization routines, there may also be problems associated with finding a global minimum, rather than a local one. Situations such as inverting RS parameters can often be helped by assigning bounds to the domain in which we search for a solution. When using the Powell algorithm, this is achieved by adding a gradually-increasing penalty to the error function if we step outside defined bounds. This does sometimes cause some problems with the inversion getting trapped on a boundary, especially if the penalty increase is too sharp on the boundary.

Figure 3. Error surface from non-linear model

1.2.3 Weighting of the Error Function

Another issue to be considered at this point is how to phrase the error function. In the Method of Least Squares, one is supposed to weight each residual by the uncertainty associated with the observations. In this way, the RMSE is similar to a Z-score, in being a normalized function relating to the average association of the observations with the model values.

In the previous lecture, we weighted all samples equally. This is equivalent to assuming that the error (noise) associated with each measurement is equal. This means that, for example, if we were trying to invert chlorophyll concentration from visible and near infrared hyperspectral data, we would be weighting the residuals in the near infrared (where there is no sensitivity to the required parameter) as strongly as those at visible wavelengths.

If the aim of the inversion were a robust estimation of a single parameter, we could chose instead to weight the residuals by the sensitivity of the signal to that parameter (this applies to angular as well as spectral sampling). We might chose a partial derivative as the weighting term – the higher the rate of change of reflectance with respect to the parameter of interest, the higher the weighting, so, instead of using:

(1)

we could use:

(2)

Such an approach is suggested by Privette (1994)[2]. There are clearly issues with the implementation of such an approach – principally: (i) the derivative will vary as a function of the model parameters for a non-linear model (so some ‘generally appropriate’ mean value must be chosen); and (ii) weighting for multiple parameters is somewhat arbitrary (so one has to decide on a compromise weighting scheme, or invert multiple times for each model parameter). As a partial compromise, one most typically sees a ‘relative RMSE’ used in remote sensing applications:

(3)

rather than the ‘absolute’ RMSE shown in equation 1.

1.2.4 Using Additional Information

Whatever method we use for inversion, and however we phrase the error function, the potential for problems exists if any model parameters are coupled (i.e. a change in one parameter can be compensated for by change in another parameter without affecting the solution at all). One way to get around this is to use additional information in the inversion problem. The most typical example of this, for vegetation science applications, is to use some model of the expected temporal dynamics of a vegetation canopy in inversion. Moulin et al. (1998)[3] review methods for linking crop growth models to the remote sensing inversion problem. Methods generally rely on using a priori information on the expected crop type (e.g. obtained from a land cover classification) and ‘growing’ the crop model to provide an expectation of canopy state variables. The inversion problem is then much more limited in scope, and one can apply expectations e.g. of coupling between clumping, LAI, and leaf angle distributions for a given crop type. Crop or other vegetation growth models, however, tend to require a significant number of ‘internal’ parameters describing how a particular variety will grow. Most models are based on considerations of ‘integrated thermal time’ (sum of degree days since planting), so knowledge of the local daily mean temperature is a primary requirement for running such models. In addition, they may have limits applied to the light environment (so that the models can be used at different latitudes). Some models consider water and nutrient limiting factors, although there is much disagreement between models as to the impacts of these. A key use then, of a canopy growth model, is to provide expectation of growth under non-limiting conditions – any large disagreement between model and measure may then be used to identify the impact of some additional limit (water, nutrients etc.).

1.2.6 Scaling

The models we build of canopy reflectance or scattering will typically be defined assuming operation over a single cover type – indeed, many models assume that the canopy is of horizontal infinite extent. In reality, this may not be the case, particularly for moderate or coarse resolution sensors. The advantage of optical moderate resolution sensors for monitoring crop dynamics is the high temporal coverage that can be achieved, even taking into account data lost due to cloud cover (which may be typically as high as 75% loss in Northern Europe). The main disadvantage is that pixel resolutions of 100s of metres to several km will tend to observe more than one field at a time, particularly given North European field sizes of typically a few Ha. High spatial resolution optical data are available relatively rarely for cloudy environments, and cannot be guaranteed to provide suitable temporal monitoring for crops. Current microwave sensors tend to have a much lower information content, but higher spatial and temporal resolution.

If we wish to invert canopy parameters from moderate resolution sensors then, we must be aware of (or ideally take account of) scaling effects. To gain an insight into scaling effects, consider a mosaic of crop types viewed within a single moderate resolution pixel (figure 4). The area comprises 20% of LAI 0, 40% LAI 4, 40% LAI 1. The ‘real’ total value of LAI is thus: 0.2x0+0.4x4+0.4x1=2.0. We assume a model of reflectance as:

(4)

We now consider two typical cases for leaf and soil properties:

visible:

NIR

For the visible reflectance, the canopy reflectance over the pixel is 0.15 and 0.60 for the NIR. If we assume the model in equation 4 to hold, this equates to an LAI of 1.4. The ‘inversion’ underestimates the ‘correct’ (scaled) answer of 2.0. In addition, different results may occur for different wavebands (they do not here as we have only a single scattering model). It is worth noting that some parameters will scale more linearly than others – examples of those that do are albedo and canopy cover. Effectively any parameter we might model with a linear model scales approximately linearly (other than through multiple scattering effects with the atmosphere and potentially between cover types). Any parameter we see in a non-linear sense in a model (e.g. exp(-LAI)) is unlikely to scale linearly.

Tian et al. (2002)[4] examined a theoretical basis for studying and understanding scaling of LAI. They consider a ‘pure’ 1 km resolution dataset composed of 6 biomes (plus bare ground) and characterize the mixing of these biomes at different spatial resolutions through a ‘purity index’ (mean of the cover type percentages). They then group pixels into three classes: ‘pure’ (group 1) (purity index > 90%), ‘mixed’ (group 2) (50% to 90%), and ‘heterogeneous’ (group 3) (<50%). A measure of heterogeneity is considered as the purity index of the dominant land cover type. They then define an ‘error’ in LAI through a relative difference between the ‘real’ and ‘inferred’ LAI values. This is shown in figure 5, for different dominant land cover types. We can see that, particularly for needle-leaf forests, when the ‘purity index’ is low, we can incur error of up to 40% in LAI due simply to scaling. A dependency on spatial resolution is shown, but the dominant effect can be attributed to ‘pixel purity’.

1.2.6 Options for Numerical Inversion

In the general case, we may have a complex, non-quadratic error surface in many dimensions. We could potentially go through all combinations of parameters to find which set gives a minimum in the error value, but may take far too long. However, if the error surface is known to be well-behaved, we may chose to sample the complete parameter space, using a so-called Look-up Table (LUT) inversion. This has been seen to be an attractive option for the inversion of optical models in recent years.

Often, we tend to use iterative numerical techniques to attempt to obtain the minimum of our error function. There are a variety of methods we can use, e.g. Quasi-Newton methods (an algorithm found in NAG which is used by a number of authors), or the Powell algorithm (an algorithm found in Numerical recipes in C, also used by many authors). Although such techniques are generally quite complex in nature, they are a valuable tool for the relatively rapid inversion of non-analytical functions. Basic aspects of these are discussed in greater detail below. Our personal experience is that the Powell algorithm is relatively robust, and tends to perform better than the NAG routine, but this may be case specific. Some improvements may be possible if we can find analytical solution to 1st and 2nd order partial derivatives, but most authors seem not to bother with this; in any case we often find that other difficulties also arise, especially in highly non-linear systems of equations.

Other options for numerical inversion, which are reviewed below, include: Knowledge-based Systems (KBS); Artificial Neural Networks (ANNs); and Genetic Algorithms (GAs). Selecting the right inversion algorithm for a particular problem is of critical importance in parameter estimation – one must weigh up the advantages and disadvantages of each method for a particular task.

1.3 How can we assess the inversion performance?

Generally, in order to assess the performance of a numerical inversion algorithm, we look at the 'goodness-of-fit' measure associated with the inversion )e.g. root-mean-squared error (RMSE)). However, we should understand that this is a global average of the error function, and doesn't tell us anything about how well the model fits the data at any particular sample location. This may well be critical if the model is particularly sensitive to errors at certain sample locations (e.g. parameters describing the hot-spot will tend only to be relevant in the region of the hot-spot). Obtaining an understanding of how critical particular angular samples are to the inversion is obviously very important if we wish to derive accurate estimates of the surface parameterization from RS data.

A relevant general comment at this point is that the more parameters we define as unknowns in our model, the more difficult it will be to derive accurate estimates of each parameter. This is because there will always be some degree of coupling between model parameters, such that errors in the estimation of any particular parameter can often be compensated for by errors introduced in other model parameters. Thus, there will always be a trade-off between the accuracy with which we can hope to model canopy BRDF with a limited number of parameters, and the maximum number of parameters we can expect to derive from our RS observations.

There are several other ways in which we can attempt to assess the performance of a model and the quality of parameters we invert from it. One way is to look at the implications of using a particular angular sampling set in our inversion, for instance, looking at how the model parameters vary as we have a greater or lesser number of samples at different solar zenith angles, as well as particular viewing angles. For any given sensor, at a given time of year, we will obtain a limited number of angular viewing samples. For many sensors, the range of sun angles of these observations will also be limited. We should therefore attempt to understand how well the different sampling capabilities of various sensors affect the parameters we obtain. It is apparent from many inversion studies, for example, that inverting the directional reflectance data for different sun angles can lead to significant variations in the derived model parameters.

The 'quality' of any inversion depends of course on the purpose of the inversion. If it is to derive biometric parameters from RS data, then we need to consider how accurately such parameters are derived. If however, the purpose is to use BRDF models to obtain estimates of the surface albedo, then we might indeed use a different set of criteria in assessing model performance. Goel (1988) discusses why it is important to understand the sensitivity of a model. This has many implications for model inversion, as, for instance, if our model was particularly sensitive to errors in knowledge of the viewing or sun angles, then any errors introduced in assigning these 'fixed' parameters would lead to significant errors in model inversion. We must also consider how sensitive a model is to any (residual) 'noise' in the system, i.e. any noise due to quantization effects or errors in atmospheric correction, any errors introduced by not taking the sensor MTF into account, etc.

1.4 How can we test the quality of the derived parameters?

There are several ways we might consider to assess how well the model is performing in inverting derived parameters.

One method is to use field radiometric data, where we can attempt to measure biophysical parameters of the canopy directly (e.g. LAI, leaf angle distribution) and, at the same time, obtain measurements of canopy reflectance. We can then invert the canopy reflectance measurements against our model and compare the resulting parameters with those we measured. There are several problems with such techniques, such as the accuracy with which we can hope to measure the average leaf angle distribution of a complex canopy, or determine the 'average leaf single-scattering albedo', but such methods can at least give an indication of how well a model performs under realistic conditions.

One can also try such methods of comparison with airborne or spaceborne remote sensing instruments, although the problems then tend to be compounded by increased uncertainties in atmospheric correction, pointing accuracy, etc., as well as being limited by the sampling capabilities of the sensor, as well as the smaller scale over which one has to obtain ground-based biometric measurements.

One can also test the performance of any particular model against other models, though one is obviously limited by the assumptions made in the original model formulation. It is perhaps more useful in this regard to make use of numerical solutions to the radiative transfer equation, and to use complex 3-D models of the vegetation canopy, because of the flexibility offered by such methods. When using numerical solutions to the radiative transfer equation, one can also attempt to test how well the particular assumptions made in any analytical approximation perform at a number of different levels, e.g. by simulating a canopy with constant LAI, known leaf angle distribution, and known single-scattering albedo (i.e. assume a Lambertian reflectance model). This can then provide us with information on how well the structural and radiometric assumptions made in the analytical formulation perform.