Evaluation of Radiometric and Atmospheric Correction Algorithms for Aboveground Forest Biomass Estimation Using Landsat 5 TM Data

Solar radiation is affected by absorption and emission phenomena during its downward trajectory from the Sun to the Earth’s surface and during the upward trajectory detected by satellite sensors. This leads to distortion of the ground radiometric properties (reflectance) recorded by satellite images, used in this study to estimate aboveground forest biomass (AGB). Atmospherically-corrected remote sensing data can be used to estimate AGB on a global scale and with moderate effort. The objective of this study was to evaluate four atmospheric correction algorithms (for surface reflectance), ATCOR2 (Atmospheric Correction for Flat Terrain), COST (Cosine of the Sun Zenith Angle), FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) and 6S (Second Simulation of Satellite Signal in the Solar), and one radiometric correction algorithm (for reflectance at the sensor) ToA (Apparent Reflectance at the Top of Atmosphere) to estimate AGB in temperate forest in the northeast of the state of Durango, Mexico. The AGB was estimated from Landsat 5 TM imagery and ancillary information from a digital elevation model (DEM) using the non-parametric multivariate adaptive regression splines (MARS) technique. Field reference data for the model training were collected by systematic sampling of 99 permanent forest growth and soil research sites (SPIFyS) established during the winter of 2011. The following predictor variables were identified in the MARS model: Band 7, Band 5, slope (β), Wetness Index (WI), NDVI and MSAVI2. After cross-validation, 6S was found to be the optimal model for estimating AGB (R2 = 0.71 and RMSE = 33.5 Mg ̈ha ́1; 37.61% of the average stand biomass). We conclude that atmospheric and radiometric correction of satellite images can be used along with non-parametric techniques to estimate AGB with acceptable accuracy.


Introduction
Apart from sensor gains/offsets, solar irradiance and Sun-Earth geometry, the characteristics of electromagnetic energy detected by remote sensing optical sensors is affected by particles and other components present in the atmosphere.Hence, as solar radiation passes through the atmosphere (Sun-surface-sensor), it is affected by absorption or scattering by particles in suspension (aerosols) and other atmospheric elements, thus creating a hazy effect that distorts the radiometric properties of satellite images [1][2][3][4].Various algorithms have been developed to correct such effects in the image.These include calibrations based on sensor parameters, solar-Earth geometry, dark object subtraction and radiative transfer [5][6][7][8].Correction of atmospheric effects is important in relation to improving data quality [9,10].However, the use of correction algorithms to minimize atmospheric effects, especially the scattering and absorption caused by aerosols, remains challenging [11].This particularly applies to the parameterization of algorithms for calculating the surface reflectance for its eventual use in estimating aerial or aboveground forest biomass (AGB) [12].
AGB is an important parameter that provides information about the current status of forests and therefore facilitates forest management decisions [13,14].The use of optical remote sensors of medium spatial resolution, such as the Landsat Thematic Mapper (TM), has become increasingly common in the last few decades.The use of these sensors for large-scale monitoring of AGB dynamics, particularly in the case of multi-temporal/multi-scene analyses or unevenly-distributed atmospheric effects in single-image analyses, frequently involves the application of different atmospheric correction algorithms.Such algorithms correct distortions between or within images other than those related to real land cover differences [15][16][17][18][19].
The aim of the present study was to evaluate four different atmospheric correction algorithms (surface reflectance)-COST (Cosine of the Solar Zenith Angle), a modification of the dark object subtraction method [20] with the addition of a multiplicative correction for atmospheric transmittance [21], ATCOR2 (Atmospheric Correction for Flat Terrain), [22]), FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes, [23,24]) and 6S (Second Simulation of Satellite Signal in the Solar Spectrum, [25]), and one radiometric correction algorithm (apparent reflectance at sensor), ToA (Top of Atmosphere), for estimating AGB from spectral data captured by the Landsat 5 TM sensor and topographic variables derived from a digital elevation model (DEM) generated using data collected by sampling of the permanent forest growth and soil research sites (SPIFyS) established in a temperate forest in the northwest of the state of Durango, Mexico.
Different methods could be used for remote AGB estimation in forests, including parametric and nonparametric approaches.However, since forest structure and biomass often entail nonlinear variability, variable interaction across scales and autocorrelation, nonparametric approaches often markedly outperform parametric methods [26].
In this study, we investigate remote AGB estimation of mixed and uneven-aged forests using the multivariate adaptive regression splines (MARS).To the best of our knowledge, MARS have rarely been used for remote estimation of AGB [27][28][29], but in these cases, the method has performed best for prediction than other parametric and non-parametric approaches, such us linear models, classification, regression trees and artificial neural networks [27] or hybrid tree-based algorithms [28,29].

