For detecting anomalies or interventions in the field of forest monitoring we propose an approach based on the spatial and temporal forecast of satellite time series data. For each pixel of the satellite image three different types of forecasts are provided, namely spatial, temporal and combined spatio-temporal forecast. Spatial forecast means that a clustering algorithm is used to group the time series data based on the features normalised difference vegetation index (NDVI) and the short-wave infrared band (SWIR). For estimation of the typical temporal trajectory of the NDVI and SWIR during the vegetation period of each spatial cluster, we apply several methods of functional data analysis including functional principal component analysis, and a novel form of random regression forests with online learning (streaming) capability. The temporal forecast is carried out by means of functional time series analysis and an autoregressive integrated moving average model. The combination of the temporal forecasts, which is based on the past of the considered pixel, and spatial forecasts, which is based on highly correlated pixels within one cluster and their past, is performed by functional data analysis, and a variant of random regression forests adapted to online learning capabilities. For evaluation of the methods, the approaches are applied to a study area in Germany for monitoring forest damages caused by wind-storm, and to a study area in Spain for monitoring forest fires.

Forests and other wooded land which cover over 40% of the European Union’s land area, are multifunctional and are serving economic, social and environmental purposes. The urgent need for monitoring forest health and vitality, including forest disturbances (i.e., forest fires, storm damage, drought stress, insect and disease outbreaks) is emphasised e.g. by the Green paper (European Commission, 2010) and in the new EU Forest Strategy (European Commission, 2013).

Currently, operational monitoring of forests is mainly applied on the basis of ground sampling data collected in national forest inventories. Typically, for several thousand randomly selected plots with a size below one hectare each, detailed forest parameters are measured in the field. The field measurements are time and cost intensive and typically the measurements, data processing, statistical analysis and reporting requires several years until the final results are provided. Whereas such inventories can provide accurate estimation of forest parameters at the regional level, they do not provide spatially explicit information. Further, in case of forest disturbances such as e.g. caused by storm damage or forest fires, up-to-date information on the damages is required. Remote sensing based monitoring methods can complement traditional field based national forest inventories for both aspects. Especially with the upcoming Sentinel2 satellites, remote sensing based forest monitoring capabilities will be improved significantly, as the Sentinel2 satellites will provide wall-to-wall data on forests with high spatial resolution of 10 m, high temporal resolution of 5 days and high radiometric resolution by providing measurements covering the optical, near infrared and short-wave infrared spectral regions.

Current remote sensing methods for monitoring forests are often based on the assessment of one to several image acquisitions. However, to fully utilise the comprehensive information content of up-coming Sentinel2 satellites, methods which exploit the information provided by dense time series are required. Dense time series are currently acquired by medium resolution sensors such as moderate resolution imaging spectrodiometer (MODIS) or Proba-V vegetation instrument, but at relatively low spatial resolution. For the current work, we used MODIS as pre-processed data is easily accessible. However, the tested methods can also be used with higher spatial resolution time series imagery, e.g. from Sentinel2 for future forest monitoring applications.

To allow qualitative and quantitative evaluations, the forecasting methods were applied to a test site in Germany for monitoring forest damages caused by wind storm, and to a test site in Spain for monitoring forest fires. The achieved accuracies were evaluated for both test sites for all methods on the basis of a sampling approach.

As basis for the methodological tests, we used time series of the MODIS product MCD43A4, Version 5, which provides 500 m reflectance data adjusted using a bidirectional reflectance distribution function. To reduce noise and to fill data gaps in the original MCD43A4 MODIS data (e.g. cloud covered areas), we applied Savitzky-Golay filtering (Savitzky and Golay, 1964) over n=15 MODIS MCD43A4 8 day mosaics. This leads to a significant noise reduction and provides as result a spatial and temporal consistent wall-to-wall cover also in areas with frequent cloud-cover such as, e.g. in the Alps. This filtering leads to a smoothed signal which is taken into account in the change detection process as described in the results section below. Of the time series data we used band 6, which covers the shortwave infrared spectral region (1628-1652 µm), and the normalised difference vegetation index (NDVI) (Kriegler et al., 1969) with:

where NIR is the reflectance in the near infrared spectral region and RED is the reflectance in the red spectral region. We selected the NDVI as index because of the long term proven robustness in various applications and ecosystem regions. Further, this index only requires data in the red and near infrared spectral regions, which is provided by a large number of operational satellites. We used the data over forested areas over the vegetation period from end of May until end of August in 8-day steps with temporal coverage from 2000 to 2013.

Two study areas were selected for testing and evaluation of the methods: Study area 1 is located in the south of Germany where storm damages occurred in the monitoring period. The MODIS time series covers this area as a raster with 541 x 541 pixels with a spatial resolution of 500 x 500 m on the ground. Study area 2 is located in northern Spain where forest fires frequently occur. Also this study area is covered with MODIS time series raster data of 541 x 541 pixels with a spatial resolution of 500 x 500 m.

Actually, functional data analysis (FDA) is an infinite dimensional (Hall and Hosseini-Nasab, 2006) analysis method. Main characteristics of functional data are repeated measurements, high frequency and multi dimensions, and that it is taken as function of an independent variable. An advantage of FDA is the usage of additional information of functions, e.g. the slope, curvature and other characteristics of curves.

If the measurements are dependent on time, we call them functional time series. With the technique of functional time series analysis (FTSA) it is possible to capture the underlying dynamic of seasonality in the data (Shang, 2013). Let Xt (x) be measurements, e.g. the vegetation index in year t with discrete observation number x. Let ft (x) denote functional data, σt (x) describes the quantity of noise varying with x, and ωt,i is a standard normal variable with zero mean and unit variance σω2. Then, the smooth function ft (x) can be extracted from the following equation:

where time t=1,…, n with n as the number of years considered, and i=1,…, p with p as the number of observations during one year. In Germany and Spain, the number of observations p equals 33. The smooth function ft (x) can be calculated using functional principal component analysis (FPCA), which reduces the infinite dimension of functional data to a finite level, and points to the most significant components of the data (Hall and Hosseini-Nasab, 2006). In the literature, there are numerous examples for FPCA with different application fields: sea surface temperature forecasting (Shang and Hyndman, 2011), call center arrivals forecasting (Shen and Huang, 2008), age-specific breast cancer mortality forecasting (Erbas et al., 2007) and fertility rates forecasting (Hyndman and Shang, 2009).

The functional principal component decomposition can be described through the following equation:

where

is the estimated mean function,

is the kth estimated orthonormal eigenfunction of the empirical covariance operator:

which maximises the variance,

is the kth principal component score for year t. This score is given by the projection of

in the direction of the kth eigenfunction

, that is

. The residual is denoted by

as iid random functions with zero mean, and the optimal number of principal components is denoted by K<n (Hall and Hosseini-Nasab, 2006; Shang, 2013). The estimated mean function can also be weighted unequally with geometrically decreasing weights wt=λ(1-λ)n-t with 0<λ<1, and can be written in the form

. Since we assume temporal autocorrelation, a reason for decreasing weights is the probability that more recent data affect the results more than data in the distant past (Hyndman and Shang, 2009).

The main aim of our calculations is to receive h-step-ahead forecasts following the equation:

by conditioning on the smooth functions f(x)=|f1 (x),…,fn(x)|T and the functional principal components Ф(x)=|Ф1(x),…,ФK(x)|T. The h-step-ahead forecasts of the principal component scores βn+h,k are denoted by

and are obtained by using an univariate time series model, e.g. an exponential smoothing (ETS) model. The purpose of the ETS technique is the estimation of a trend and seasonal component with additive and/or multiplicative terms (or none of them).

The error, trend and seasonal component follow the equations (Hyndman et al., 2002):

where m denotes the number of seasons each year, α, β, δ, and γ are constants. The values Lt, Kt, St, and Ct differ depending on the trend and seasonal component terms (Hyndman et al., 2002; Hyndman and Khandakar, 2008). From the h-step-forecast

and the variance of the prediction error

assuming Normal distribution, the lower and upper limits of the (1- α)-prediction interval can be calculated by

with qα/2 as the

th quantile of the Standard Normal distribution (Hyndman et al., 2008).

