Training Area Concept in a Two-Phase Biomass Inventory Using Airborne Laser Scanning and RapidEye Satellite Data

This study evaluated the accuracy of boreal forest above-ground biomass (AGB) and volume estimates obtained using airborne laser scanning (ALS) and RapidEye data in a two-phase sampling method. Linear regression-based estimation was employed using an independent validation dataset and the performance was evaluated by assessing the bias and the root mean square error (RMSE). In the phase I, ALS data from 50 field plots were used to predict AGB and volume for the 200 surrogate plots. In the phase II, the ALS-simulated surrogate plots were used as a ground-truth to estimate AGB and volume from the RapidEye data for the study area. The resulting RapidEye models were validated against a separate set of 28 plots. The RapidEye models showed a promising accuracy with a relative RMSE of 19%–20% for both volume and AGB. The evaluated concept of biomass inventory would be useful to support future forest monitoring and decision making for sustainable use of forest resources.


Introduction
Forests play an important role in global carbon cycling, since the world's forests sequestrate and conserve more carbon than all other terrestrial ecosystems, and account for 90% of the annual carbon flux between the atmosphere and the Earth's land surface [1].Accurate estimation of forest biomass is required for countries ratifying the Kyoto Protocol to the United Nations Framework Convention on Climate change [2].Assessment of biomass within all compartments of forest ecosystems can be achieved by choosing the appropriate scientific techniques.The combination of remote sensing with ground surveys potentially offers way to guarantee accurate monitoring, reporting and verification of terrestrial forest biomass [3][4][5].Southworth and Gibbes [6] reviewed the past, present and future directions of remote sensing applications in the forestry sector.They concluded that the fusion of different remote sensing techniques is promising for the monitoring and assessing of forest ecosystems.
The estimation of forest carbon is still relatively uncertain considering the errors in regional carbon stock estimates.Houghton [7] mentioned that 89% of the carbon losses are due to the loss of living biomass; therefore, attempts have been centered on estimating the above-ground biomass (AGB) of vegetation [8].Accurate large-scale estimation of AGB using active or passive remote sensing has been a significant challenge to forest scientists; however, it is crucial for future implementation of carbon credit verification in the Land-Use Change and Forestry (LUCF) sector [9].The purpose of verifying forest biomass inventory is to establish their reliability and to monitor the accuracy of the numbers reported by independent means.
Finland, the northernmost European country, has 73% of its land area forested [10].In Finnish boreal forests, coniferous trees represent over 80% of the growing stock.Globally, boreal forests are of particular interest because among all the biomes, they may undergo the greatest climatically induced changes during the 21st century [11].Boreal coniferous forests are the largest terrestrial biome on Earth and store 31 × 10 12 kg of carbon in the trees alone [12].Research indicates that the boreal forests are more significant as carbon sinks than has been thought previously [12,13].
A wide range of approaches have been proposed for quantifying forest biomass using active and passive remote sensing systems.Although many alternative remote sensing techniques have been suggested for the estimation of forest attributes such as AGB, the most promising one seems to be the airborne laser scanning (ALS) [14].ALS sensors provide accurate information on the 3D structure of the vegetation, as they directly measure distances to the targets below.ALS has potential to increase accuracy and reduce costs in large-scale forest inventories [13][14][15][16][17].The two established approaches for predicting forest biomass using ALS are the area-based method (ABM) and the individual tree delineation (ITD).Both approaches have been proven to provide competent results, but large-scale applications are usually conducted using the ABM [18].In this approach, field-measured data and ALS metrics are used to make statistical models that describe the relationship between ALS metrics (used as predictors) and field-measured forest properties (e.g., AGB).
The RapidEye is a multispectral satellite sensor launched on 2008 with an ability to acquire images with high spatial resolution (5 m) on five spectral bands.The RapidEye images have been used in several forestry operations, including cost-effective monitoring, harvesting and mapping [19,20].The RapidEye data might be easier to obtain in many countries than other remotely sensed data, such as ALS data.The drawback of all optical satellite images in the estimation of biomass is the saturation of predicted values at dense leaf canopies, which restricts the estimates to low biomass levels [21,22].
The most commonly used features derived from the optical data to predict forest attributes are the spectral and textural features.Spectral features describe the tonal variation in portions of the electromagnetic spectrum.Textural features contain information about the spatial distribution of tonal variations within an image.Texture has qualities such as periodicity and scale; it can therefore describe, for example, the direction, coarseness, and contrast of image components [23].The combination of spectral and textural features has provided better accuracy in the estimation of forest attributes than the use of any feature alone.Haralick et al. [24] presented the use of grey-level co-occurrence matrices in quantifying the texture.Their method has been widely used in remote sensing-based forestry applications [25,26].
Albeit the costs of ALS data have decreased considerably, it has relatively high cost compared to the optical sensors.Decreasing the number of field reference plots to minimize the costs of data collection is desirable, as the cost of field measurements was estimated to be around 100 € per sample plots (9 m radius) in Finland [27].Tegel [28] mentioned that if the number of field reference plots is not sufficient, the accuracy of large-sale inventory results can be improved by combining satellite imagery with ALS data and field measurements, instead of using just the field plots.
The objective of this study was to assess the AGB and volume using ALS data and RapidEye satellite data in a two-phase sampling procedure.ALS-estimated sample plots (called surrogate plots) were used as a simulated ground-truth instead of more expensive field sample plots.A small number of field sample plots were collected to create ALS-based models, which were applied to predict the AGB and volume for the surrogate plots.Finally, RapidEye image and the surrogate plots were used to generalize the predictions for the study area.
The accuracy of the predictions was evaluated by means of root mean square error (RMSE) and bias.To our knowledge, there are few studies focused on using ALS data as a simulated ground-truth for the satellite-based forest inventory (e.g., [22]).This study is aimed at obtaining an overall view on the accuracy of two-phase biomass inventory using ALS and RapidEye satellite data.