Study Area
The study area is located in the northwest of the state of Durango (Figure 1).The forest vegetation comprises pine (Pinus spp.), oak (Quercus spp.), fir (Abies spp.) and plurispecific stands with different proportions of species of the genera Pinus and Quercus, in accordance with the Land and Vegetation Use Map, scale 1:250,000, Series V [30].The weather is cold in the canyons (with temperatures ranging from ´20 ˝C in winter to 20 ˝C in summer) and mild or warm in the valleys (with temperatures ranging from 10 ˝C in winter to 40 ˝C in summer).

Field Data
The field data were AGB values taken from a database of stand variables measured in 99 permanent forest growth and soil research sitesSPIFyS established during the winter of 2011 in accordance with the method developed by Corral-Rivas et al. [31].The SPIFyS are square plots of a size of 50 × 50 m (surface area, 2500 m 2 ).The location of the sites was carried out by systematic sampling at the regional level using a squared grid of five kilometers, although a small squared grid (3 km) was used on very steep areas to include their variability.
In each sample plot, the tree species were recorded, and the diameter at breast height (cm) and total height (m) of all standing trees were measured.Species-specific models developed for the study area were used to estimate individual tree aboveground biomass [32].The goodness of fit for such statistical models ranged from 0.87-0.99(R 2 ), and the root mean square error (RMSE) varied from 22.8-95.2kg.Once the tree aboveground biomass was estimated, all values from each sampling plot were summed and expressed on a per hectare basis.Table 1 shows the descriptive statistics of the main stand variables.

Spectral Data from the Landsat 5 TM Sensor
The image used in the study, which was captured by the Landsat 5 TM sensor in April 2011, is available at the U.S. Geological Service (USGS) website [33].The sensor operates in seven bands of the electromagnetic spectrum: Band 1 (blue, 450-520 nm), Band 2 (green, 520-600 nm), Band 3 (red,

Field Data
The field data were AGB values taken from a database of stand variables measured in 99 permanent forest growth and soil research sitesSPIFyS established during the winter of 2011 in accordance with the method developed by Corral-Rivas et al. [31].The SPIFyS are square plots of a size of 50 ˆ50 m (surface area, 2500 m 2 ).The location of the sites was carried out by systematic sampling at the regional level using a squared grid of five kilometers, although a small squared grid (3 km) was used on very steep areas to include their variability.
In each sample plot, the tree species were recorded, and the diameter at breast height (cm) and total height (m) of all standing trees were measured.Species-specific models developed for the study area were used to estimate individual tree aboveground biomass [32].The goodness of fit for such statistical models ranged from 0.87-0.99(R 2 ), and the root mean square error (RMSE) varied from 22.8-95.2kg.Once the tree aboveground biomass was estimated, all values from each sampling plot were summed and expressed on a per hectare basis.Table 1 shows the descriptive statistics of the main stand variables.The image used in the study, which was captured by the Landsat 5 TM sensor in April 2011, is available at the U.S. Geological Service (USGS) website [33].The sensor operates in seven bands of the electromagnetic spectrum: Band 1 (blue, 450-520 nm), Band 2 (green, 520-600 nm), Band 3 (red, 630-690 nm), Band 4 (near-infrared, 760-900 nm), Bands 5 and 7 (mid-infrared region, 1550-2350 nm) and Band 6 (which provides information in the thermal infrared region and is not used in this type of study).Likewise, the spectral indices Normalized Difference Vegetation Index (NDVI, [34]) and Modified Soil Adjusted Index (MSAVI2, [35]) were used.Such indexes are potentially less sensitive than single band values to artefacts due to differences in light conditions, exposed soil or terrain slope and orientation, which might eventually affect AGB estimation.

Radiometric Correction Algorithms
A useful prior step to the interpretation of satellite images is converting the digital numbers (DNs) stored in the original image into biophysical variables of standard significance (reflectance).These variables are comparable in the same sensor over time and over scenes, as well as between different sensors and between remote sensing and other methods of detecting electromagnetic energy.The correction is also advisable in the case of unevenly-distributed atmospheric effects in the image and also when vegetation indexes based on band ratios are included in the analyses [36][37][38][39][40][41].In order to calculate surface reflectance, the atmospheric effects must be removed.This involves estimating the atmospheric transmissivity (both up and down welling), diffuse irradiance and atmospheric radiance due to scattering [42].In the present study, we evaluated five different radiometric correction methods: the first four methods (namely ATCOR2, COST, FLAASH and 6S) including atmospheric corrections (surface absolute reflectance), while the fifth (ToA) without considering atmospheric effects (apparent reflectance at the sensor) for the final estimation of AGB.

Atmospheric Correction for Flat Terrain (ATCOR2)
The aim of ATCOR2 correction [22] is to remove the atmospheric effects in order to recover the physical parameters of the Earth's surface, including the surface reflectance, soil visibility and the temperature.The correction is carried out with the ATCOR module [43] included in ERDAS IMAGINE software, Version 2013 [44].

Cosine of the Sun Zenith Angle (COST)
This is a radiometric calibration method that considers the atmospheric effect and is based entirely on the characteristics of the satellite image, in contrast with other methods of atmospheric correction, like ATCOR2, FLAASH or 6S, requiring some extra parameters, such as atmospheric profiles, the aerosol models or visibility [21,45].COST applies Dark Object Subtraction (DOS, [20]) to compensate for the additive components of the atmosphere, which mainly affect the shortest wavelengths.DOS does not take into account the multiplicative effect on the longer wavelengths.For initial estimation of the multiplicative effect, the value of the atmospheric transmittance along the path from the ground to the sensor (TAUz) was computed from the cosine of the solar zenith angle [46].This correction was carried out by implementing the algorithm in the Model Maker ® module of ERDAS IMAGINE software, Version 2013 [44].

Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH)
This is an advanced atmospheric correction module based on the MODTRAN4 algorithm of radiative transfer developed by Spectral Sciences Inc. (Burlington, MA, USA) under the sponsorship of the U.S. Air Force Research Laboratory [23,24].The technique is initially based on the standard equation of spectral radiance for each pixel of the sensor, which applies to the range of wavelengths of solar radiation (thermal emission is omitted) and a Lambertian and flat surface or their equivalents.It considers the radiance reflected from the Earth's surface and scattered by the atmosphere towards the sensor.The difference between these two radiances is due to the adjacency effect (spatial mixture of radiance among nearby pixels) caused by atmospheric scattering.The correction was carried out with the FLAASH module in ENVI ® software, Version 5.1 [47].