A classical and widely used approach to model time series is the autoregressive integrated moving average (ARIMA) model. It is a flexible model, which combines the autoregressive (AR) model, the moving average (MA) model and the autoregressive moving average (ARMA) model. These three main model types can be used if the time series data are stationary. But in practice time series often have a trend over time or a non-constant mean and hence they are non-stationary. Such non-stationary time series can be handled with an integrated ARMA(p,q) model, the ARIMA(p,d,q) model. In the ARIMA approach the non-stationary time series is differenced d times to obtain a stationary process. Generally an ARIMA model (Shumway, 2006) can be written as:

where Xt describes the time series, ωt is a Gaussian white noise series with mean zero and variance

and B is the backshift operator with Bk Xt = Xt–k. ф(B) and θ(B) are the autoregressive and moving average operators, respectively, with:

The ARIMA(p,d,q) model can be extended to an multiplicative seasonal autoregressive integrated moving average model to account for seasonal effects – so called SARIMA model, where autoregressive and moving average terms for seasonal effects are included in the same manner. The structure of this model – denoted as ARIMA(p,d,q)(P,D,Q) – is:

with autoregressive and moving average operators:

The model selection is done via Akaike information criterion (AIC).

Let the model forecast h steps ahead be

which is expressed through the conditional expectation of xn+h given all past information up to n. The variance of the prediction error is

. With the calculation of a prediction interval the precision of the forecast can be assessed. The lower and upper limits of the (1- α)-prediction interval are given by:

where qα/2 denotes the

quantile of the standard Normal distribution (Shumway, 2006).

The ARIMA model approach based on remote sensing data is applied for various questions. For example, Jiang et al. (2010) use a SARIMA model approach to model the leaf area index (LAI) time series of the years 2000 to 2006 and predict LAI values for the year 2007. The model and forecast results of the SARIMA approach are compared with two other methods with respect to different land cover types. Whereas the model performance of these three model approaches varies with land cover type, the SARIMA model and forecast perform very well overall. Han et al. (2010) use ARIMA models to predict drought in the Guanzhong Plain in China based on the vegetation temperature condition index (VTCI).

Random forest classifiers (RFC) (Breiman, 2001) are well known and studied classification methods. They show better or at least comparable performance in comparison to other state of the art classification algorithms such as support vector machines (SVM) (Vapnik, 2000) or boosting technologies (Freund and Schapire, 1997). They have been successfully applied to a number of applications. Moreover, they have also recently been adapted to density estimation, manifold learning, semi-supervised learning and regression tasks in a very successful manner. They have also shown their practical feasibility in famous commercial products such as e.g. Microsoft’s KINECT® (Microsoft Corp., Redmond, WA, USA). Criminisi et al. (2011) proposed a unified approach of random decision forests, which has been applied to a number of machine learning, computer vision and medical image analysis.

Analysing the change of a dependent variable (e.g. a certain forest monitoring parameter) with respect to a set of (multi-dimensional) independent variables (e.g. observations from satellite data) is a typical regression task. Thus it is natural to take into account regression variants of randomised decision trees as an alternative model generation technique. Since satellite data typically arrive sequentially and the processing time of model prediction becomes an important issue, multiple standard offline calculation of the regression is not feasible. Thus streaming or online learning capabilities become an essential requirement for that type of practical application.

We propose a novel regression method termed online random regression forests (ORRF) to build a proper, adaptive model for satellite data. The method follows to some extend the generic approach of Criminisi’s framework (2011) for transferring standard classification forests to regression problems. But coevally some of the concepts introduced by Saffari et al. (2009) have to be used to change structure, layout and learning rules of the forest in order to enable the algorithm for having proper online learning capabilities.

Random (classification) forests are an ensemble combination of several binary decision trees where each binary decision tree is treated independently. The final output decision of the forest is obtained by e.g. a simple majority voting of all the individual leaf node predictions or other combinations with respect to probabilities.

Each tree itself consists of a singular root node, which subsequently splits up into two child nodes in a hierarchical manner. A simple test function on a training sample is applied to decide about the path the sample moves down the tree. During the offline training, the best split for each individual tree node is calculated by globally optimising several random test functions with respect to the overall information gain obtained by each split.

Taking into account a large number of samples continuously arriving over time, the main disadvantages of such offline optimised classification forests are the increasing calculation time and the necessity for storing all samples. If boundary conditions change or if additional retraining is done the required execution time typically exceeds the processing capabilities of a system.