Study Area and Field Data
The study area is located in eastern Finland (62°31'N, 30°10'E) (Figure 1).The boreal forests of the area are managed for timber production and ecological sustainability.Scots pine (Pinus sylvestris L.) is the dominant tree species, representing 73% of the volume.Norway spruce (Picea abies L.) represents 16% of the volume, and rest come from deciduous trees such as downy birch (Betula pubescens Ehrh.) and silver birch (Betula pendula Roth.), which usually occur as minor species.
The field measurements were carried out in May to June 2010.Altogether, 78 field plots were placed subjectively into different stands in an attempt to record the species and size variation over the area.The field plots were placed into the area subjectively based on the development class and dominant tree species.The sizes of the field plots ranged from 20 × 20 to 30 × 30 m. Field sample plots were divided into training (n = 50) (phase I) and validation (n = 28) data sets; the training set comprised the plots whose size was 25 × 25 m (0.065 ha).The other plots with varying sizes were used for validation.It is worth noting that the larger sample plots maintain a greater amount of spatial overlap, minimize the edge effects, and increase the sample variances [29].A total of 200 surrogate plots (35 × 35 m) were placed to cover the study area using an ortho-rectified aerial photograph with the green, red and near-infrared portions of the spectrum.These plots were used as a simulated ground-truth and while training the RapidEye models in the phase II (described in Section 3.4).We used visual interpretation of aerial photographs in the placement of the surrogate plots.The idea was that the surrogate plots should cover all the variation (i.e., vegetation type, geographic location, species composition) in the study area.
All trees with either diameter at breast height (DHB) ≥ 4 cm or height ≥ 4 m were measured in the field.The volumes of the individual trees were calculated as a function of diameter at breast height (DBH) and tree height using species-specific models [30].The AGB of individuals trees was calculated by using the biomass Equations (1-3) developed by Repola [31,32], respectively, for Scots pine, Norway spruce and deciduous.Repola [31,32] where y ki is the total AGB (kg) of tree i in stand k, b 0 , b 1 and b 2 are the vectors of fixed effects parameters, d ski is (2 + 1.25 × DBH), h ki is the height (m) of tree, u k is the variance of random parameters, and e ki is the residual error.
Table 1 shows the multivariate model fixed parameters and residual errors values.Finally, total volume and AGB were calculated for each plot and stands per hectare.The characteristics of the training plots and validation stands are presented in Table 2.
Table 1.Estimates of multivariate model fixed parameters, and variances of random stand parameters (u k ) and residual errors (e ki ) [31,32].