Second Simulation of Satellite Signal in the Solar Spectrum (6S)
This procedure eliminates atmospheric effects on the reflectance values in images captured by sensors on board satellite or aircraft platforms.It is based on an advanced code of radiative transfer, designed for simulating the interaction between solar radiation and an atmospheric-surface system, together with a wide range of atmospheric, spectral and geometric conditions [25].The code acts on the basis of the successive order of scattering (SOS) method and explains the polarization of radiation in the atmosphere by calculating the different components.The scenes used to belong to the National Landsat Archive Processing System (NLAPS) and correspond to the product obtained by Landsat 4-5 Thematic Mapper Level 1 of surface reflectance processed by the Standard Landsat Product Generation System (LPGS) and using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (available at the U.S. Geological Service website [33]).

Apparent Reflectance at the Top of Atmosphere (ToA)
This technique enables calculation of the apparent reflectance in a satellite image and consists of converting the DNs to radiance values and then to reflectance values.The word "apparent" refers to the fact that the reflectance has not been corrected for atmospheric effects and represents an initial normalization of the image [48].The correction was carried out with the "apparent reflectance" method implemented in the "landsat" package [49] in the R software [50].

Parameters Derived from the Digital Elevation Model (DEM)
The primary and secondary terrain parameters were calculated from the DEM (Table 2), after the application of a low bandpass filter with a moving window of 5 ˆ5, with the aim of reducing the banding effect present in the original archive [51].These parameters are considered as predictor variables for estimating AGB, because they are directly related to the distribution and development of forest species and can subsequently be evaluated [52,53].

Variable Formula Reference
Elevation Digital Elevation Model where: Z, average elevation; R, point radio altitude units; As, drainage area specified; tanβ, local slope angle; D, F, G and H were derived according to the equation of [56].

Generation of a Database
The database used to estimate AGB was generated by extracting the mean values of the satellite image pixels after each of the radiometric correction techniques and of the different DEM parameters by considering a buffer of 25 m for the geolocalization of the SPIFyS plots.The values were extracted using the "raster" package [57], implemented in R statistical software [50].

Analysis of Variance
In order to establish whether there were any significant differences between the radiometric correction algorithms considered, analysis of variance (ANOVA) was used to compare the values of each spectral band of the sensor obtained by processing the images of the SPIFyS and applying the different correction algorithms.The correction algorithm was used as a fixed factor: where y i are the values of each spectral band of the sensor corrected by the algorithm i; CA i is the correction algorithm used (fixed factor) and ε i are the associated errors.
Prior to the application of the ANOVA, the dependent variables normality and homogeneity of variance were checked using respectively Shapiro-Wilk's test [58] and Levene's test [59].When the results of these tests indicated that the ANOVA assumptions were not satisfied, the dependent variables were transformed by a natural log transformation.When the ANOVA indicated significant differences between correction algorithms (α = 0.05), Tukey's multiple comparison test was used to identify which algorithms were different.Likewise, box plots were elaborated for each variable (spectral band) with the aim of facilitating graphical interpretation of the results.All analyses were carried out with R statistical software [50].