Online adaptation of classification forests has been introduced by Saffari et al. (2009). The authors combined the ideas from online bagging and extremely randomised forests. They proposed a novel procedure (algorithm) for growing a decision tree in an online fashion for visual tracking and interactive real-time segmentation tasks.

For our approach we use the basic idea of online bagging in a similar way. In particular we use the basic principle of modelling the sequential arrival of data by Poisson distribution sampling according to Oza and Russell (2001). This allows continuous growing and updating the tree structure.

Although the main difference between classification and regression trees is the change of the output labels from discrete class labelling to continuous prediction values, the scheme of Criminisi’s framework (2011) cannot directly be applied due to numerous reasons. In the following we mention the most important ones for our problem and describe the solutions to deal with.

First, the strict online learning constraints (strict online learning means, that it is not allowed/possible to store any data samples for later usage. Data has to be processed on the fly with no buffering as e.g. used in incremental learning) require for an incremental update of the individual probabilities and statistics for each node of the forest. This is no problem for the classification case, because class-statistics are discrete histograms and their update is trivial.

As regression trees use continuous densities it is necessary to store them and to enable incremental update without having access to previous data points. Although it would be possible to sample the density function and store as a regular raster, this approximation of the density function turns out to be not feasible because of a high memory amount required. We approximate the density probability function by a simple Gaussian distribution with mean µ and variance σ². The two parameters can be incrementally updated via the following formulas:

where wn is the weight of the actual sample xn in iteration n of the node, µn denotes the mean of sample xn and σn is the respective standard deviation.

A second important issue of changing online classification to regression behavior is the proper choice of a prediction model for each leaf node. While each node in higher levels of the tree directs an incoming sample down the tree, the leaf node has to deliver a proper prediction of the continuous output as well as the respective confidence. The forms of prediction functions of each leaf node can vary e.g. polynomial functions as proposed in Criminisi et al. (2011). It makes sense to adapt that function to the expected target function. In general a simple constant model related to μn is sufficient, if the need for higher forest complexity is less important. The regression forest’s overall posterior is then calculated by the median of all individual tree posteriors.

Another important aspect is the proper choice of an objective function for maximising information when splitting intermediate nodes during online training of the forest/tree. In the offline case, the optimal split can be determined during the global training step.

In contrast to the offline case an early split decision has to be made during the samples arriving. The criteria for performing a split are similar to the ones in Saffari’s approach (2009). In particular a node only splits if i) a minimum number of samples has already passed the node (ensures statistical significance), ii) the depth of the tree has not exceeded the predetermined maximum model complexity (ensures a final size model), and iii) the minimal information gain required by a split is reached (avoids early growing).

The objective function for the information gain is the sum of all n label’s differences between actual (y) and the node’s labels (μ), for all parent (P) and left/right (L, R) child (C) nodes.

Another challenge introduced by stream processing tasks is the unknown range of feature dimensions. Extremely random forest approaches select random dimensions and random thresholds of the samples for decision candidates (weak learners). In the offline case, feature ranges can be trivially calculated from the whole dataset, but in the online case this is usually not possible. It is reasonable to carefully select and adapt the feature dimension’s boundaries, to avoid useless calculations and a waste of memory. We allow changing the feature dimension boundaries during the growth of the tree, starting with very conservative feature ranges and extending them dynamically to the highest/lowest values seen so far.

One problem introduced by this strategy is the fact, that it changes the tree configuration over time and the actual splitting criterions might no longer be optimal. Thus it is necessary to remove certain trees from time to time. Additionally we allow a novel tree to learn focusing optimally on novel examples with updated tree configuration. Therefore we introduce a jig-saw criterion, calculated on the out-of-bag error value calculation (Oza and Russell, 2001). Thus we also are able to react on slightly changing boundary constraints of the problem (e.g. changing geological conditions when drilling down to the earth). It is also possible to derive some measure for the overall prediction quality of the forest by averaging or taking the mean of the individual trees’ out-of-bag errors already calculated for the jig-saw criterion introduced above.