Remote Sensing Data
The ALS data was collected on 18 July 2009 using an Optech ALTM Gemini laser scanning system.The nominal pulse density was about 0.65 per square meter.The test site was scanned from an altitude that is approximately 2000 m above ground level with a field view of 30 degrees and side overlap between transects of 20%.Pulse repetition frequency was set to 50 kHz.The swath width was approximately 1,050 m.
The RapidEye satellite images were collected on 19 May 2012 for the test area.The RapidEye image index numbers were 2012-05-19T102327_RE5_3C-N05_9429301_137127 and 2012-05-19T102330_RE5_3C-N05_9429336_137127.RapidEye imagery provides multispectral optical imagery of five bands (blue 440-510 nm, green 520-590 nm, red 630-685 nm, red-edge 690-730 nm, and near infrared 760-850 nm).A total of two RapidEye images were collected with spatial resolution of five meters to cover the study area.All the RapidEye images were radiometrically and geometrically corrected (overall standard error was 0.53 m) according to the standard of RapidEye [33] and aligned to a cartographic map projection.The RapidEye satellite orbit altitude was 630 km in sun-synchronous orbit with a swath width of 77 km.

Preprocessing and Estimation of ALS Predictors
First points were classified as ground and non-ground hits according to the approach described by Axelsson [34].The height (z) values of the laser echoes were changed to the altitude above the DTM (Digital Terrain Model) surface instead of the altitude from the base geoid of the projection.A raster DTM was then obtained by interpolation using Delaunay triangulation.The spatial resolution of the DTM was 0.5 m.Corresponding DTM elevations were subtracted from the ellipsoidal heights of the ALS points to scale their elevations to the above ground heights used for analysis.
The area-based method was used to model the relationships between field-measured variables (e.g., AGB) and canopy height/density metrics from the ALS data [35].The first echo data included first-of-many and single echoes whereas the last echo data included last-of-many and single echoes, while all intermediate echoes were ignored.A total of 46 ALS explanatory variables were extracted based on ALS variables defined by Junttila et al. [36] and Naesset [37] with slight modification and are described in Table 3. Ratio of last and single echoes with height lower than 1.5 m + i × 3 m for i = 0.7 and the total number of last and single echoes.

39...41 D fp10,30,50
Ratio of first and single echoes with intensity I ≤ 0.5+i for i = 10, 30, 50 and the total number of first and single echoes.

42...44 D lp10,30,50
Ratio of last and single echoes with intensity I ≤ 0.5+i for i = 10, 30, 50 and the total number of last and single echoes.

D flog
Logarithm of the ratio of the number of first and single echoes below 5 m (low vegetation) and the total number of first and single echoes.

H f3mean
Mean of the largest three heights within first and single echoes.

RapidEye Image Preprocessing
Before using the RapidEye images in the final calculation, the necessity of radiometric correction was examined according to the Ridge method [38].The main idea is the linear relationship of the digital numbers (DNs) over pseudo-invariant features (PIFs) across these two images [39,40].We made five feature space (for five RapidEye bands) images within these two images in the overlapping area.We used one band as the X axis and the same band in the other image as a Y axis to create the feature space image using the ERDAS IMAGINE software (version 10.1) and with the different colors representing histogram frequencies.We found that the slope and intercepts in the resulting feature space images (Figure 2a-e) were always one and zero, respectively.So, we determined that the radiometric correction was unnecessary.We also created the histograms and found a similar distribution in the overlapping area of these two images, which confirmed our previous results from the feature space images (e.g., blue band, Figure 3).

Estimation of RapidEye predictors
The spectral and textural features were calculated from the RapidEye images and used as predictors for modelling.The extraction was done based on the field plot size and RapidEye image pixel size so that the extracted image value could properly represent the forest attributes (e.g., AGB) at the plot level.The spectral predictors were derived from each band of RapidEye image by taking the mean DN values based on the plot size.Three spectral vegetation indexes were computed from the RapidEye images.Normalized Difference Vegetation Index (NDVI) (Equation ( 4)) is strongly related to the photo-synthetically active radiation intercepted by live vegetation, and is thus a robust method to identify live canopies in multispectral images [41].NDVI has been widely used for estimating forest biomass [17].NDVI is calculated from the visible (red) and near-infrared (NIR) bands of the RapidEye sensor (Equation ( 4)); NDVI thus varies between −1.0 and +1.0, i.e., green vegetation gives a positive value, snow and clouds tends to give negative values, whereas water and soil give values close to zero.
Also, the red-edge and green band ratio in the RapidEye images has a good response to forest biomass.Thus, the second NDVI [42] was calculated from these individual measurements (Equation ( 5)).In addition, the third NDVI [43] was calculated from the red-edge and red channel of the RapidEye image (Equation ( 6)) and is well-known as Normalized Difference Red Edge Index (NDRE).Given the sensitivity of the NDRE, this index is used for a variety of forest applications including biomass mapping and forest health.
We ran the Pearson correlation test to confirm which NDVI has the better correlation with AGB and volume.The Pearson correlation showed that the first NDVI and third NDVI had better correlations (r = −0.32 and −0.34, respectively) compared to the second NDVI (r = −0.12).A total of 14 textural features [24] were computed from the first and the third NDVI, separately.
The Haralick textural features were derived from the RapidEye images by considering the mean DNs based on the plot size.The co-occurrence matrix was used as an input to compute the Haralick's features.We used single spectral bands (mean), three vegetation indices (e.g., normalized difference vegetation index) and Haralick textural features from RapidEye satellite data.Finally, a total of 36 RapidEye explanatory predictors (Table 4) including both spectral and textural features were available for the prediction of the volume and AGB.