Multivariate Adaptive Regression Spline (MARS)
The statistical analysis was carried out by the non-parametric multivariate adaptive regression splines (MARS) proposed by Friedman [60], using the "earth" package [61] implemented in the R software [50].This method involves constructing a non-linear regression model based on a product of functions denominated "smoothed basis functions" (splines).The predictors are incorporated in their structure as part of a function generating a model for the dependent variable, which may be continuous or binary, and at the same time automatically selects the predictors in the final model.The general form of the MARS non-parametric regression model, formulated on the dependent "y" variable and the "x" predictors, is as follows: where ε is the error and f (x) is the unknown regression function, derived as follows: where β 0 is the intercept of the model, B m (x) is the m-th basis function, βm is the coefficient of the m-th basis function and M is the number of basis functions in the model.Each basis function B m (x) takes one of the following two forms: (i) a hinge function of the form max (0, x-k) or max (0, k-x), where k is a constant value; (ii) a product of two or more hinge functions that, therefore, can model the interaction between two or more predictors (x).The optimal model was selected using the backward procedure.An overfitted model was generated with all possible predictor variables, and the model was pruned by removing terms one by one, deleting the least effective term at each step until the best submodel was found.Model subsets were compared using the generalized crossed validation (GCV) criterion.GCV is an approximation to the error that would be determined by leave-one-out validation and considers both the residual error and the model complexity.Therefore, the optimal model will be that yielding the lowest GCV.GCV pMq " where y i are the observed values of the dependent variable, n is the number of observations, f M px i q is the MARS model with M basis functions, x i are the observed values of the predictors included in the MARS model and p M is the number of parameters of the MARS model.
To analyze the importance of the independent variables that contribute most to the final model, three different criteria were used: (i) nsubsets, i.e., the number of times that each variable is included in a subset (in the final model); (ii) sqr_rss, i.e., the decrease in the sum of square errors (RSS) for each subset relative to the previous subset; and (iii) sqr_gcv, i.e., the decrease in the GCV value for each subset relative to the previous subset.These criteria are normalized on a scale of 100, to facilitate the interpretation of the contribution of each predictive variable in the model.
With the aim of evaluating the performance of the fitted model for each of the radiometric correction algorithms, the following statistics were determined to compare the AGB reference data for the monitoring plots with the estimates based on the images subjected to different corrections: root mean squared error (RMSE), the coefficient of determination (R 2 ) and the generalized coefficient of determination (GR 2 ): GR 2 " 1 ´GCV GCV null (7) where y i , ŷi and y are, respectively, the observed, estimated and mean values of the dependent variable; n is the total number of observations used to fit the model; p is the number of parameters to estimate; and GCV null is the GCV of a model with a single independent term (β 0 ).Because the quality of the fit does not necessarily reflect the quality of the prediction, assessment of the validity of the model with an independent dataset is desirable [62].Therefore, cross-validation was carried out by splitting the dataset into different numbers of folds (2-10), and the goodness-of-fit statistics (Equations ( 5)-( 7)) were calculated for each fold (test data) by using the models fitted to the other folds (training data).

Generation of Thematic Maps
The final stage of the overall process illustrated was the calculation of AGB maps.These maps were derived by implementing the basic functions resulting from the MARS model fit for each of the atmospheric correction algorithms by using the raster calculator in ArcGIS ® software [63].
An overview of the workflow described in the previous sections is presented in Figure 2.
Remote Sens. 2016, 8, 369 7 of 18 subset relative to the previous subset.These criteria are normalized on a scale of 100, to facilitate the interpretation of the contribution of each predictive variable in the model.With the aim of evaluating the performance of the fitted model for each of the radiometric correction algorithms, the following statistics were determined to compare the AGB reference data for the monitoring plots with the estimates based on the images subjected to different corrections: root mean squared error (RMSE), the coefficient of determination (R 2 ) and the generalized coefficient of determination (GR 2 ): ( ) where yi, ŷi and are, respectively, the observed, estimated and mean values of the dependent variable; n is the total number of observations used to fit the model; p is the number of parameters to estimate; and GCVnull is the GCV of a model with a single independent term (β0).
Because the quality of the fit does not necessarily reflect the quality of the prediction, assessment of the validity of the model with an independent dataset is desirable [62].Therefore, cross-validation was carried out by splitting the dataset into different numbers of folds (2-10), and the goodness-of-fit statistics (Equations ( 5)-( 7)) were calculated for each fold (test data) by using the models fitted to the other folds (training data).

Generation of Thematic Maps
The final stage of the overall process illustrated was the calculation of AGB maps.These maps were derived by implementing the basic functions resulting from the MARS model fit for each of the atmospheric correction algorithms by using the raster calculator in ArcGIS ® software [63].
An overview of the workflow described in the previous sections is presented in Figure 2.

Results
In the first step, a visual comparison of the application of the different radiometric correction algorithms to the Landsat 5 TM image of the study area was done.All atmospheric correction algorithms were parameterized using the same type of "rural" aerosol, for zones not affected by urban or industrial sources and a visibility larger than 40 km.In comparison with the original image in DNs, the brightness of the other images differs due to the parameterization used in each algorithm to correct the radiometric and atmospheric effects (Figure 3).