The detection of unusualness of the vegetation index and short-wave infrared band with above described methods follows a common scheme. First of all the NDVI and SWIR are modelled from the satellite data and the expected time series for a certain year is forecasted using the FTSA, ARIMA or ORRF method. Therefore the 80, 95 and 99% prediction intervals, as already mentioned in the respective method’s descriptions above, are calculated for the FTSA and ARIMA method. Since ORRF is a non-parametric method, it turns out, that a much simpler estimation of the prediction intervals as the one proposed, leads to more stable and better results. In particular, we validate the ORRF’s predictions on all the training/update data of the preceding year by comparing all real training values to the ORRF’s model prediction for all pixels and time points. The upper and lower boundaries for including 80, 95 and 99% of the training values (with respect to the ORRF prediction value at each observation time point in the year) can then be treated as estimation for the prediction boundaries in the actual year of observation. More formal, let

be the h-step forecast and

the estimation error, then the empirical α-quantiles of the prediction error qr,α are calculated. Hence, the lower and upper boundaries of the (1-α)-prediction interval follow the equations

approximately.

Subsequently, for all methods unusualness is detected, if the original time series is outside the prediction boundaries. The forecast models are done solely spatial (FTSA and ORRF), solely temporal (FTSA and ARIMA) and combined spatio-temporal (FTSA and ORRF).

The results of above methods provide for each pixel for each observation the probability that the reflectance in the SWIR spectral region or NDVI has changed significantly over time. To reduce noise which is inherent in the time series measurements (such as e.g. caused by locally varying atmospheric conditions during image acquisition, sensor noise, etc.) change detection is not based on single observations but on the aggregation of three observations, as the forest changes to be observed are in general long term processes (e.g. crown cover recovery needs several years after a storm damage). A pixel was therefore classified as changed only in case of at least three consecutive results that the respective pixel is an outlier at the 95% prediction interval.

Figure 1 shows Spain study site where forest fire happened in 2012 (red colored pixels in the black cycle). Pixels with green color mean that no unusualness was detected with the forecasting methods whereby yellow, orange and red mean that the original vegetation index time series was significantly different from unusualness. From that region we consider one certain pixel corresponding to the region indicated by the yellow bounding box in Figure 2.

For this certain pixel, Figures 3-5 show the forecasts of each of the three model types (spatial, temporal and spatio-temporal) using the corresponding methods (FTSA, ARIMA, ORRF) for 2012. Additionally, the prediction intervals are plotted as dotted lines. Unusualness happened if the original time series is above the upper or below the lower prediction boundaries. As the time series data are pre-processed for noise reduction and gap-filling over time (e.g. gap-filling in case of cloud cover at the satellite image acquisition) by Savitzky-Golay filtering (Savitzky and Golay, 1964), the change of the signal is not abrupt. We can see that phenomenon in all three plots starting at about observation number 25 when forest fire occurred. The difference between FTSA and ORRF in the case of the spatio-temporal forecast (Figure 5) is caused by the different adaptation sensitivities of the forecast methods. Figure 6 shows the original time series and the FTSA prediction with the 80, 95 and 99% prediction intervals. The temporal forecast and the three prediction intervals based on ARIMA are shown in Figure 7. The spatio-temporal forecast and its prediction intervals based on ORRF can be seen in Figure 8.

A qualitative evaluation of the monitoring results was performed by comparing the detected changes with time series of higher resolution Landsat thematic mapper (TM) satellite imagery with spatial resolution of 30 m by 30 m per pixel. These comparisons have shown that the derived change detection results correspond well with the changes visible in the TM imagery.

For quantitative evaluation of the different monitoring methods in addition a sampling based approach was applied. As only small parts of the study areas are affected by damages, stratification was performed as a first step. As basis for the stratification, a forest loss map generated by Hansen et al. (2013) who mapped forest loss at the Landsat pixel scale was used. These maps were aggregated according to the MODIS spatial resolution (e.g. from 30m by 30m pixels to the MODIS resolution). As the qualitative evaluation has shown that the change detection methods are not only capable of indicating changes at the pixel level but also at the sub-pixel level, stratification was applied for the strata no forest loss, forest loss below 20% and forest loss above 20%. Within the stratum no forest loss and the stratum forest loss above 20% 25 grid cells were selected in both strata randomly. As in the German study site only small areas are affected by storm-damages, the stratum forest loss above 20% covered only 23 MODIS pixels. For this stratum, no random sampling was therefore performed but all 23 grid cells were selected for the accuracy assessment. For each of the selected MODIS grid cells, the percentage of forest loss mapped by Hansen et al. (2013) was calculated by aggregation. Then, the grid cells were visually checked on the basis of Landsat TM time-series data and the MODIS time series and the Hansen et al. (2013) estimate was revised in cases of clear deviations. Such deviations occurred e.g. in case of misclassification in the Hansen et al. (2013) map. For both study areas then contingency tables were derived for all selected change detection methods (Table 1). As a measure of the achieved accuracy, the percentage of grid cells, which were classified correct was calculated and is given in percent.