Two-Phase Sampling Method
Figure 4 shows the flowchart of the two-phase sampling method.In the first phase of this approach, a regression model is generated based on the relationship between ALS-metrics and field-measured sample plots.In the second phase, the forest characteristics that are estimated for the "surrogate plots" from ALS data are applied as simulated ground-truth to generate a regression model between forest parameters (e.g., AGB) and features derived from RapidEye satellite imagery.

Statistical Modeling
Ordinary least squares (OLS) is the most commonly employed method for estimating the unknown parameters of a linear regression model.The OLS can be written as follows in Equation ( 7); where = field values of volume/AGB, is a vector of regression coefficients; ˊ is a matrix of the explanatory variables from ALS data as well as from RapidEye data; is an unobserved error which account for the difference between the actually observed responses and the predicted values ˊ .We did not employ non-parametric method due to the insufficient number of training plots.
As we had a large number of predictors, we used the leaps package (regsubsets algorithm) in R (Version 2.15.2) environment to select the best combination of predictors for the models [44].It performs an exhaustive search for the best subsets of the explanatory variables in X for predicting Y (dependent variable) in linear regression, using an efficient branch-and-bound algorithm.The selection criteria in the algorithm were Bayesian Information Criterion (BIC) and adjusted R 2 (coefficient of determination) value.Since the algorithm returns the best model of each size (number of predictors), the results do not depend on a penalty model for model size.In addition to stepwise variable selection, we used the Pearson correlation techniques and the maximum R 2 improvement variable selection techniques to select ALS/RapidEye-derived variables to be included in the models.During the variable selection, no explanatory variable was left in the models with a partial F statistic with a significance level greater than 0.05.It is worth noting that regsubsets algorithms might have many independent variables that could introduce multicollinearity problems, i.e., one independent variable is a linear combination of the other independent variables.The problem is that, as the independent variables (e.g., ALS height percentiles) become more highly correlated to the dependent variables (e.g., AGB), it becomes more and more difficult to determine which independent variable is actually producing the effect on the dependent variables.In this circumstance, we used maximum R 2 improvement techniques to search for the "best" one-variable model, the "best" two-variable model, and so forth, although it is not guaranteed that in finds the model with the highest R 2 for each size.We repeated the selection 10 times by randomly selecting 50% of the observations for each repetition to produce stable models for our final regression models.Our final models contained the most frequently occurring ALS/ RapidEye-derived variables.Naesset et al. [45] and Latifi et al. [46] used similar predictor selection techniques for model building to predict forest biomass and volume.

Model Accuracy Assessment
The RMSE and relative RMSE (Equation ( 8)) have been widely used for validating the accuracy of forest parameters estimation.The accuracy of fitted models was measured by using leave-one-out validation (LOOV) techniques.Bias is the mean residual error between predicted and estimated model values (Equation ( 9)).Relative bias gives the percentage of mean residual error (Equation ( 9)).
where y i is the predicted value for sample plot i; is the observed value for sample plot i; is the average value of measured sample plots, and n is the total number of plots.The statistical significance of the bias was estimated by the t-test (Equation ( 10)).
where SD is the standard deviation of the residuals ( − ).The bias was considered to be significant if the absolute value of the t was greater than t corresponding with a probability of 0.05 [47].The coefficients of determination (R 2 and adjusted R 2 ) describe the model fit.A higher adjusted R 2 indicates a better model.The equation of adjusted R 2 is following (Equation ( 11)); where, n is the sample size and k is total number of regressors in the linear model (not counting the constant term).

