The polynomial least squares operation (PoLeS): A method for reducing noise in NDVI time series data
Profile extraction: The polynomial least squares operation
The premise of data smoothing is that one is measuring a variable that is both slowly varying and also
corrupted by random noise. Rather than having their properties defined in the Fourier domain, and then
translated to the time domain, PoLeS derive directly from a particular formulation of data smoothing
problem in the time domain. PoLeS filters were initially used to render visible the relative widths and
heights of spectral lines in noisy chemical spectrometric data.
In principle, a digital filter is applied to a series of equally spaced data values fi= f(ti), where ti=to +i
for some constant sample spacing and i = …-2, -1, 0, 1, 2, …A digital filter replaces each data
value fi by a linear combinationgi of itself and some number of nearby neighbors,
(1)
Here nL is the number of points used “to the left” or “backward in time” of a data point i, i.e. earlier than it,
while nR is the number used “to the right” or “forward in time”, i.e. later or after the data point.
In order to understand the PoLeS, consider the simplest possible averaging procedure: For some fixed
nL = nR, compute each gI as the average of data points from fi – nL to fi+ nR. This is sometimes called the
moving window averaging and corresponds to equation 1 with constant cn = 1 / (nL + nR + 1). If the
underlying function is constant , or is changing linearly with time (increasing or decreasing), then no bias
is introduced into the result. Higher points at the one end of the averaging interval are on the average
balanced by lower points at the other end. A bias is introduced , however, if the underlying function has a
nonzero second derivative. At a local maximum, for example, moving window averaging always reduces
the function value. In the spectrometric application, a narrow spectral line has its height reduced and its
width increased. Since these parameters are themselves of physical interest, the biased introduced is
distinctly undesirable.
Note, however, that the moving window averaging does preserve the area under the spectral line,
which is its
zeroth moment, and also (if the window is symmetric with nL = nR ) its mean position in time , which is
its first moment. What is violated is the first moment, equivalent to the line width.
The idea of PoLeS is to find the filter coefficients Cn that can preserve higher moments. Equivalently,
the idea is to approximate the underlying function within the moving window not by a constant (whose
estimate is the average), but by a polynomial of higher order, typically quadratic or quartic: For each
point fi, we least squares fit a polynomial to all nL + nR + 1 points in the moving window, and then set gi
to be the value of the polynomial at position i. The value of the polynomial is not used at any other point.
Moving to the next point, fi + 1, the least squares process is applied using a shifted window.
The process of least squares fitting involves a linear matrix inversion, the coefficients of a fitted
polynomial are themselves linear in the values of the data, which would mean that coefficients can be
derived which could be applied to different sets of data as some kind of filter coefficients.
To derive such coefficients, consider how go might be obtained: By fitting a polynomial of degree M in i,
namely ao + a1i + . . . + aMi to the values f-nL . . ., fnR. Then go will be the value of that polynomial at i =
0, namelya0. The design matrix for this problem is
and the normal equations for the vector of aj’s in terms of the vector of fi’s is in the matrix notation.
we also have the specific forms
k j fk.Since the coefficients cn is the component a0 when f is replaced by the unit vector en, - nL < n < nR, we
have
NOAA-AVHRR Pathfinder NDVI 10-day time series data for the year 1984 were taken and were
geometrically corrected. 3X3 Sample data different land cover types namely, for Central Africa (desert),
Eastern Canada (deciduous), Philippines (two-time cropping), and Central Brazil (tropical) were collected
and annual NDV profiles were constructed. Agbu et al (1994) describes the Pathfinder format.
The typically noisy profile can be seen in the original data shown in the figure1. The PoLeS algorithm
assumes (1) that cloud and cloud shadow and poor atmospheric conditions will usually depress NDVI
values (Arino et al. 1991, Simpson et al 2000); (2) that there exists vestiges (relics) of various
systematic and data transmissions errors which tend to increase the NDVI values. Lovell et al (2001)
had explicitly indicated that the 10-day compositing period (based on the MVC method) did not remove
all cloud and spurious high values (e.g. data transmission errors). Other systematic errors include
anistropic (describes a medium which has physical properties different in different directions) sources
such as the Bi-directional Reflectance Distribution Function (BRDF) (Macmanus 2001). Likewise,
scattering within the atmosphere is said to be anistrophic (CCRS, 2002). Similar effects of increased
value NDVI due to changes in satellite view angle have been reported (METECH, 2001).
The surface effects such as BRDF and Bi-directional Emission Distribution Function (BEDF) can all be
looked at as varying scene brightness that does not indicate changes in the surface type and condition
(phenology) but rather the varying sun condition and view geometry of the data. AVHRR and airborne
data both share the problems of widely varying scene brightness. Canadian researches have indicated
about 30% of the variation in AVHRR based NDVI can be due to sun and view angle effects (Jupp,
1997). Similar anomalies have been found in analyses by Epiphanio et al (1995) and Eastman and Fulk
(1993). While the effects of changing view and illumination are sometimes mentioned in studies using
AVHRR, the reference is usually made in passing and other disturbance factors are sought to try to
explain the anomalies apparent .
In view of the fact that most of the previously proposed NDVI profile extraction methodologies have the
intrinsic bias of considering only low frequency NDVI values as noise data and, with the assumption
that the surface characteristics are non-isotropic. Previous NDVI profile extraction algorithms also
assume that the peaks interpolated from and which conform to the conditions set by the various
algorithm (i.e. BISE and MVI) correspond to the absolute values of NDVI curve at those particular
points. The PoLeS method, on the other hand, assumes that both high frequency as well as low
frequency noise exist in the PAL NDVI data, and should therefore be considered in the interpolation of
the NDVI profile.