A main factor influencing the achieved accuracy is the nomenclature definition of the category forest loss. In the current study we selected a very ambitious threshold for the definition of this category with forest loss indicated already in case that only 20% of the area covered by one pixel are affected e.g. by a storm damage or forest fire. The current methods are therefore performing change detection at the sub-pixel level. In the following, the results are shortly discussed for each of the monitoring methods.

Exceptional high accuracy was achieved for both study areas with the combined spatio-temporal as well as with the temporal approach. No significant errors could be identified. For the spatial approach, the achieved accuracies are in the range of current operational monitoring methods with overall accuracies between 72 and 92 percent. Processing time is a drawback for functional time series analysis methods, which requires processing time of about 8 hours per study area on a usual PC (i7 processor). Software is currently implemented in the R packages ftsa (Shang, 2013) and forecast (Hyndman et al., 2008), where each pixel is treated separately. For large area applicability a pre-requisite would be a significant optimisation of the implementation, e.g. via graphics processing unit (GPU) processing.

The qualitative evaluation by comparison of the monitoring result with the Landsat time series showed that the ARIMA method is less sensitive to sub-pixel changes compared to the functional time series analysis approach. However, the quantitative accuracies with 86% in the Spain study area and 79% in the German study area are high as we perform change detection at the sub-pixel level. In case that the nomenclature specification of the forest loss category is defined with a minimum change of 50% instead of 20%, overall accuracies for both study sites are above 96%. Main drawback for operational application is the processing time of about 3 hours per study area on a usual PC (i7 processor). The R-package forecast is used (Hyndman et al., 2008). As with the functional time series analysis method, significant optimisation with (GPU) processing would be required for operational application.

Accuracy of the spatial approach is 70% for the German and 92% for the Spain study site. Qualitative evaluation has shown a local bias of the change detection result (overestimation of changes in the southern area of the Spain study site). For the spatio-temporal approach, high accuracies were achieved with overall accuracies of 76% in the Spain and 94% in the German study site. As with the ARIMA method, accuracy increases significantly in case of changed nomenclature specification. In case that the nomenclature specification of the forest loss category is defined with a minimum change of 50% instead of 20%, overall accuracies for both study sites are above 94%. A main advantage of the online regression forest method is that the required processing time is low, in the range of several seconds per study area on a usual PC. This method allows therefore near real time processing.

Exceptional high accuracies were achieved with the temporal and spatio-temporal functional time series analysis methods, which however need significant optimisation of the software implementation to allow large area applicability (GPU processing). It is expected that near real time monitoring based on upcoming Sentinel2 satellite data can be achieved in this case at the local to regional level. Implementation can be recommended for such applications but also for applications based on medium resolution satellite data, e.g. MODIS or upcoming Sentinel3 time series data. For near real time, continuous forest monitoring at the national to continental level based on upcoming Sentinel2 satellite data the implementation of spatio-temporal online random forest method is recommended, however because of sensitivity issues, monitoring should focus on severe changes only (e.g. monitoring at the pixel level with change fraction above 50%).

The results demonstrate the applicability of dense time series of remote sensing data for large area forest monitoring with spatial-explicit wall-to-wall coverage. Compared to more simple methods, which are often based on bi-temporal comparison of satellite imagery, the methods applied in this publication are based on wall-to-wall, fully consistent time series and allow fully automatic processing. This automation is a specific requirement for future forest monitoring applications based on satellite imagery with high spatial and high temporal resolution such as e.g. Sentinel2 time series.