Model Building and Accuracy at Phase I (ALS data)
The precisions of the ALS models are detailed in Table 5.Table A1 shows the model parameters and their statistical significance.Scatter plots are presented in Figure 5 to illustrate the residual plots for AGB.Although, we had a total of 46 ALS predictors initially in the set, the final AGB and volume models include only one ALS predictor.We also tested the transformation of the ALS predictors instead of the original value but it did not improve our model accuracy.The ratio of last and single echoes with height lower than 13.5 m and the total number of last and single echoes constituted the best model for AGB and volume.Table 5 shows the accuracy of the ALS prediction against field estimation using three combinations of training plots (50, 25 and 10 training plots).We reduced the number of training plots from 50 to 25 and 10 training plots (phase I) to assess the prediction accuracy after reducing the number of training plots.After sorting the training plots based on the tree height and AGB distribution, we selected every second plot to obtain 25 training plots, and every fifth plot to obtain 10 training plots out of 50 plots.The results in Table 5 proved that the models based on all three combinations of training plots had promising prediction accuracy.The linear model that was calibrated to estimate AGB showed high correlation when validated against field estimation using LOOV.A lowest RMSE value of 13 tons/ha (12%) was observed for AGB, while volume had similar RMSE value of 25 m 3 /ha (12%) for the models based on 10 training plots.In contrary, the highest 20% RMSE of AGB was employed for the models based on 25 training plots.In addition, the R 2 of the model was 0.83 based on 10 training plots.A slightly better R 2 value (0.85) was observed for volume.However, the models based on 25 training plots had a lowest R 2 value for both AGB (0.65) and volume (0.61).The statistical outliers were not frequent in the each residual plot (Figure 5).Although our final AGB and volume models each contain one ALS predictor, tests were also conducted by adding more ALS height and density predictors to the models.However, this did not improve the model accuracy significantly.For instance, we added the 80% ALS height percentile in the AGB model in addition to the used ALS predictor.However, the R 2 value dropped from 0.83 to 0.50.Similarly, we added the ratio of the number of last and single echoes below 5 m (low vegetation) and the total number of last and single echoes in the AGB model.Similarly, the R 2 value dropped to 0.66.Besides, the Pearson correlation between the used density metric and tree AGB, volume, tree height, basal area were respectively 0.84, 0.85, 0.82, 0.75, which confirms that this density metric explained very well the structure of study area.In contrast, the correlations between the ALS height percentiles of 60%, 70%, 80%, 90% and AGB were 0.60, 0.61, 0.62, and 0.64, respectively.In the case of volume, the correlations were higher for the ALS height percentiles and were respectively 0.68, 0.68, 0.69 and 0.70, which indicates that the ALS height percentiles correlate better with volume than AGB.The key of the used density metric must be the 13.5 m height threshold; it is much higher than what is typically used.13.5 m seems to be close to the smallest mean height at the plots.We think that it is a combination of a small variation in the plot data and the unusually high threshold which somehow makes this variable special for this particular data set.It might also be possible to achieve these results in a managed commercial forest, where an increase in height would also result in an increase in density, so one density variable would be enough.

Model Building and Accuracy at Phase II (RapidEye Data)
The error statistics of the linear models are given in Table 6.Table A1 shows the details of the regression models for RapidEye data.Figure 6 shows the residual plots for AGB.The predictor selection techniques picked up six and five predictors, respectively, for AGB and volume from a set of 36 independent RapidEye predictors.The best model for AGB constitute with the mean of four RapidEye spectral bands (green, red, red-edge, NIR), mean NDVI (based on red-edge and green) and the Haralick feature of inverse difference moment (HR5).Volume model consist of all AGB predictors except red band (mean).The mean NDVI and inverse difference moment accounted for much more of the variability in biomass and volume than the individual RapidEye bands.Notably, the RapidEye spectral bands were very helpful in increasing the variance explained.We calibrated the RapidEye models against ALS-trained set of surrogate plots using LOOV.The mean and standard deviation of AGB for the ALS-trained set of surrogate plots were 103 tons/ha and 34 tons/ha, respectively.We had a similar accuracy for all the combination of training plots number (10, 25, and 50).However, RapidEye prediction based on the 10 training plots against ALS estimation gave a lowest RMSE of 19 tons/ha (19%) for AGB.In addition, volume had a similar RMSE 38-41 m 3 /ha (20%-21%) value from all the combination of training plots (Table 6).The RapidEye predicted value showed a high correlation against the ALS estimates.The model for AGB explained 59% of the variability whereas the model for volume explained 58% (Table 6).The RapidEye model was not affected by the saturation effect (suppression of variance).Albeit, the AGB model did not show any heteroscedasticity in the model, the volume model had some affect but was not statistically significant.We used the logarithmic and quadratic transformation in an attempt to minimize this problem but they did not improve the model significantly.