Results
In the first step, a visual comparison of the application of the different radiometric correction algorithms to the Landsat 5 TM image of the study area was done.All atmospheric correction algorithms were parameterized using the same type of "rural" aerosol, for zones not affected by urban or industrial sources and a visibility larger than 40 km.In comparison with the original image in DNs, the brightness of the other images differs due to the parameterization used in each algorithm to correct the radiometric and atmospheric effects (Figure 3).The mean reflectances of the 99 SPIFyS plots obtained for each band of the Landsat 5 TM sensor and radiometric correction algorithm were graphically compared (Figure 4).The differences vary according to the atmospheric effects corrected by each particular algorithm.In general, the spectral signature of the vegetation showed in all cases typical behavior: low reflectance in the bands of the visible spectrum (450-690 nm) and high reflectance in the near-and mid-infrared bands (760-2350 nm).However, a clear difference appears depending on whether the radiometric correction algorithms were applied to variables with physical meaning (e.g., ground reflectance) rather than to DNs.In this respect, the ATCOR2 and FLAASH algorithms overcorrected the atmospheric effects (underestimation of transmissivity) in the short wavelength regions and underestimated them in the near-and mid-infrared regions, while COST and ToA produced the inverse spectral pattern.Finally, the performance of 6S was intermediate between that of the other algorithms evaluated.The mean reflectances of the 99 SPIFyS plots obtained for each band of the Landsat 5 TM sensor and radiometric correction algorithm were graphically compared (Figure 4).The differences vary according to the atmospheric effects corrected by each particular algorithm.In general, the spectral signature of the vegetation showed in all cases typical behavior: low reflectance in the bands of the visible spectrum (450-690 nm) and high reflectance in the near-and mid-infrared bands (760-2350 nm).However, a clear difference appears depending on whether the radiometric correction algorithms were applied to variables with physical meaning (e.g., ground reflectance) rather than to DNs.In this respect, the ATCOR2 and FLAASH algorithms overcorrected the atmospheric effects (underestimation of transmissivity) in the short wavelength regions and underestimated them in the near-and mid-infrared regions, while COST and ToA produced the inverse spectral pattern.Finally, the performance of 6S was intermediate between that of the other algorithms evaluated.As the Shapiro-Wilk's test indicated that the corrected values of the spectral bands of the sensor were not normally distributed, therefore a log transformation was used.The ANOVA applied to the data from the 99 SPIFyS plots detected significant differences for all algorithms in the visible region (Figure 5), except for the COST and FLAASH algorithms in Band 3. Significant differences were observed in the near-infrared region for all algorithms tested, except ATCOR2 and 6S, which were included in the same population group (Tukey's test); they also shared characteristics with FLAASH and ToA.Finally, in the bands corresponding to mid-infrared regions, the algorithms ATCOR2, FLAASH and 6S did not provide significantly different results (p > 0.05).
The importance of the predictor variables (spectral bands and indices and topographic variables) for each of the MARS models fitted to the different radiometric correction algorithms was calculated (Figure 6).Results pointed out the near-infrared as the most relevant spectral region for AGB prediction.In fact, among the predictor variables, Band 7 made the greatest contribution to the capacity for predicting AGB for the different algorithms analyzed in almost all the cases except in algorithm 6S, where the variable Band 5 appeared as the most important in accordance with the number of times that each variable is included in a subset.
Table 3 shows the results of applying the MARS technique for estimating the AGB on the basis of different correction algorithms.The algorithm that yielded the best predictions was ToA (R 2 = 0.89, RMSE = 29.82Mg•ha −1 ; 33.49% of the average stand biomass).However, after the application of cross-validation, the generalized coefficient of determination (GCV) decreased to 0.68.By contrast, the algorithm with the lowest predictive power was ATCOR2 (R 2 = 0.75 and GR 2 = 0.59).The algorithm with the highest generalized coefficient of determination was 6S (GR 2 = 0.71).As the Shapiro-Wilk's test indicated that the corrected values of the spectral bands of the sensor were not normally distributed, therefore a log transformation was used.The ANOVA applied to the data from the 99 SPIFyS plots detected significant differences for all algorithms in the visible region (Figure 5), except for the COST and FLAASH algorithms in Band 3. Significant differences were observed in the near-infrared region for all algorithms tested, except ATCOR2 and 6S, which were included in the same population group (Tukey's test); they also shared characteristics with FLAASH and ToA.Finally, in the bands corresponding to mid-infrared regions, the algorithms ATCOR2, FLAASH and 6S did not provide significantly different results (p > 0.05).
The importance of the predictor variables (spectral bands and indices and topographic variables) for each of the MARS models fitted to the different radiometric correction algorithms was calculated (Figure 6).Results pointed out the near-infrared as the most relevant spectral region for AGB prediction.In fact, among the predictor variables, Band 7 made the greatest contribution to the capacity for predicting AGB for the different algorithms analyzed in almost all the cases except in algorithm 6S, where the variable Band 5 appeared as the most important in accordance with the number of times that each variable is included in a subset.
Table 3 shows the results of applying the MARS technique for estimating the AGB on the basis of different correction algorithms.The algorithm that yielded the best predictions was ToA (R 2 = 0.89, RMSE = 29.82Mg¨ha ´1; 33.49% of the average stand biomass).However, after the application of cross-validation, the generalized coefficient of determination (GCV) decreased to 0.68.By contrast, the algorithm with the lowest predictive power was ATCOR2 (R 2 = 0.75 and GR 2 = 0.59).The algorithm with the highest generalized coefficient of determination was 6S (GR 2 = 0.71).Figure 7 shows the variations in the goodness of fit statistics (R 2 and GR 2 ) obtained by application of the cross-validation technique (nfold = 2-10) to the MARS model fits to each of the  Figure 7 shows the variations in the goodness of fit statistics (R 2 and GR 2 ) obtained by application of the cross-validation technique (nfold = 2-10) to the MARS model fits to each of the algorithms considered.In this case, the COST algorithm produced the best results and the highest stability, followed by ToA and 6S.
Remote Sens. 2016, 8, 369 12 of 18 algorithms considered.In this case, the COST algorithm produced the best results and the highest stability, followed by ToA and 6S.Finally, AGB maps of the study area were generated by implementing the fitted MARS models for each radiometric correction algorithm on the satellite image using the raster calculator with ArcGIS ® software [63] and then intersected with the vector layer of vegetation cover [30] (Figure 8).Finally, AGB maps of the study area were generated by implementing the fitted MARS models for each radiometric correction algorithm on the satellite image using the raster calculator with ArcGIS ® software [63] and then intersected with the vector layer of vegetation cover [30] (Figure 8).algorithms considered.In this case, the COST algorithm produced the best results and the highest stability, followed by ToA and 6S.Finally, AGB maps of the study area were generated by implementing the fitted MARS models for each radiometric correction algorithm on the satellite image using the raster calculator with ArcGIS ® software [63] and then intersected with the vector layer of vegetation cover [30] (Figure 8).