Accuracy at the Independent Validation Plots
The accuracies at the validation plots are shown in Table 7. Figure 7 shows the residual plots comparing the RapidEye prediction against field estimation for AGB.The residual plots for AGB are reasonably good (Figure 7) despite a slight non-linearity in the case of volume.Table 7 shows that there was relatively similar accuracy for the validation plots using three combinations of training plots.However, the models based on the 10 training plots had slightly better accuracy for AGB, while volume prediction accuracy was highest in the models based on 50 training plots.A RMSE value of 23-24 tons/ha (20%) was observed AGB, while volume had a RMSE value of 43-44 m 3 /ha (19%-20%).We found that these accuracies at the independent validation plots can be considered acceptable and are even better than in studies concerning conventional compartments-based field inventory in Finland [48,49].All the biases were insignificant at the independent validation plots.It is worth noting that due to the small number of independent validation plots (n = 28), the imprecision in terms of the relative bias is no longer observed when predictions are plotted versus measured observations because the absolute scatter is approximately the same.In contrast, we also fitted the RapidEye model directly using 50 field plots, and used it to predict 28 validation plots.We found that the accuracy was slightly lower compared to the prediction with the use of surrogate plots (phase II).The relative RMSE and the adjusted R 2 were 21% and 44%, respectively, for AGB.

Discussion and Conclusion
This paper presented and tested forest AGB and volume estimation using ALS and RapidEye sensor data in a mixed-species boreal forest of eastern Finland.We found that the outcomes of the study were encouraging and we felt that they need further validation in other forest types along with development of the methodology.The overall strength of the ALS-RapidEye fusion revealed here is a promising accuracy to characterize biomass accounting in the coniferous forest ecosystems.The analysis of linear regression showed that the ALS data had good prediction accuracy at phase I.It was promising that the ALS models explained 83% (R 2 value of 83%) of the variability for AGB.Furthermore, the RapidEye models provided a relatively good accuracy at phase II for AGB and as well as for volume.RapidEye data had a relative RMSE of 20% at independent validation plots.Such accuracy indicates that the combination of ALS and RapidEye would be a promising fusion for the estimation of boreal forest attributes.Nevertheless, more testing and validation should be done in different forest landscapes at large-scale.