Discussion
The use of radiometric correction algorithms enabled the transformation of DNs to reflectance values (i.e., variables with biophysical meaning).The spectral reflectance of the forest cover in the study area was low (<0.10) in the visible region being typical of vegetated land (as chlorophyll absorbs most of the light received from the visible spectrum).Accordingly, the reflectance was higher (around 19%) in the near-infrared region, indicating the contrast between these regions of the electromagnetic spectrum [42,64] typical of green vegetation.
In accordance with [42], reflectance in the blue band (450-520 nm) was lower relative to the DNs and was distorted by the ToA algorithm, which assumes a transparent atmosphere over flat land and perfectly Lambertian surfaces.However, with the ATCOR2 algorithm, overcorrection at the shortest wavelengths in the visible region (i.e., green or red bands) is a consequence of scattering caused by particles of ozone and water vapor present in the atmosphere [65,66].According to Broszeit and Ashraf [67], the images corrected using the COST algorithm are less accurate than those corrected with the ATCOR2 algorithm for vegetation cover, because COST automatically selects the lowest pixel values for each band to eliminate atmospheric haze from the data, whereas ATCOR2 uses specific parameters of atmospheric geometry and of the sensor for improved conversion of the reflectance data.
In the present study, the 6S algorithm performed more consistently than the other algorithms evaluated (ATCOR2, COST, FLAASH and ToA), because of the capacity of 6S to correct the effects of water vapor, high temperature and, therefore, high aerosol levels.This finding is similar to that reported by Nazeer et al. [68], who observed that, among five algorithms tested, 6S produced the smallest difference in field-measured surface reflectance and that obtained using the Landsat ETM sensor.
The analysis of variance applied to data from the 99 SPIFyS plots revealed similar behavior of the population means for the ATCOR2, FLAASH and 6S algorithms in the bands of the mid-infrared region.However, the bands in the visible and near-infrared regions indicated a statistically-significant difference in the configuration of the groups of algorithms.These findings are consistent with those reported by Mahiny and Turner [69], who compared the first four bands of the Landsat TM sensor under five radiometric correction methods (two relative approaches: pseudo invariant features (PIF) and radiometric control sets (RCS); and three absolute approaches: COST, 6S and ToA), showing that in most cases, the four bands produced significantly different results.
Regarding the use of the MARS technique for estimating AGB in the present study, the ToA algorithm initially showed the greatest predictive power (R 2 = 0.89, RMSE = 29.8Mg¨ha ´1; 33.49% of the mean biomass in the stand).This result showed a slightly better performance than that reported by Hamdan et al. [70] who estimated AGB by fitting regression models to reflectance data (corrected by ToA) from the SPOT-5 and ALOS PALSAR (R 2 = 0.803; RMSE = 32.6 Mg¨ha ´1).Furthermore, in a study in southwest Colorado (USA), a lower correlation for biomass prediction (r = 0.86, RMSE = 45.6 Mg¨ha ´1) with data corresponding to vegetation indices and texture analysis in Landsat TM images corrected by ToA radiometric correction was obtained [18].However, further analyses in the present study, such as the use of the generalized error of cross-validation (GCV), indicated a GR 2 value of 0.68, pointing out a decrease in the predictive capacity of the model.Therefore, after application of GCV, the 6S algorithm proved to be optimal (GR 2 = 0.71) despite the potentially negative effect of considering a larger number of predictors (12 of 17) than in ToA correction.It should be taken into account that the GCV criterion considers both the residual error and also the model complexity, penalizing those models with a high number of parameters.The results are similar to those obtained by Nguyen et al. [71] in a study carried out in South Korea in which it was confirmed that the 6S algorithm fitted by the kNN technique (RMSE = 22.5 Mg¨ha ´1) produced better results for estimating AGB than the other algorithms tested (DOS, FLAASH and ToA), especially for Landsat ETM bands in the infrared region.Furthermore, various studies have concluded that the MARS technique is a flexible method that yields robust predictions [72,73].The scatterplot of observed versus predicted aboveground forest biomass obtained with the MARS model fitted to the spectral data corrected using the proposed algorithm (6S) is shown in Figure 9.No trends to under-or over-estimate were observed.
The scatterplot of observed versus predicted aboveground forest biomass obtained with the MARS model fitted to the spectral data corrected using the proposed algorithm (6S) is shown in Figure 9.No trends to under-or over-estimate were observed.This statistical approach also enables identifying the importance of the predictive variables, in this case the near-infrared being the most important.Hence, in accordance with the nsubsets criterion, the spectral variable Band 7 (2080-2350 nm) contributed most to the fitted MARS models for the ATCOR2, COST, FLAASH and ToA algorithms, in contrast to algorithm 6S, in which Band 5 (1550-1750 nm) was the most important variable.However, the importance of the variables in algorithm 6S, under the criteria sqr_rss and sqr_gcv, indicates that Band 7 makes the greatest contribution to the MARS model, as also occurred with the other algorithms tested.In this respect, the variable Band 7 is associated with the moisture content of the forest stand, so that the greater the moisture content of the land cover, the greater absorption in this band of the electromagnetic spectrum; in other words, the surface reflectance captured by the sensor tends to decrease [74,75].
Among the topographical variables considered, the slope (β) and the wetness index (WI) contributed most to defining the MARS models and is also directly related to the moisture content in the field sites as a factor that controls run-off; the lower the value of β, the higher the moisture content and, therefore, the higher the AGB [18,76,77].Hence, it is quite clear that radiometric and topographic variables related to water availability played a key role in the prediction of AGB in Durango's temperate forests.
Regarding the inclusion of vegetation indices, NDVI and MSAVI2 were the variables that contributed most to the MARS models, given their potential usefulness in estimating AGB.Inclusion of these indices improved the prediction of the biophysical variables in forest ecosystems: in the case of NDVI because of it high sensitivity to the chlorophyll content of the vegetation [78,79] and in the case of MSAVI2 because of its capacity to discriminate soil from vegetation, even in zones with scarce vegetation cover [80].

Conclusions
The study findings demonstrate the potential of the combined application of Landsat Thematic Mapper (TM) sensor satellite imagery and the consideration of topographical variables for estimating AGB, as well as the effects of different radiometric correction algorithms on the estimates obtained.The algorithm that showed the greatest spectral stability and that best estimated the AGB This statistical approach also enables identifying the importance of the predictive variables, in this case the near-infrared being the most important.Hence, in accordance with the nsubsets criterion, the spectral variable Band 7 (2080-2350 nm) contributed most to the fitted MARS models for the ATCOR2, COST, FLAASH and ToA algorithms, in contrast to algorithm 6S, in which Band 5 (1550-1750 nm) was the most important variable.However, the importance of the variables in algorithm 6S, under the criteria sqr_rss and sqr_gcv, indicates that Band 7 makes the greatest contribution to the MARS model, as also occurred with the other algorithms tested.In this respect, the variable Band 7 is associated with the moisture content of the forest stand, so that the greater the moisture content of the land cover, the greater absorption in this band of the electromagnetic spectrum; in other words, the surface reflectance captured by the sensor tends to decrease [74,75].
Among the topographical variables considered, the slope (β) and the wetness index (WI) contributed most to defining the MARS models and is also directly related to the moisture content in the field sites as a factor that controls run-off; the lower the value of β, the higher the moisture content and, therefore, the higher the AGB [18,76,77].Hence, it is quite clear that radiometric and topographic variables related to water availability played a key role in the prediction of AGB in Durango's temperate forests.
Regarding the inclusion of vegetation indices, NDVI and MSAVI2 were the variables that contributed most to the MARS models, given their potential usefulness in estimating AGB.Inclusion of these indices improved the prediction of the biophysical variables in forest ecosystems: in the case of NDVI because of it high sensitivity to the chlorophyll content of the vegetation [78,79] and in the case of MSAVI2 because of its capacity to discriminate soil from vegetation, even in zones with scarce vegetation cover [80].