Statistical Modelling
We used the simple linear regression to predict biomass and volume which is one of the most commonly used methods in remote sensing-based regression modeling (e.g., [17,[50][51][52]).Naesset et al. [53] made a conclusion that employing a multivariate (multi-variable) method would not have had any significant impact on the biomass and volume estimation among the tested OLS regression analysis, seemingly unrelated regression (SUR), and partial least-squares (PLS) regression from two different inventories using ALS data.
Though we used regression subset selection techniques (leaps package) using an efficient branch-and-bound algorithm in R statistics, Packalén et al. [54] mentioned that different predictor variables were identified during different runs of the automated variable selection method.For instance, Latifi et al. [46], Packalén et al. [55] reported that genetic algorithm (GA) and simulated annealing (SA) were also superior variable selection methods.These algorithms should nevertheless be tested as alternatives to see whether a different explanatory variable selection method would improve the result or at least make the search for an eventual solution more effective.Nevertheless, the ALS predictor (density metric) that we employed in the biomass and volume models was closely matched to the predictors selected in other studies (e.g., Rana [56]).The RapidEye spectral and textural features used in this study are also closely matched with other studies (e.g., Packalén et al. [26] and Latifi et al. [46]).

Method Pros and Cons
Our presented forest biomass inventory method naturally places some advantages and constraints on the applicability of the method.Albeit, it involves a sample ALS data at phase I and requires only a few sample plots for model calibration (fitting).In addition, the surrogate plots (phase II) could be placed systematically in the whole area based on vegetation types, geographic location, climatic condition and tree species compositions.Therefore, our calibrated models at phase II represent the whole study area and have less chance of missing the vegetation types and tree species' compositions.The surrogate plots should be placed so that they reflect the full range of variation in biomass over the study area.In addition, the surrogate plots should cover also the rare forest types.Systematic selection or selection based on predicted AGB would be a natural choice.The locations of the surrogate plots could be selected through weighted random sampling or stratified sampling.Traditional forest inventory depends on a large set of ground-truth data for model calibration.However, this study showed a promising outcome where a sample ALS data was employed as a ground-truth data for model calibration with the satellite image for mapping the whole area of interest.
Another issue is that the satellite image acquisition dates may be different, introducing seasonal effects.Therefore, it would require the DN (digital number) values to be comparable between datasets which in turn would require the use of absolute reflectance instead of relative DN values.It would be recommended to use the same season for image acquisition.The result from Figures 2 and 3 confirmed that we did not need image normalization here.However, image normalization would be required for all satellite images relative to each other.In addition, the field data should always be collected at the same time period/season.In our study, the RapidEye images were collected almost two years after the field and ALS data acquisitions.However, it is acceptable in our study because there was no harvesting and logging activities and no significant naturally occurring changes within that time period.
Our study showed that ALS-based forest inventory have produced very accurate estimates (phase I).ALS data have high accuracy to predict forest biomass when regressing ALS height/density metrics with data from field measured plots [57].Subsequently, ALS estimates for the surrogate plots were used as simulated ground-truth (phase II) for the interpretation of optical satellite images, which improved the results compared with using models that were directly based on the field plots.Covering the whole area of interest with ALS is relatively expensive, thus we used a two-phase estimation approach that requires ALS data only from a sample of the study area.This is known as the ALS-assisted multi-source programme that combines ALS information with field plots and satellite data to develop a forest biomass map [58].
Reducing the number of training plots is of great importance to minimize field inventory cost.Our findings show that a reduced number of sample plots provides similar accuracy compared to the full dataset.When we reduced the sample plot from 50 to 10, there was no significant reduction in the prediction accuracy of AGB and volume (Table 7).Junttila et al. [59] reported similar performances in the boreal managed forest area, Finland.They concluded that the effect of a reduction in sample-plot number had only a slight decrease in accuracy and a reduction of 2.9% RMSE for the reduced sample-plot number compared to the full dataset.Latifi and Koch [60] thoroughly studied the integration of off-site samples into small-scale forest inventory in central Europe, in which they reported the substantial effectiveness of low-to-medium sampling intensity on improving forest structure models by means of ALS and multispectral aerial images.Dalponte et al. [61] concluded that it is possible to reduce the number of field sample plots without losing model accuracy.A bottleneck of their approach was that ALS data is necessary as prior information before collection of field reference data.They reduced the number of sample plots based on genetic algorithms, whereas our selections were based on tree height and AGB distribution.

Model Calibration and Validation
The ALS models were calibrated against the field sample plots in phase I of the model building.The R 2 value of our predicted forest AGB is close to the others studies ranged between 0.74 and 0.88 in the studies by Hall et al. [62], Bright et al. [63], Fu et al. [64], Gobakken et al. [65], Naesset and Gobakken [45], Hawbaker et al. [66] and Ferster et al. [67].Subsequently, in the phase II, we calibrated the RapidEye models against the ALS-simulated estimation.When the resulting models at phase II were tested against an independent dataset, the error statistics of testing were similar to those of the model fitting (Tables 6 and 7).Our models are stable and more promising considering the error statistics.Naesset [68], Gómez et al. [69], Eckert [70], and Gautam [22] mentioned that optical imagery based estimation may lead to significant underestimation of biomass stock in areas with high (>200 ton/ha) AGB concentrations.In western Finland, Muukkonen and Heiskanen [71] had faced a similar problem in mixed special boreal forest employing optical sensor alone where biomass was underestimated above 200 ton/ha.Our models and vegetation index were not affected by the saturation problem (underestimation of variance) because the biomass stock at the study area is mostly less than 200 ton/ha, although many authors have reported this problem for optical sensors.However, it would be interesting to explore the accuracy of our presented approach in other landscapes for wider acceptance based on the geographic variation and different forest types and environments.
The RMSE at the independent validation plots was lower than in most previous studies of boreal coniferous and mixed forest.Our study had a relative RMSE of 19%-20% for volume, whereas previous studies using the optical sensor had a relative RMSE of 42%-82% for volume, for instance, studies by Tokola and Heikkilä [72] (82%), Tomppo et al. [73] (59%), Kilpeläinen and Tokola [74] (56%), Hyyppä et al. [75] (50%), Mäkelä and Pekkarinen [76] (47%), and Hyvönen [77] (42%).Similar to the volume, the relative RMSE (20%) for AGB was also noticeably lower in the present study than the corresponding statistics of Tomppo et al. [73] (36%), and Muukkonen and Heiskanen [71] (41%).It is worth noting that the above studies did not follow exactly the same approach of our study but are the closest available comparison.The biases were considerably higher for testing datasets compared to the model fitting, though the biases were not statistically significant.Due to the small number of testing plots we could not make a clear conclusion about biases here.

Sensor Pros and Cons
The relatively small errors in this study are probably due to the type of biomass being measured and the relatively good spatial and spectral resolution of the RapidEye data with high geometric accuracy.Hyyppä et al. [75] also mentioned that the estimation error of the forest characteristics were smaller for higher resolution data compared to the lower resolution data.Good spatial resolution of the RapidEye bands makes it is a good data source for studies of forest characteristics.Koch [78], Treitz et al. [79], and Kankare et al. [80] mentioned that the mapping of forests biomass, from local to global in terms of their status and development, is of increasing importance for REDD.

Concluding Remarks
We examined here the AGB and volume estimation employing ALS and RapidEye data in a two-phase sampling method.Linear regression was employed to predict forest characteristics in a managed boreal forest in eastern Finland.To conclude, the present study has confirmed that the accuracy of AGB and volume comparable to the forest inventory by compartments (or stand) in Finland.The approach presented here is a promising alternative for use in forest management in Finland with enough accuracy for the purpose of forest resource inventory.In addition, calibrating the RapidEye data using the ALS-simulated surrogate plots and different models at phases I and II make this approach promising and stable for biomass and carbon accounting in the boreal forest.It could also offer valuable methodology for inventories need in the REDD program.As the tested area is rather small, further validation is needed in a larger study area for better justification of the presented approach.Finally, we feel that future studies should be carried out using a combination of ALS and other optical sensors to confirm the validity of this approach in other forests landscape.This approach could be tested by tracking changes in forest biomass at the regional and national levels.

Figure 1 .
Figure 1.Training, surrogate and validation plots location and administrative map of Finland (left side), and RapidEye image with study area marked (right side).

Figure 4 .
Figure 4. Flowchart of the two-phase sampling design.

Figure 5 .
Figure 5.The residuals plots of ALS predicted AGB at phase I.

Figure 6 .
Figure 6.The residuals plots of RapidEye for AGB at phase II.

Figure 7 .
Figure 7.The residual plots of RapidEye for AGB at validation plots.
applied a multivariate mixed linear model in biomass predictions and the equations are as follows:

Table 2 .
The mean characteristics of the training plots and validation plots.

Table 3 .
List of airborne laser scanning (ALS) explanatory predictors.
Height for which the cumulative sum of ordered first and single echo heights is closest to 10%, 20%, 30%...100% of the total height sum.11...20 H lp10-100 Height for which the cumulative sum of ordered last and single echo heights is closest to 10%, 20%, 30%...100% of the total height sum.21...23 Ifp 30-90 Intensity for which the cumulative sum of ordered first and single echo intensities is closest to 30%, 60% and 90% of the total intensity sum.24...26 Ilp 30-90 Intensity for which the cumulative sum of ordered last and single echo intensities is closest to 30%, 60% and 90% of the total intensity sum.27 H mean Mean height of first and single echo vegetation points (points over high vegetation threshold 5 m).28 H std Standard deviation of first and single echo heights.29 Dfp Ratio of the number of first and single echoes below 5 m (low vegetation) and the total number of first and single echoes.30 Dlp Ratio of the number of last and single echoes below 5 m (low vegetation) and the total number of last and single echoes.

Table 5 .
Accuracy of the ALS prediction against field estimation.

Table 6 .
Accuracy of RapidEye prediction against ALS estimation at phase II.
AGB *: above ground biomass; ǂ The number of training plots from phase I used in the estimation of 200 ALS surrogate plots for phase II.R 2 Adj: R 2 adjusted.

Table 7 .
Accuracy of AGB and volume obtained from RapidEye at validation plots.
AGB *: above ground biomass; ǂ The number of training plots used in the estimation of 200 ALS surrogate plots at phase II.R 2 Adj: R 2 adjusted.

Table A1 .
Cont. : above ground biomass; D hlp4 : the ratio of last and single echoes with height lower than 13.5 m and the total number of last and single echoes; B2: green, mean; B3: red, mean; B4: red-edge, mean; B5: NIR, mean; NDVI 2 : red-edge and green band based NDVI, mean; HR5: inverse difference moment.©2013 by the authors; licensee MDPI, Basel, Switzerland.This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/). AGB