Conclusions
The study findings demonstrate the potential of the combined application of Landsat Thematic Mapper (TM) sensor satellite imagery and the consideration of topographical variables for estimating AGB, as well as the effects of different radiometric correction algorithms on the estimates obtained.
The algorithm that showed the greatest spectral stability and that best estimated the AGB after application of the cross-validation technique was the 6S algorithm.However, each algorithm has advantages and disadvantages, and the user must parameterize each algorithm on the basis of the specific objectives.The MARS statistical method proved suitable for AGB estimation by using easy-to-obtain variables from remote sensing techniques.For all algorithms considered, the spectral band in the mid-infrared region (Band 7), the slope (β) and the wetness index (WI) of the land, both related to water availability, along with the vegetation indices MSAVI2 and NDVI were the most important predictor variables for estimating AGB.Generation of an AGB map from the fitted models may be used in local and regional forest management, thus facilitating the location and precise delimitation of zones with different levels of AGB.Further research could aim at testing the effects of radiometric corrections on biomass modelling on a multi-scene and/or multi-temporal basis.

Figure 1 .
Figure 1.Location of the study area and of the permanent forest growth and soil research sites (SPIFyS).

Figure 1 .
Figure 1.Location of the study area and of the permanent forest growth and soil research sites (SPIFyS).

Figure 2 .
Figure 2. Flow diagram of the processes involved in estimating AGB.

Figure 2 .
Figure 2. Flow diagram of the processes involved in estimating AGB.

Figure 4 .
Figure 4. Spectral performance of each of the radiometric correction algorithms considered.

Figure 4 .
Figure 4. Spectral performance of each of the radiometric correction algorithms considered.

Figure 5 .
Figure 5. Box plots for each band of the Landsat 5 TM sensor and corresponding groupings of the radiometric correction algorithms (different letters indicate significant differences between algorithm performances, at p < 0.05).

Figure 5 .
Figure 5. Box plots for each band of the Landsat 5 TM sensor and corresponding groupings of the radiometric correction algorithms (different letters indicate significant differences between algorithm performances, at p < 0.05).

Figure 6 .
Figure 6.Importance and selection of predictor variables in the multivariate adaptive regression splines (MARS) models for each radiometric correction algorithm considered.

Figure 6 .
Figure 6.Importance and selection of predictor variables in the multivariate adaptive regression splines (MARS) models for each radiometric correction algorithm considered.

Figure 7 .
Figure 7. Variations in the goodness of fit statistics obtained by applying the cross-validation technique (nfolds = 2-10) to the MARS models.

Figure 8 .
Figure 8. Maps of AGB in the study area generated from images corrected using the different radiometric correction algorithms considered.(a) ATCOR2; (b) COST; (c) FLAASH; (d) 6S; and (e) ToA.

Figure 7 .
Figure 7. Variations in the goodness of fit statistics obtained by applying the cross-validation technique (nfolds = 2-10) to the MARS models.

Figure 7 .
Figure 7. Variations in the goodness of fit statistics obtained by applying the cross-validation technique (nfolds = 2-10) to the MARS models.

Figure 8 .
Figure 8. Maps of AGB in the study area generated from images corrected using the different radiometric correction algorithms considered.(a) ATCOR2; (b) COST; (c) FLAASH; (d) 6S; and (e) ToA.

Figure 8 .
Figure 8. Maps of AGB in the study area generated from images corrected using the different radiometric correction algorithms considered.(a) ATCOR2; (b) COST; (c) FLAASH; (d) 6S; and (e) ToA.

Figure 9 .
Figure 9. Observed versus predicted AGB obtained using the MARS model fitted to the spectral data corrected using the 6S algorithm.The broken line corresponds to the diagonal.

Figure 9 .
Figure 9. Observed versus predicted AGB obtained using the MARS model fitted to the spectral data corrected using the 6S algorithm.The broken line corresponds to the diagonal.

Table 1 .
Descriptive statistics of the main stand variables estimated from the 99 permanent sample plots.The dominant stand height was calculated as an average value from the 100 thickest trees per hectare.

Table 1 .
Descriptive statistics of the main stand variables estimated from the 99 permanent sample plots.The dominant stand height was calculated as an average value from the 100 thickest trees per hectare.

Table 3 .
MARS model selection criteria for AGB estimation and the different radiometric correction algorithms considered.

Table 3 .
MARS model selection criteria for AGB estimation and the different radiometric correction algorithms considered.