Manizheh Rajab Pourrahmati

DR, IRD, Cameroun Pierre Couteron

Examiner, President of jury Mehrez Zribi

Engineer IGN, HDR Clément Mallet

HDR, MCF Pierre-Louis Frison

HDR, MCF Olivier Strauss

HDR Jean Stéphane Bailly

DR, Irstea Nicolas Baghdadi

Irstea Supervisor

Prof Ali A Darvishsefat

HDR Valéry Gond

Cirad Advisor

président Cameroun Rapporteur

Mallet, Ingénieur IGN, HDR Lastig Examinateur

Agroparistech Invité

Investigating the possibility of forest height/volume estimation using lidar, radar and optical images (Case study: Nowshahr Forests in Mazindaran, Iran

Keywords: Lidar, ICESat GLAS, Alos PALSAR, Optical images, Maximum canopy height, Lorey's height, Forest volume, Iran

Iran)

Résumé

L'importance de mesurer les paramètres biophysiques de la forêt pour la surveillance de la santé des écosystèmes et la gestion forestière encourage les chercheurs à trouver des méthodes précises et à faible coût en particulier sur les zones étendues et montagneuses.

Dans la présente étude, Le lidar satellitaire GLAS (Geoscience Laser Altimeter System) embarqué à bord du satellite ICESat (Ice Cloud and land Elevation Satellite) a été utilisé pour estimer trois caractéristiques biophysiques des forêts situées dans le nord de l'Iran: 1) hauteur maximale de la canopée (H max ), 2) hauteur de Lorey (H Lorey ), et 3) le volume du bois (V). Des régressions linéaires multiples (RLM), des modèles basés sur les Forêts Aléatoires (FA : Random Forest) et aussi des réseaux de neurones (ANN) ont été développés à l'aide de deux ensembles différents de variables incluant des métriques obtenues à partir des formes d'onde GLAS et des composantes principales (CP) produites à partir de l'analyse en composantes principales (ACP) des données GLAS. Pour valider et comparer les modèles, des critères statistiques ont été calculées sur la base d'une validation croisée. Le meilleur modèle pour l'estimation de la hauteur maximale a été obtenu avec une régression RLM (RMSE = 5.0 m) qui combine deux métriques extraites des formes d'onde GLAS (étendue et hauteur pour une énergie à 50%, respectivement W ext et H 50 ), et un paramètre issu du modèle numérique d'élévation (Indice de relief TI). L'erreur moyenne absolue en pourcentage (MAPE) sur les estimations de la hauteur maximale est de 16.4%. Pour la hauteur de Lorey, un modèle basé sur les réseaux de neurones et utilisant des CPs et le W ext fournit le meilleur résultat avec RMSE = 3.4 m et MAPE = 12.3%. Afin d'estimer le volume du bois, deux approches ont été utilisées: (1) estimation du volume à l'aide d'une relation volume-hauteur avec une hauteur estimée à partir de données GLAS et (2) estimation du volume du bois directement à partir des données GLAS en développant des régressions entre le volume in situ et les métriques GLAS. Le résultat de la première approche (RMSE=116.3 m 3 /ha) était légèrement meilleur que ceux obtenus avec la seconde approche. Par exemple, le réseau de neurones basé sur les PCs donnait un RMSE de 119.9 m 3 /ha mais avec des meilleurs résultats que l'approche basée sur la relation volume-hauteur pour les faibles (<10 m 3 /ha) et les forts (> 800 m 3 /ha) volumes. Au total, l'erreur relative sur le volume de bois est estimée à environ 26%. En général, les modèles RLM et ANN avaient des meilleures performances par rapport aux modèles de FA. En outre, la précision sur l'estimation de la hauteur à l'aide de métriques issues des formes d'onde GLAS est meilleure que celles basées sur les CPs.

Compte tenu des bons résultats obtenus avec les modèles de hauteur GLAS (hauteurs maximale et de Lorey), la production de la carte des hauteurs d'étude par une utilisation combinée de données de télédétection lidar, radar et optique (GLAS, PALSAR, SPOT-5 et Landsat-TM) et de données environnementales (pente, aspect, et altitude du terrain ainsi que la carte géologique) a été effectuée à l'intérieur de notre zone. Ainsi, des régressions RLM et FA ont été construites entre toutes les hauteurs dérivées des données GLAS, à l'intérieur de la zone d'étude, et les indices extraits des données de télédétection et des paramètres environnementaux. Les meilleurs modèles entrainés pour estimer H max (RMSE = 7.4 m et R ! " =0.52) et H Lorey (RMSE = 5.5 m et R ! " =0.59) ont été utilisées pour produire les cartes de hauteurs. La comparaison des H max de la carte obtenue avec les valeurs de H max in situ à l'endroit de 32 parcelles produit un RMSE de 5.3 m et un R 2 de 0.71. Une telle comparaison pour H Lorey conduit à un RMSE de 4.3m et un R 2 de 0.50. Une méthode de régression-krigeage a également été utilisée pour produire une carte des hauteurs en considérant la corrélation spatiale entre les hauteurs. Cette approche, testée dans le but d'améliorer la précision de la carte de la hauteur du couvert fournie par la méthode nonspatiale, a échouée due à l'hétérogénéité de la zone d'étude en termes de la structure forestière et de la topographie. Land-use change and soil degradation are major processes for the release of CO 2 to the atmosphere. The increase in greenhouse gases (GHGs) in the atmosphere is now recognized to contribute to climate change (IPCC, 2001). Although uncertainties remain regarding the causes, consequences and extent of climate change, it is believed that human activities are having an impact on the energy balance of the earth. Its influence on the climate is a major concern in the twenty-first century. This concern has led to the 1997 international agreement in Kyoto (the so-called Kyoto Protocol) made by the UNFCCC (United Nations Framework Convention on Climate Change), whereby most countries are committed to reducing their GHG emissions to the atmosphere. Furthermore, at the Paris climate conference (COP21) in December 2015, 195 countries adopted the first-ever universal, legally binding global climate deal. The agreement sets out a global action plan to put the world on track to avoid dangerous climate change by limiting global warming to well below 2°C. Forests' role in combating climate change was formally recognized in this agreement. This recognition includes formal mention of the Reducing Emissions from Deforestation and Forest Degradation (REDD+) in the agreement text as well as new commitments to increase forest financing [START_REF] Metzel | Forests in the cop21 climate change agreement: momentum on and mention of forests in the climate change text[END_REF].

List of figures

List of tables

Introduction

Issues mentioned above highlight the increasing importance of sustainable forest management. Measuring biophysical parameters of forest is of primary steps for forest ecosystem management. Tree's height has a primary and fundamental importance among all other parameters. In fact the information about vertical structure of forest specifically height is important for ecosystem health assessment, site fertility, volume, biomass and carbon cycle measurement and monitoring [START_REF] Namiranian | Measurement of Tree and Forest Biometry, 1st Edition[END_REF][START_REF] Cairns | Carbon dynamics of Mexican tropical evergreen forests: influence of forestry mitigation options and refinement of carbon-flux estimates[END_REF].

Different heights are defined and measured for a tree including total height, commercial height, trunk height, crown height etc. In forestry, total height is defined as vertical distance between base of tree and top of it [START_REF] Husch | Forest Mensuration[END_REF][START_REF] Namiranian | Measurement of Tree and Forest Biometry, 1st Edition[END_REF]. Generally, indirect methods are used to measure forest height in the field. The most common way is a tangent method which tree's height (h) is computed through measuring angles of top (α) and base (β) of the tree by a clinometer, and horizontal distance of measurer to the tree (a

): h = a[tan α -tan β].
It is not practical to measure the height of all trees in a forest stand. But the height of some trees are measured and then based on relationship between height and diameter at breast height (at 1.3 m above ground) (DBH), the height of all stand trees is estimated (Knowing that, measuring diameter is much easier and cheaper.). Depending on the subject and the objective, the maximum canopy height, mean height, mean Lorey's height or dominant height is measured or computed. The maximum canopy height in a stand is important for assessing the site quality and trees growth rate. The mean Lorey's height as representative of mean height in uneven-aged stands is also a valuable parameter for forest ecosystem management.

Forest volume, measured in cubic meters per hectare, is considered for forest quantification. Stand volume at a nominated age is related to the site quality. Volume can also be used to estimate biomass quantities (dry weight of forest) and levels of carbon sequestered in the forest. In other words, the data for the forest biomass quantities depend importantly on the ability to measure forest volumes and conversion factors. Scientific researchers study relationship between forest biomass and biodiversity [START_REF] Vance-Chalcraft | Relationship between aboveground biomass and multiple measures of biodiversity in subtropical forest of Puerto rico[END_REF][START_REF] Lasky | The relationship between tree biodiversity and biomass dynamics changes with tropical forest succession[END_REF]. Forest carbon estimates are of scientific importance to understand the quantitative role of forest carbon sequestration in earth's climate system [START_REF] Hamburg | Estimating the carbon content of Russian forests: a comparison of phytomass/volume and allometric projections[END_REF][START_REF] Ding | European forests and carbon sequestration services: an economic assessment of climate change impacts[END_REF]. Changes in forest volume can be a good proxy for changes in forest carbon [START_REF] Cheng | Eeffect of forest structural change on carbon storage in a coastal metasequoia glyptostroboides stand[END_REF]. Hence, volume may ultimately provide the most reliable estimates of deforestation and forest carbon changes [START_REF] Macauley | Forest Measurement and Monitoring: Technical Capacity and How Good Is Good Enough?[END_REF].

The most accurate method of measuring standing forest volume is to measure the DBH and the height of each tree [START_REF] Namiranian | Measurement of Tree and Forest Biometry, 1st Edition[END_REF][START_REF] Macauley | Forest Measurement and Monitoring: Technical Capacity and How Good Is Good Enough?[END_REF]. For a large area of forest, sampling methods are used along with complex equations derived from regression models to estimate the forest volumes [START_REF] Schreuder | Mensurational aspects of forest inventory[END_REF][START_REF] Namiranian | Measurement of Tree and Forest Biometry, 1st Edition[END_REF].

For very large heterogeneous forests, measuring forest on the ground, quantitatively and qualitatively, would be prohibitively expensive, time consuming and laborious. This highlights the importance of remotely sensed data in global estimation of forest biophysical parameters. Digital large-scale remote sensing data provide a less expensive option for estimation of forest biophysical parameters over a large area, while potentially also providing accurate and unbiased estimates. Developed remote sensing techniques including photogrammetry (e.g. [START_REF] Miller | An investigation of the potential of digital photogrammetry to provide measurements of forest characteristics and abiotic damage[END_REF][START_REF] Gobakken | Comparing biophysical forest characteristics estimated from photogrammetric matching of aerial images and airborne laser scanning data[END_REF], Synthetic Aperture Radar (SAR) interferometry (InSAR) and polarimetric interferometry (PolInSAR) (e.g. [START_REF] Balzter | Forest canopy height and carbon estimation at Monks wood national nature reserve, UK, using dual-wavelength SAR interferometry[END_REF][START_REF] Garestier | Forest height inversion using highresolution P-band Pol-InSAR data[END_REF], and lidar (light detection and ranging) (e.g. [START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF][START_REF] Andersen | A rigorous assessment of tree height measurements obtained using airborne lidar and conventional field methods[END_REF]Chen, 2010;[START_REF] Chen | A support vector regression approach to estimate forest biophysical parameters at the object level using airborne lidar transects and QuickBird data[END_REF][START_REF] Duncanson | Aboveground biomass estimation using spaceborne lidar in managed conifer forests in south central British Columbia[END_REF] make us able to measure three-dimension objects. For instance, [START_REF] Garestier | Forest height inversion using highresolution P-band Pol-InSAR data[END_REF] estimated forest canopy height over pine forest stands using P-band PolInSAR data with an RMSE of 2m. [START_REF] Balzter | Forest canopy height and carbon estimation at Monks wood national nature reserve, UK, using dual-wavelength SAR interferometry[END_REF] provided a map of canopy height in Monks wood national nature reserve with special heterogeneity of vegetation type and density using dual-wavelength InSAR at X-and L-band. They evaluate the result using airborne imaging lidar data. The RMSE of the lidar canopy height estimates compared to theodolite data was 2.15 m (relative error 17.6%). The RMSE of the dual-wavelength InSAR-derived canopy height model compared to lidar (light detection and ranging) was 3.49 m (relative error 28.5%).

Concerning volume/biomass estimation, methods using radar and optical data are successful in forests with low to medium levels of biomass. Passive optical sensing have shown limited sensitivity to biomass using medium to high resolution imagery when the biomass reaches intermediate levels (150-200 Mg/ha) (Ploton et al., 2011;[START_REF] Lu | Aboveground forest biomass estimation with Landsat and lidar data and uncertainty analysis of the estimates[END_REF]. This is due to inability of optical data to detect variation in biomass density after complete closure of the canopy top, which can occur from low or intermediate biomass values (depending on forest characteristics). In SAR system, the saturation threshold of radar backscatters increases by the radar wavelength. For instance, L-band SAR systems (wavelength about 25 cm) are limited to low and intermediate biomass levels, with maximum values reaching 150t/ha (e.g. [START_REF] Sandberg | L-and P-band backscatter intensity for biomass retrieval in hemiboreal forest[END_REF]Baghdadi et al., 2014;[START_REF] Attarchi | Improving the estimation of above ground biomass using dual polarimetric PALSAR and ETM+ data in the hyrcanian mountain forest (Iran)[END_REF]. It also depends on the forest characteristics. Improvements in altimetry technology especially lidar led to the most direct measurements of forest structure including height of canopy and forest volume/biomass. To this date, canopy height estimation over large areas is best achieved using lidar data.

lidar is an active remote sensing system not limited to the canopy height, basal area, leaf area index, and canopy cover. It produces and transmits short laser pulses to the surface and objects. The returned pulses are captured by a telescope. Time delay between laser transmission and reception (t) is converted to distance (D) considering the speed of light (c= 3×10 8 m/s) (D = c × t/2). Since this ability is used for detecting height of objects, the lidar system is also called laser Altimeter. Scanning laser systems may be mounted on different platform; on tripod (terrestrial lidar system), on airplane (airborne lidar system) or on satellite (spaceborne lidar system). lidar was first developed as fixed-position terrestrial instrument for investigating atmospheric composition, clouds and aerosols. These systems produce dense point data with centimeter accuracy and are often used for localized terrainmapping applications that require frequent surveys. Another type of terrestrial lidar systems uses mobile platforms (water-based and land-based). Data collected from these platforms are highly accurate and are used extensively to map railroads, roadways, airports, buildings, harbors, and shorelines. Airplanes are the most common and costeffective platforms for acquiring lidar data over vast areas. Most airborne platforms can cover about 50 km 2 per hour and collect a large amount of detailed information for applications that require high-accuracy data (NOAA Coastal Services Center, 2012).

Airborne lidar systems are able to produce centimeter accuracy high resolution digital elevation model (DEM) in a relatively small area [START_REF] Hodgson | An evaluation of lidar-and IFSAR-derived digital elevation models in leaf-on conditions with USGS Level 1 and Level 2 DEMs[END_REF][START_REF] Mount | Representing, Modeling, and Visualizing the Natural Environment[END_REF]. In spaceborne lidar systems, the lidar instrument is mounted on satellite operating in orbits of 700-800 km altitude and collect data over large area. An advantage of spaceborne lidar over airborne is providing global coverage of earth.

Many studies have been performed using airborne lidar to estimate different forest parameters like tree and forest height [START_REF] Andersen | A rigorous assessment of tree height measurements obtained using airborne lidar and conventional field methods[END_REF][START_REF] Chen | A support vector regression approach to estimate forest biophysical parameters at the object level using airborne lidar transects and QuickBird data[END_REF][START_REF] Khorrami | Potential of lidar data for estimation of individual tree height of Acer velutinum and Carpinus betulus[END_REF], volume [START_REF] Tonolli | Mapping and modeling forest tree volume using forest inventory and airborne laser scanning[END_REF][START_REF] Mohammadi | Estimation of quantitative characteristics of forest structure using a combination of airborne lidar and digital aerial photos (Case Study: broadleaf forests Shastkalate Gorgan). A dissertation of doctor of philosophy[END_REF], basal area [START_REF] Drake | Estimation of tropical forest structural characteristics using largefootprint lidar[END_REF][START_REF] Mohammadi | Estimation of quantitative characteristics of forest structure using a combination of airborne lidar and digital aerial photos (Case Study: broadleaf forests Shastkalate Gorgan). A dissertation of doctor of philosophy[END_REF], Leaf area index [START_REF] Zhao | lidar-based mapping of leaf area index and its use for validating GLOBCARBON satellite LAI product in a temperate forest of the southern USA[END_REF][START_REF] Sabol | Usage of lidar data for leaf area index estimation[END_REF], and biomass [START_REF] Gleason | Forest biomass estimation from airborne lidar data using machine learning approaches[END_REF][START_REF] Takagia | Forest biomass and volume estimation using airborne lidar in a cool-temperate forest of northern Hokkaido, Japan[END_REF].

Generally, there are many researches employing airborne lidar for studying different forest sciences which commonly achieved suitable results but they are site, instrument and species specific [START_REF] Iqbal | Evaluating the Potential of Icesat/GLAS Data to Estimate Canopy Height in the New Forest National Park, UK[END_REF]. In Iran, which is subject of the current study, only two studies were performed using airborne lidar in broadleaf mountainous forests: Mohammadi (2013) combined airborne lidar and UltraCam-D digital images to estimate standing volume, basal area and number of trees per hectare. He predicted standing volume and number of trees per hectare using non-parametric statistic method of support vector machine (SVM) with a relative RMSE of 31.4 and 35.5, respectively, and basal area using non-parametric method of random forest (RF) with a relative RMSE of 27.2. [START_REF] Khorrami | Potential of lidar data for estimation of individual tree height of Acer velutinum and Carpinus betulus[END_REF] estimated individual tree's height using airborne lidar with sampling density of 4.5 points/m 2 . They obtained R 2 of 0.96 and 0.95 and RMSE of 1.05 m and 1.48 m for Acer velutinum (with height range of 6.8-36 m) and Carpinus betulus (with height range of 11-36 m), respectively.

The critical point about airborne lidar is that it is expensive and also the capacity to collect annual data over whole countries does not currently exist. There are some governmental restrictions that prevent access to airspace of sensitive areas or of foreign countries and also physical restrictions in employing aircraft over the remote ice caps and polar regions. While satellite orbits are not subject to these restrictions. Employing spaceborne lidar for estimation of forest biophysical parameters over large extent area was investigated since ICESat (Ice, Clouds, and Land Elevation Satellite) was launched into the space in 2003.

The Geoscience Laser Altimeter System (GLAS) onboard ICESat operated for a total of 18 missions during its operational years (2003)(2004)(2005)(2006)(2007)(2008)(2009). GLAS illuminated surface or footprint has a diameter of 70 m in diameter on average, and waveforms were acquired every 170 m along the track. GLAS has been used to retrieve forest canopy height and biomass since 2005 over planted (e.g. Rosette et al., 2008a;Baghdadi et al., 2014) or natural forests including coniferous (e.g. [START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF][START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF]Chen, 2010;[START_REF] Duncanson | Estimating forest canopy height and terrain relief from GLAS waveform metrics[END_REF][START_REF] Saatchi | Benchmark map of forest carbon stocks in tropical regions across three continents[END_REF], deciduous broadleaf (e.g. [START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF][START_REF] Mitchard | Mapping tropical forest biomass with radar and spaceborne lidar in Lop´e National Park, Gabon: overcoming problems of high biomass and persistent cloud[END_REF][START_REF] Los | Vegetation height and cover fraction between 60˚S and 60˚N from ICESat GLAS data[END_REF][START_REF] Khalefa | Retrieval of Savanna vegetation canopy height from ICESat-GLAS spaceborne lidar with terrain correction[END_REF] and mixed coniferous-broadleaf forests (e.g. [START_REF] Sun | Forest vertical structure from GLAS: An evaluation using LVIS and SRTM data[END_REF][START_REF] Xing | An improved method for estimating forest canopy height using ICESat-GLAS full waveform data over sloping terrain: A case study in Changbai mountains, China[END_REF][START_REF] Los | Vegetation height and cover fraction between 60˚S and 60˚N from ICESat GLAS data[END_REF]. The most concerning point about GLAS data is waveform extent broadening over sloped area (mainly because of the large footprint size, about 70 m), and difficulties of canopy top and ground peak identification due to mixed vegetation and ground returns [START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF][START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF]Chen, 2010). Chen (2010) has illustrated terrain slope effects and also plant size and distribution on canopy height estimation. It may occur that the base of tallest tree over a sloped terrain locates above or below the ground elevation peak that is recorded by lidar as representative of the tree base. On the other hand, top of the tallest tree over a sloped terrain may be below the maximum elevation (canopy top peak) recorded by lidar because of the presence of a shorter tree over higher slope. Another possibility is that maximum elevation recorded by lidar which is supposed to be the canopy top, is greater than the tree top. The later condition happens when very short sparse stands are located over sloped terrain. Even if the terrain is simple with approximately constant slope, a non-flat terrain might cause canopy height to be overestimated or underestimated, depending on the spatial distribution of plants within footprints. This is very challenging in our study since Hyrcanian forests in the north of Iran are mostly mountainous with considerable slopes.

The terrain information can be derived from ancillary DEMs [START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF]Rosette et al., 2008;Chen, 2010;[START_REF] Xing | An improved method for estimating forest canopy height using ICESat-GLAS full waveform data over sloping terrain: A case study in Changbai mountains, China[END_REF] or from the waveform itself [START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF][START_REF] Pang | Temperate forest height estimation performance using ICESat GLAS data from different observation periods[END_REF]. [START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF] proposed multiple transformations of three waveform metrics and then used stepwise regression to develop correction factors for broadened waveform extent (distance between two signals assigned to signal start and signal end) to estimate mean canopy height. Their algorithm estimated forest canopy height with an RMSE of 5 m (R 2 = 0.83). [START_REF] Duncanson | Aboveground biomass estimation using spaceborne lidar in managed conifer forests in south central British Columbia[END_REF] improved estimation of canopy height by modeling topography directly from GLAS waveforms and consequently, inclusion of terrain relief in canopy height estimation. They developed a model to estimate maximum relief (R 2 = 0.76) and used it for classification of the maximum relief of the area sensed by GLAS. Forest canopy height model was also developed from waveform metrics for three separate relief classes: 0-7 m (R 2 = 0.83), 7-15 m (R 2 = 0.88) and >15 m (R 2 = 0.75). The moderate relief class model resulted better predictions of forest height than the low (increasing in waveform metrics variability by terrain relief) and high (mixing of vegetation and terrain signals in waveforms) relief classes.

Chen (2010) adopted Lefsky et al.'s (2005;2007) methods to retrieved maximum canopy height over mountainous areas (average slope=20˚) including two conifers sites of tall and closed canopy and one broadleaf woodland site of shorter and sparse canopy. Three regressions (edge-extent linear and non-linear models and DEM-linear model) were used to remove slope effect on GLAS waveforms. The results showed a better performance for DEM-linear model (the best result was for the broadleaf woodland site with an RMSE of 4.88 m) rather than two edge-extent models. The author stated this result is reasonable because DEM index (difference between maximum and minimum elevation of airborne lidar DEM within each GLAS footprint) was derived from airborne lidar data, which provides the most direct and precise information about terrain variability. [START_REF] Lee | Physically based vertical vegetation structure retrieval from ICESat data: validation using LVIS in White Mountain National Forest, New Hampshire, USA[END_REF] calculated height difference between real ground elevation and last GLAS waveform return based on size of footprint and terrain slope, and corrected canopy height retrieved from ICESat/GLAS data over slope terrains. They compared GLAS with airborne LVIS (Laser Vegetation Imaging Sensor) lidar heights. LVIS height was calculated as the mean of three tallest LVIS heights (slope-corrected with the same approach of GLAS) within each GLAS footprint. Slope-corrected GLAS vegetation heights matched well with top three LVIS (Laser Vegetation Imaging Sensor) mean (slope-corrected) heights (R 2 = 0.64, RMSE = 3.7 m). They found vegetation height can be overestimated by 3 m over a 15° slope without slope correction.

Concerning forest volume/biomass estimation, researchers attempted to estimate height from lidar data and then predicting volume and biomass using volume/biomass-height relationships in large extent areas (e.g. [START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF][START_REF] Saatchi | Benchmark map of forest carbon stocks in tropical regions across three continents[END_REF][START_REF] Mitchard | Mapping tropical forest biomass with radar and spaceborne lidar in Lop´e National Park, Gabon: overcoming problems of high biomass and persistent cloud[END_REF][START_REF] Healey | A sample design for globally consistent biomass estimation using lidar data from the Geoscience Laser Altimeter System (GLAS)[END_REF]Baghdadi et al., 2014;[START_REF] Asner | Mapping tropical forest carbon: Calibrating plot estimates to a simple lidar metric[END_REF]. It was also considered to retrieve forest volume/biomass directly from waveform metrics. [START_REF] Boudreau | Regional aboveground forest biomass using airborne and spaceborne lidar in Québec[END_REF], [START_REF] Duncanson | Aboveground biomass estimation using spaceborne lidar in managed conifer forests in south central British Columbia[END_REF] and [START_REF] Zhifeng | Estimating forest aboveground biomass using HJ-1 satellite CCD and ICESat GLAS waveform data[END_REF] estimated above ground biomass (AGB) using Multiple Linear Regression (MLR) between AGB and metrics extracted from GLAS waveforms. [START_REF] Fu | Estimating forest biomass with GLAS samples and MODIS imagery in northeastern China[END_REF] and [START_REF] Nelson | Estimating Siberian timber volume using MODIS and ICESat/GLAS[END_REF] practiced the same approach in their research where non-parametric technique of neural network was employed.

Although lidar is a promising technique for forest structure measurement, it does not provide wall-to-wall coverage except for small footprint lidar for a small area. Synergistic use of multiple sensors has been used for mapping forest volume/biomass accurately with remote sensing data (e.g. [START_REF] Lefsky | Global forest canopy height map from the moderate resolution imaging spectroradiometer and the geoscience laser altimeter system[END_REF][START_REF] Peterson | Mapping forest height in Alaska using GLAS, Landsat composites, and airborne lidar[END_REF][START_REF] Mitchard | Mapping tropical forest biomass with radar and spaceborne lidar in Lop´e National Park, Gabon: overcoming problems of high biomass and persistent cloud[END_REF][START_REF] Zhifeng | Estimating forest aboveground biomass using HJ-1 satellite CCD and ICESat GLAS waveform data[END_REF][START_REF] Quiñones | Above ground biomass map of Kalimantan[END_REF]. [START_REF] Peterson | Mapping forest height in Alaska using GLAS, Landsat composites, and airborne lidar[END_REF] produced a forest height map using a combination of spaceborne lidar (ICESat GLAS), airborne lidar, Landsat ETM+ images and field data for Alaska. To do this, forest height was estimated using a multiple linear regression based on waveform metrics extracted from GLAS waveforms and airborne lidar data. The estimated heights were evaluated using field measurement data. The resulting model was applied to all waveforms covering the study area. To spatially extrapolate the GLAS-based canopy height estimates and generate a continuous forest layer, a regression tree (RT) approach was used. GLAS-derived canopy height values were used as dependent variables while extracted values from Landsat composite bands, national elevation dataset DEM and derived slope and aspect, as well as existing vegetation type map were used as independent variables to build the RT models.

These models were then applied to the input geospatial layers to generate spatially continuous maps (30 m resolution) of forest height for allover Alaska. [START_REF] Mitchard | Mapping tropical forest biomass with radar and spaceborne lidar in Lop´e National Park, Gabon: overcoming problems of high biomass and persistent cloud[END_REF] produced AGB map for Gabon's Lope National Park (5000 km 2 ) using a combination of terrain-corrected L-band radar (ALOS PALSAR), spaceborne lidar (ICESat GLAS) and ground based data. They estimated Lorey's height based on GLAS waveform features, and predict AGB from AGB-Lorey's height equation developed based on in situ data. An unsupervised classification was performed on original and synthetic bands extracted from PALSAR data to provide a map of vegetation structures (40 classes) for the study area. Then the AGB values derived from GLAS footprints were averaged within each class to produce a 100 m resolution AGB map.

Complex structure of forests in the north of Iran, vertically and horizontally, even and uneven aged stands, presence of diverse broadleaf species, severe topography, etc., brought into question the capability of GLAS data to estimate the forest canopy height and volume in such complexity. However a few literatures evaluated capability of optical images, airborne lidar and radar data for estimating forest volume/biomass in Iran [START_REF] Khorrami | Investigation on the capability of Landsat7 ETM+ data for standing volume estimation of beech stands (Case study: Sangdeh forests)[END_REF][START_REF] Kalbi | Estimation of forest structural attributes using ASTER data[END_REF][START_REF] Mohammadi | Estimation of quantitative characteristics of forest structure using a combination of airborne lidar and digital aerial photos (Case Study: broadleaf forests Shastkalate Gorgan). A dissertation of doctor of philosophy[END_REF] and two studies employed the synergy of optical and radar data to provide biomass map [START_REF] Attarchi | Improving the estimation of above ground biomass using dual polarimetric PALSAR and ETM+ data in the hyrcanian mountain forest (Iran)[END_REF][START_REF] Amini | SAR and optical images for forest biomass estimation[END_REF], there was no investigation using spaceborne lidar or combination of that with other sources of remote sensing data up to now. This study aimed first to investigate capability of GLAS data for estimation of forest canopy height and volume in part of mountainous forests of Iran. To overcome slope effects, numerous parametric and non-parametric regressions were developed based on metrics derived from GLAS waveforms (user defined metrics and metrics derived using statistical method of principal component analysis (PCA)), and also terrain index extracted from 10 and 90 m DEM. Consequently, providing a wall to wall map of forest height was on the agenda.

Hence, the best GLAS height models (maximum and Lorey's heights) were applied to all GLAS shots over study area. GLAS estimated heights were then used as reference heights to develop new height models based on indices extracted from radar and optical images and also environmental data. Canopy height maps (maximum and Lorey's heights) were produced using new height models and also regression-kriging method.

In summary, this research pursues following objectives: -Estimation of maximum canopy height and Lorey's height using ICESat/GLAS data by developing parametric and non-parametric statistical methods between forest height and metrics extracted from GLAS waveforms and DEM.

-Estimation of forest volume using ICESat/GLAS data by developing: 1) volume-height relationship, and 2) parametric and non-parametric statistical methods between forest volume and metrics extracted from GLAS waveforms and DEM -Producing forest height/volume map using a combination of ICESat/GLAS, ALOS/PALSAR, optical images and environmental data (aspect, slope and geology maps).

To address above objectives following questions were raised:

-Does GLAS estimate forest height/volume in mountainous forests of Iran with a suitable accuracy?

-What is preferable statistical method (multiple linear regression, artificial neural network, random; forest) for forest height/volume estimation using GLAS in our study area? -Whether terrain index will reduce topography effects on GLAS waveform and improve the predictions? Does the resolution of DEM affect performance of the developed model? -Is there possibility of providing forest height or volume map with a suitable accuracy in such high heterogeneity of horizontal and vertical structure? -Does regression-kriging improves maps provided from regression models?

Materials and data processing 2.1. Study area

This research was performed in Nowshahr forests, a part of Hyrcanian forests in the north of Iran (Fig. 2.1), located between 36.26 to 36.68 degrees N latitudes and 51.32 to 51.94 degrees E longitudes. It contains temperate deciduous broadleaved forests extended from 100 to 2200 meters altitude above sea level with slopes ranging from flat to greater than 80%. Covering even and uneven aged stands with various species led to a diverse structure across the study site. Depending on the site, the dominant species are oriental beech (Fagus orientalis), European hornbeam (Carpinus betulus), chestnut-leaved oak (Quercus castanifolia), Persian ironwood (Parotia persica), oriental hornbeam (Carpinus orientalis), and Persian oak (Quercus macranthera). Annual mean precipitation is 1200 mm, and average maximum and minimum temperature are 6˚C and 25˚C, respectively. We focused on part of Nowshahr forest that is covered by lidar footprints.

Description of data

Field measurements

In order to estimate forest biophysical parameters (heights and volume) using GLAS, field data were collected during leaf-on seasons as lidar data acquisition dates. In total, 60 GLAS footprints were located on the ground using global positioning system (GPS; Garmin Colorado 300), 33 plots in September 2013 and 27 plots in May 2014. Since the study area includes climax community comprised of low growth deciduous species, trees have almost no or very small growth in height and DBH (diameter at breast height). On the other hand, spaceborne lidar is less sensitive to little changes in forest volume in comparison with very high spatial resolution data like airborne lidar. Therefore, 5-6 year interval between in situ measurement and lidar acquisition was not considered as an important factor influencing the analysis.

As the first step, point layer of GLAS footprints and also roads were added to topographic and slope map of the study area in ArcGIS to consider conditions and accessibility of each plot before going on the field mission. Then geographic coordinates of the plots (location of center of GLAS footprints) were entered in GPS and were next navigated on the ground. DBH of all trees (DBH > 7.5 cm) within a 70 m diameter circle (as large as lidar footprints) were measured using caliper or tape measure in case of large diameter trees. As laser energy decreases towards the margins of the footprint and, consequently, the returned waveform is most representative of the features closest to the footprint center (Harding & Carabajal, 2005;Rosette et al., 2008a;[START_REF] Duong | Processing and Application of Icesat Large Footprint Full Waveform Laser Range Data. A dissertation of doctor of philosophy[END_REF], this was taken into account through field measurements. So totally 10 dominant heights, 5 within a 36 m diameter circle and 5 in a co-center 70 m diameter circle (outer margin of smaller circle), were measured using a clinometer. The height of the nearest tree to the plot center (or tree with lowest DBH if the nearest is a measured dominant tree), was measured in order to have sufficient trees with different DBH for building height-DBH relationships which will be discussed later in subsection 2.3.1. Figure 2.2 shows a schematic of height measurement plan in a plot.

Obtaining tree height requires the use of basic trigonometry: h = a[tg α -(tg β)] where h is the tree height, d is the horizontal distance from tree, α is the angle to the top of tree and β is the angle to the base of tree (Fig. 2.3a). On steep terrain it is almost impossible to accurately determine the horizontal distance from the tree. In situations where the ground is sloped, it is needed to measure slope distance. Once measure slope angle and slope distance was measured, horizontal distance can be calculated (Fig. 2.3b). All required information including environmental parameters such as slope, aspect and elevation level and also biophysical parameters were recorded in inventory forms as presented in table 2.1 and 2.2.

Digital elevation model

Digital elevation model was provided using two sources of data. The first one, Shuttle Radar Topography Mission (SRTM) data sampled at 3 arc-second (about 90 meters).

Elevations were measured in meters referenced to the WGS84/EGM96 geoid. As all data used in a research project should have the same coordinate system, including both horizontal and vertical aspects, geoidal heights were transferred to ellipsoidal heights by adding the geoid undulations to geoidal heights (DEM 90 ) (Equation 2.1). Figure 2.5 illustrates the relationship between geoidal heights and ellipsoidal heights. The second source of data was digital 1:25000 topographic maps with counter interval of 10 meters and used to produce DEM with 10 meter resolution (DEM 10 ).

h = H + N (2.1)
Where h, N and H stands for ellipsoidal height, geoid undulation and geoidal height, respectively.

Fig. 2.5. Schematic of the relationship between geoidal heights and ellipsoidal heights

Geological map

Rock largely controls the physico-chemical properties of resulting soils which affects tree growth and forest parameters. A geological map produced by geological survey and mineral exploitation of Iran, at the scale of 1:100000, was therefore used in this study. Two sheets of geological map inside the study area were merged and converted to raster layer with different resolutions to match other sources of data. As it is seen in figure 2.6, all classes except one belong to major category of the sedimentary rocks.

GLAS/ICESat

Lidar is an active remote sensing system that produces and transmits short laser pulses to the surface and objects. The returned pulses are captured by a telescope. Time delay between laser transmission and reception (t) is converted to distance (D) considering the speed of light (c = 3×108m/s) (Equation 2.2). Since this ability is used for detecting height of objects, the LiDar system is also called laser Altimeter.

D = c × t/2 (2.2)
Scanning laser systems may be mounted on different platform; on tripod (terrestrial lidar system), on airplane (airborne lidar system) or on satellite (spaceborne lidar system). In spaceborne lidar systems, the lidar instrument is mounted on satellite operating in orbits of 700-800 km altitude and providing global coverage of earth.

ICESat is an experimental scientific satellite launched by NASA in 2003 to measure mainly ice sheet elevations and its changes over the time and also to provide measurements of cloud and aerosol height profiles, land elevation, and vegetation cover. ICESat moves 26,000 km per hour on orbits at 600 km altitude and 94 degrees inclination to the equator. The Geoscience Laser Altimeter System (GLAS) onboard ICESat consists of three lasers that operate exclusively to measure distance, a Global Positioning System (GPS) receiver, and a star-tracker attitude determination system. The laser transmits short pulses (4 ns) of infrared light (1064 nm) for measuring the elevation of surfaces and dense clouds and visible green light (532 nm) for measuring the vertical distribution of clouds and aerosols [START_REF] Aronoff | Remote Sensing for GIS Managers[END_REF][START_REF] Pflugmacher | Regional applicability of forest height and aboveground biomass models for the Geoscience Laser Altimeter System[END_REF]. Laser pulses at 40 times per second illuminate 70 meter diameter footprints on average, spaced at 170-meter intervals along Earth's surface. Separation of the tracks is 15 km at the equator and 2.5 km at 80 degrees latitude [START_REF] Abshire | Geoscience laser altimeter system (GLAS) on the ICESat mission: On-orbit measurement performance[END_REF]NSIDC, 2012). The mean horizontal geolocation accuracy of the ground footprints is less than 5 m for all ICESat missions (NSIDC, 2014). Mean vertical geolocation accuracy was also reported by NSIDC between 0 and 3.2 cm over flat surfaces (NSIDC, 2014). GLAS operated for a total of 18 missions during its operational years (2003)(2004)(2005)(2006)(2007)(2008)(2009). Figure 2.7 illustrates operation of GLAS from ICESat orbiting the Earth. A list of ICESat missions and their operational period is seen in table 2.3 (NSIDC, 2014). It should be noted footprint shape is not fully circular. Laser 3 footprints are moderately elliptical, Laser 2 very elliptical and Laser 1 very elliptical with side-lobe (NSIDC, 2012). The different footprint size and shape make it difficult to describe the surface covered by footprints when all missions' data are employed in a project. This may causes uncertainty in estimations. Within each footprint, laser reflected energy by all intercepting objects and surfaces are collected by a telescope of 1 meter diameter and results a waveform that represents a vertical profile of laser-illuminated surfaces. In the early GLAS missions, energy of each returned pulse was telemetered in 544 bins over ice sheet and land, corresponding to a height of 81.6 m (each bin corresponds to one nanosecond) [START_REF] Brenner | Derivation of range and range distributions from laser pulse waveform analysis for surface elevations, roughness, slope, and vegetation heights[END_REF]. In highly sloped area or area where feature heights exceed 81.6 m, GLAS waveform would truncate, making it impossible to derive range information. So, in later operations height extent was increased to 150 m (1000 bins) over land, using a "waveform compression scheme" (Harding & Carabajal, 2005).

GLAS collected data were processed by National Snow and Ice Data Center (NSIDC), and 15 products at three levels of corrections; L1A, L1B and L2 were provided (Table 2.4). Since this research has two main parts (Estimation of forest biophysical parameters using GLAS data and providing height/volume map using synergy of GLAS, PALSAR, optical images and environmental data), GLAS data acquired on October 2007 and October 2008, corresponding to L3I and L3K missions respectively, were used for the first part of this research (estimation of forest biophysical parameters using ICESat GLAS data). All GLAS missions over the study area were also used for the second part of thesis to provide height map using combination of remote sensing data (lidar, radar and optical images) and environmental data. GLA01 and GLA14 among 15 products produced by NSIDC, release 33 1 , were employed to drive forest height and volume. GLA01 contains intensities of transmitted and received waveforms. GLAS digitizes these intensities as counts (0-255) which were converted to volts using calibration table. It is worth to mention the ordering of the transmitted pulse is in time order. The value of the first sample (bin) is for the sample closest to the spacecraft in time, and the value of the last sample is for the sample farthest from the spacecraft in time. In reverse, the received echo is in time-reversed order. So the value of the first sample is for the sample farthest from the spacecraft in time, and vice versa (NSIDC, 2012).

1-At the moment, there are release 34 data. A new release is created when changes occur in the input data or when improvements are made to the processing algorithms.

GLA14 is a level-2 elevation product derived from level-1 products GLA05 and GLA06. This product provides information about surface elevations. It also includes the laser footprint geolocation and reflectance, as well as geodetic, instrument, and atmospheric correction flags to filter out bad data (more explanation in subsection 2.3.3.1) for range measurements (NSIDC, 2012;[START_REF] Wang | Earth science applications of ICESat/GLAS: a review[END_REF]. As the transmitted and received pulses are assumed to have a Gaussian shape [START_REF] Brenner | Derivation of range and range distributions from laser pulse waveform analysis for surface elevations, roughness, slope, and vegetation heights[END_REF], up to 6 Gaussian peaks were fitted to the waveforms. Hence, the GLA14 contains parameters for these peaks including amplitude, area and standard deviation.

ALOS/PALSAR

The word radar stands for radio detection and ranging. In general, radar systems use modulated waveforms and directive antennas to transmit electromagnetic energy into a specific volume in space to search for targets. Objects (targets) within a search volume will reflect portions of this energy (radar returns or echoes) back to the radar. These echoes are then processed by the radar receiver to extract target information such as range, velocity, angular position, and other target identifying characteristics. In this research spaceborne radar data acquired by ALOS-1 (Advanced Land Observing Satellite) was used along with other source of remotely sensed data.

ALOS-1 (Advanced Land Observing Satellite) was launched on January 24, 2006 from the Tanegashima Space Center and completed its operation on 12 May 2011 (Japan Space Systems, 2012). ALOS satellite has three remote-sensing instruments: the Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) and for digital elevation models. The Advanced Visible and Near Infrared Radiometer type 2 (AVNIR-2) for precise land coverage observation, and the Phased Array type L-band Synthetic Aperture Radar (PALSAR). PALSAR is an active microwave sensor using L-band frequency (center frequency is 1270 MHz) to achieve cloud-free and day-and-night land observation. It was improved based on Synthetic Aperature Radar (SAR) onboard the first observation satellite, JERS-1. Four operation modes are defined for PALSAR; fine resolution mode, direct downlink mode, scanSAR mode, and polarimetric mode. FB (Fine resolution Beam) mode comprises 18 selections in the off-nadir angle range between 9.9º and 50.8º, each with 4 alternative polarizations: single polarization HH or VV, and dual polarization HH+HV or VV+VH. Out of the 72 possible FB modes, two have been selected for operational use. The direct transmission (or downlink) mode is a contingency backup mode which allows the downlink of the FB mode data to local ground stations in case the highspeed DRTS (Data Relay and Test Satellite) becomes unavailable. ScanSAR is available at a single polarization only (HH or VV) and can be operated with 3, 4, or 5 sub-beams transmitted in short (14 MHz) or long bursts (28 MHz). Out of the 12 ScanSAR modes available, the sort-burst, HH polarization, 5-beam mode has been selected for operational support. It features a 350 km swath width with an incidence angle range of 18-43º. The 14 MHz polarimetric mode provides the full quad-polarization (HH+HV+VH+VV) scattering matrix with 12 alternative off-nadir angles between 9.7º and 26.2º. Polarization is changed in every pulse of transmission signal, and dual polarization signals are simultaneously received. The operation is limited in lower incident angle in order to achieve higher performances. At the nominal off-nadir angle (21.5º), the swath width is 30 km with 30 m spatial resolution under the maximum data rate condition (240 Mbit/s) [START_REF] Ito | ALOS/PALSAR characteristics and status[END_REF] Full polarimetry (multi-polarization), off nadir pointing function and other functions of PALSAR improved the accuracy of analyzing geological structure, distribution of rocks and so on, and acquired a lot of effective data for resource exploration and other purposes. At the same time, multi-polarization was effective in acquiring vegetation information, which encouraged the data use in fields such as global and regional observation of vegetation, distinguishing feature on the ground, classification of land use and other purposes [START_REF] Polychronaki | Evaluation of ALOS PALSAR imagery for burned area mapping in Greece using object-based classification[END_REF][START_REF] Attarchi | Improving the estimation of above ground biomass using dual polarimetric PALSAR and ETM+ data in the hyrcanian mountain forest (Iran)[END_REF][START_REF] Mermoz | Biomass assessment in the Cameroon savanna using ALOS PALSAR data[END_REF]. Main characteristics of PALSAR were presented in table 2.5. PALSAR mosaic data, supplied by JAXA 1 , with 25 m resolution from dates 2007, 2008, 2009 and 2010 were used in this research. It included four layers "HH" and "HV" polarization, "Local incidence angle" and "date" layer from the date of satellite launch (January 24, 2006).

1-Japan Aerospace Exploration Agency (JAXA)

Passive optical remote sensing data

Optical remote sensing makes use of visible, near infrared and short-wave infrared sensors to form images of the earth's surface by detecting the solar radiation reflected from targets on the ground. Different materials reflect and absorb differently at different wavelengths. Thus, the targets can be differentiated by their spectral reflectance signatures in the remotely sensed images. Vegetation has a unique spectral signature which enables it to be distinguished readily from other types of land cover in an optical/near-infrared image. The reflectance is low in both the blue and red regions of the spectrum, due to absorption by chlorophyll for photosynthesis. It has a peak at the green region which gives rise to the green color of vegetation. In the near infrared (NIR) region, the reflectance is much higher than that in the visible band due to the cellular structure in the leaves. Hence, vegetation can be identified by the high NIR but generally low visible reflectance.

In this research, available cloud free images of two multispectral remote sensing systems of Landsat-TM and spot5, from 2003 onwards (since GLAS data belongs to different mission periods and do not reflect a static moment in time), were selected and analyzed.

The Landsat Thematic Mapper (TM) sensor was carried onboard Landsats 4 and 5 from July 1982 to May 2012 with a 16-day repeat cycle. It possesses seven spectral bands with a spatial resolution of 30 meters for bands 1 to 5 and 7. Spatial resolution for Band 6 (thermal infrared) is 120 meters, but is resampled to 30-meter pixels. For this study, four dates of cloud free orthorectified images of Landsat5-TM (Table 2.6), held in the USGS archives, were downloaded from EarthExplorer: http://earthexplorer.usgs.gov. SPOT-5 is the fifth satellite in the SPOT series of CNES (Space Agency of France) which was launched on May 2002 and its commercial mission was ended on 27 March 2015. On 2 April 2015, an experimental phase for the mission, SPOT-5 (Take 5), involved the satellite being tasked in a new orbit (lowered by 3 km with 5 day repeat cycle). From that date, and for 5 months until the 15th of September, SPOT5 (Take 5) observed 150 sites every five days with constant observation angles. The data were processed and distributed at CNES. Cloud free orthorectified multispectral bands of SPOT 5-HRG acquired on April and June 2015 were downloaded from the website of Centre d'Etudes Spatiales de la BIOsphère (CESBIO) (http://www.cesbio.ups-tlse.fr) and employed in this research (Table 2.6). Figure 2.8 shows a color composite of SPOT data over study area.

Data processing and information extraction

Analysis of in situ data

In situ measurements were carried out in two phases. Phase one includes 60 plots for developing GLAS height and volume models and their validation. It was intended to collect data in all elevation and slope classes, but it was not practical as a reason of lack of GLAS data or lack of forest cover in some elevation range, and also inaccessibility to some area. Phase two includes 32 plots to validate height maps produced from combination of remotely sensed and environmental data. Table 2.7 shows their conditions in terms of elevation and slope. [START_REF] Burnham | Information and likelihood theory: a basis for model selection and inference, in Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach[END_REF]. Besides statistical criteria, biological behavior of models was considered to select the best model. Table 2.9 presents six best non-linear height-DBH models and their statistical performances developed for six groups of species mentioned above. Figure 2.9 shows the height curve depicted for the six height growth models. Next, the Lorey's height was calculated using equation 2.3. Lorey's height as a mean height of a stand weights the contribution of trees to the stand height by their basal area. Therefore, it is more stable than arithmetic height specifically in uneven-aged stands.

U VWXYZ = [ \] ^×_ âb [ \] âb = [ c\_ ^d×_ âb [ c\_ ^d àb (2.3)
Where U VWXYZ , ef g , heU g and H i are Lorey's height (m), basal area (cm 2 ), diameter at breast height (cm), and height (m) of tree i, respectively, and n is total number of trees in each plot.

Volume is usually expressed quantitatively as a function of DBH and height [START_REF] Macauley | Forest Measurement and Monitoring: Technical Capacity and How Good Is Good Enough?[END_REF][START_REF] Namiranian | Measurement of Tree and Forest Biometry, 1st Edition[END_REF]. So the selected height-DBH relationships were next used to estimate the height of all trees. Local species level volume equations based on DBH and height developed by FRWO (Table 2.10) were also used to calculate per tree stem volume.

Extraction of terrain index, slope, aspect and elevation class maps from DEM

Terrain Index (TI) was calculated at the location of lidar footprints using: 1) a fine resolution DEM 10 produced based on 1:25000 topographic maps (called TI 10 ), and 2) SRTM DEM with 90 meter resolution (called TI 90 ). The elevation range within a 7×7 neighborhood of 10m-DEM (Rosette et al., 2008a;Chen, 2010b) and 3×3 neighborhood of 90m SRTM DEM (Baghdadi et al., 2014) at location of each GLAS footprint was considered as TI. The effect of using higher resolution DEM on model performance was also investigated.

Terrain index, slope in degrees, aspect in nine classes (flat, north, north-east, east, southeast, south, south-west, west and north-west) and elevation classification map in six classes ( < 100, 100-300, 300-700, 700-1200, 1200-1600 and 1600 < ) were also produced using DEM 10 for the entire study area. Figure 2.12 shows DEM 10 and four extracted maps from it. The GLA01 (Global Altimetry data) and GLA14 (Global Land Surface Altimetry data) data products were converted from binary to the ASCII format using IDLreadGLAS provided by NSIDC. Required information such as latitude, longitude, elevation, centroid elevation, and fitted Gaussian peaks were extracted from GLA14 data, and row waveforms were extracted from GLA01 data. The most important information derived from these products was summarized in table 2.12. ICESat/GLAS elevations are referenced to the TOPEX/Poseidon ellipsoid which is 70 cm smaller than WGS84 ellipsoid. For comparison between ICESat, SRTM DEM and in situ data, datasets needed to be available in the same coordinate system. So, GLAS elevations were transformed to the WGS84 ellipsoid by adding 70 cm.

i_lat GLA14

The geodetic latitude of the forty laser spots in the 1 second time frame in micro degree

i_lon GLA14

The longitude of the forty laser spots in the 1 second time frame in micro degree 6 i_elev GLA14

Surface elevation with respect to the ellipsoid at the spot location determined in mm 7 i_campaign GLA14 The campaign. i.e.: for campaign L3K, it will be "3K".

8 i_Gamp GLA14 Amplitude of each Gaussian solved for (up to six) in 0.01 volts 9 i_Garea GLA14

Area under each of the Gaussians solved for (up to six) in 0.01 volts × ns 10 i_Gsigma GLA14 Width (sigma) of each Gaussian solved for (up to six) in 0.001ns

i_satNdx GLA14

The count of the number of gates in a waveform which have an amplitude greater than or equal to saturation index threshold 12 i_FRir_qaFlag GLA14

Indicates the presence of clouds (0-15); 15 = no cloud, 14=likely presence of low clouds, etc.

13 i_4nsBgMean GLA01 Background Noise Mean Value in 0.01 counts 14 i_4nsBgSDEV GLA01 The standard deviation of the background noise in 0.01 counts Some pre-processes were applied to remove inappropriate and useless waveforms (Chen, 2010b;[START_REF] Hilbert | Influence of surface topography on ICESat/GLAS forest height estimation and waveform shape[END_REF]Baghdadi et al., 2014):

1-Eliminating waveforms affected by cloud; Flag i_FRir_qaFlag in GLA14 data indicates the estimated atmospheric conditions over each GLAS footprint using a cloud detection algorithm. To do this, waveforms with i_FRir_qaFlag=15 were kept and the rest were removed (Chen, 2010b;[START_REF] Duncanson | Estimating forest canopy height and terrain relief from GLAS waveform metrics[END_REF]. 2-Eliminating saturated waveforms; i_satNdx in GLA14 presents the count of the number of gates in a waveform which have an amplitude greater than or equal to saturation index threshold (i_satNdxTh). So only waveforms with i_satNdx=0 were used for analysis in this study (Chen, 2010b;[START_REF] Hilbert | Influence of surface topography on ICESat/GLAS forest height estimation and waveform shape[END_REF].

3-Removing noisy waveforms with a signal to noise ratio (SNR) lower than 15 (Baghdadi et al., 2014); To calculate SNR, maximum energy of samples from GLA01 (i_maxRecAmp) was divided to standard deviation of the background noise (i_4nsBgSDEV) recorded in GLA14 data. 4-Removing waveforms in which difference between centroid elevation (i_elev from GLA14) and corresponding SRTM DEM is greater than 100 meters (|ICESat -SRTM)001>| (Baghdadi et al., 2014).

Waveform metrics extraction

GLAS provides a full waveform of illuminated objects on the surface. Each waveform is a function of vertical structure of illuminated surface and their reflection properties within the footprint. Figure 2.13 illustrates a schematic of transmitted pulses over a vegetated area and a returned waveform. As mentioned before, 544 or 1000 bins (1 bin= 1 ns=15 cm) is recorded for each received waveform. Signal start and end are defined as first and last bins in the waveform where the waveform intensity exceeds background noise threshold, nσ+μ, where σ and μ recorded in GLA01 product are standard deviation and mean background noise respectively, and n=0.5,1,…,5. Different thresholds including 3σ+μ [START_REF] Sun | Forest vertical structure from GLAS: An evaluation using LVIS and SRTM data[END_REF], 4σ+μ [START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF], 4.5σ+μ (Baghdadi et al., 2014;[START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF][START_REF] Lee | Physically based vertical vegetation structure retrieval from ICESat data: validation using LVIS in White Mountain National Forest, New Hampshire, USA[END_REF] were applied in previous studies. Chen (2010b) tested different thresholds for signal start and end for each three sites from 2.5σ+μ to 5σ+μ. He found that the value of n for optimal threshold is 3.5 for signal start and 5 for signal end. Therefore, an average value of 4.5 was used in this study. [START_REF] Hilbert | Influence of surface topography on ICESat/GLAS forest height estimation and waveform shape[END_REF] stated that the optimal thresholds might differ according to the waveform types, laser periods or footprint structure. In this research the threshold was set to 4.5σ+μ, the optimum threshold used in the most studies.

The vertical distance between signal start and signal end of a waveform was computed as waveform extent (W ext ) which could be affected by terrain slope, canopy height and canopy density [START_REF] Yang | Assessment of the impacts of surface topography, offnadir pointing and vegetation structure on vegetation lidar waveforms using an extended geometric optical and radiative transfer model[END_REF].

To identify ground peak as an important feature of waveforms to extract maximum height over flat area (defined as vertical difference between signal start and ground peak) or derive some metrics in waveforms, Gaussian components (up to 6 peaks) are fitted to the row waveform [START_REF] Brenner | Derivation of range and range distributions from laser pulse waveform analysis for surface elevations, roughness, slope, and vegetation heights[END_REF]. [START_REF] Duong | ICESat full-waveform altimetry compared to airborne laser scanning altimetry over the Netherlands[END_REF], [START_REF] Sun | Forest vertical structure from GLAS: An evaluation using LVIS and SRTM data[END_REF] and [START_REF] Xing | An improved method for estimating forest canopy height using ICESat-GLAS full waveform data over sloping terrain: A case study in Changbai mountains, China[END_REF] have considered the last peak as ground peak. But in dense vegetated area, ground peak may have lower amplitude than adjacent peaks. In this case, it was suggested to consider the peak with highest amplitude as the ground peak even if it is not the last peak [START_REF] Boudreau | Regional aboveground forest biomass using airborne and spaceborne lidar in Québec[END_REF]. Rosette el al. (2008b) considered the centroid of the Gaussian peak with greatest amplitude out of the last two peaks as ground surface. [START_REF] Hilbert | Influence of surface topography on ICESat/GLAS forest height estimation and waveform shape[END_REF] applied two approaches to find the ground return; the first was the maximum of the last two Gaussian peaks of a waveform, and the second approach was identifying ground return based on local maxima within the GLA01 waveform. They found both approaches represented the terrain accurately and could be used as basis to calculate the tree height. However the usage of the original waveform (GLA01) is more flexible and meets the ground return more precisely. Chen (2010b) found the stronger peak of the two last one is a better representative of the ground elevation in the studied coniferous site, whereas for the studied woodland site the strongest peak of the five last peaks matched well with the ground elevation. [START_REF] Iqbal | Evaluating the Potential of Icesat/GLAS Data to Estimate Canopy Height in the New Forest National Park, UK[END_REF] detected the peak with maximum amplitude between signal start and end of a waveform. Then by moving from this point to the signal end, the peak with amplitude one-fifth the maximum amplitude was considered as ground return. If this condition was not met, the peak with maximum amplitude was assigned as ground return. In the current research, we followed last achievements (Baghdadi et al., 2014;[START_REF] Fayad | Canopy height estimation in French Guiana with lidar ICESat/GLAS data using principal component analysis and random forest regressions[END_REF]Rosette et al., 2008a;Chen, 2010b) and we chose the stronger one among two last Gaussian peaks as ground peak. The first Gaussian peak was selected as canopy top. The distance between ground peak and signal start has been defined as maximum canopy height in flat area. The vertical distance from ground peak to signal end and from canopy top to signal start was considered as trail edge and lead edge extents, respectively (Baghdadi et al., 2014;[START_REF] Yang | Assessment of the impacts of surface topography, offnadir pointing and vegetation structure on vegetation lidar waveforms using an extended geometric optical and radiative transfer model[END_REF]. H 25 , H 50 , H 75 and H 100 as quartile heights have been extracted from waveforms by calculating the vertical distance between ground peak and position of waveform at which respectively 25%, 50%, 75% and 100% of the returned energy between signal start and end occurs [START_REF] Nelson | Estimating Siberian timber volume using MODIS and ICESat/GLAS[END_REF][START_REF] Sun | Forest vertical structure from GLAS: An evaluation using LVIS and SRTM data[END_REF]. So the total waveform energy was calculated by summing all the return energies from signal start to end. Starting from the signal end, the position of the 25%, 50%, and 75% of energy were located by comparing the accumulated energy with total energy. H 100 Along track laser pulses is the maximum canopy height as defined above. Figure 2.14 illustrates a GLAS waveform from study area with Gaussian peaks and some extracted metrics. The metrics extracted from GLAS waveforms and their derivatives, used in this research, are listed in table 2.13. Over mountainous areas with large relief and complex terrain, the peaks from ground and surface objects can be broadened and signals returned from ground and vegetation be mixed. It makes that difficult to identify the ground elevation, and subsequently calculation of metrics which are dependent on the ground location (Chen, 2010b;[START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF][START_REF] Zwally | ICESat's laser measurements of polar ice, atmosphere, ocean, and land[END_REF][START_REF] Lee | Physically based vertical vegetation structure retrieval from ICESat data: validation using LVIS in White Mountain National Forest, New Hampshire, USA[END_REF]. To clarify, the two most commonly used height metrics to infer vegetation height from lidar are waveform extent (vertical distance between signal start and end) and H 100 (vertical distance between signal start and ground peak) (Fig. 2.15a). These metrics would be influenced by different factors such as surface topography, footprint size, forest density and laser pulse energy.

Over flat area H 100 equals to the canopy height (H), but over sloped terrain it may increase to "H + d×tanθ" [START_REF] Lee | Physically based vertical vegetation structure retrieval from ICESat data: validation using LVIS in White Mountain National Forest, New Hampshire, USA[END_REF] where d is lidar footprint diameter and θ is the slope angle (Fig. 2.15b). Broadening effect of terrain slope on waveform extent and also Gaussian peaks have been illustrated also in figure 2.15 (b,c,d). In addition to deterministic heuristics (user defined metrics), a non-parametric statistical technique named principal component analysis (PCA) was used to remove noises and reduce dimensionality of waveform signals.

In short, PCA finds a set of synthetic variables (the principal components) that summarizes the original set. It rotates the axis of variation to give a new set of ordered orthogonal axis that summarizes describing proportions of the variations. In fact, the principal components (PCs) are uncorrelated and ordered such that the k th PC has the k th largest variance among all PCs [START_REF] Ulfarsson | Model Based Principal Component Analysis with Application to Functional Magnetic Resonance Imaging, A dissertation of doctor of philosophy[END_REF]. The traditional approach is to use the first few PCs in data analysis since they have most of the variation in the original data set.

In this study, lidar signal intensities were used for the PCA analysis. In order to apply PCA, it is necessary to have equal number of samples in all waveforms. So, the length of largest waveform extent was considered as basis (400 samples) and other waveforms were apart from signal start toward signal end till the number of samples reach the base W ext 60) is less than the number of samples in the useful part of waveforms (400 samples), it was aimed to reduce the number of samples by selecting one among each ten samples. So PCA was performed using 41 samples as variables using package "FactoMineR" in R [START_REF] Husson | Package 'FactoMineR': Multivariate exploratory data analysis and data mining[END_REF] to find the main factors (waveform signals) determining most effects on forest canopy height. This procedure of data reduction was also used by [START_REF] Fayad | Canopy height estimation in French Guiana with lidar ICESat/GLAS data using principal component analysis and random forest regressions[END_REF] to estimate canopy height using ICESat GLAS data. As it is seen in figure 2.16, three first components had the most information, and explained 77.5% of variance in the data.

Fig. 2.16. Information explained by ten first PCs and other PCs

In overall, two sets of metrics were extracted from waveform; one, deterministic heuristics representing vertical distance between different positions of waveform and ground peak (represented in Table 2.13), and one non-parametric metrics including principal components produced from PCA. In this dissertation, these metrics are addressed as "waveform metrics" and "PCs", respectively, however PCs are also metrics derived from waveforms.

PALSAR data processing and extraction of metrics

PALSAR mosaic data which were used in this study have been processed by JAXA using a long mosaicking algorithm described by [START_REF] Shimada | Generating large-scale high-quality SAR mosaic datasets: application to PALSAR data for global monitoring[END_REF]. It includes orthorectification, slope correction and intensity tuning between neighboring strips. Absolute radiometric calibration were performed to drive sigma naught and gamma naught (backscatter coefficient) from each polarization of PALSAR mosaic data (HH and HV), using equations 2.4 and 2.5 [START_REF] Shimada | Generating large-scale high-quality SAR mosaic datasets: application to PALSAR data for global monitoring[END_REF][START_REF] Mermoz | Biomass assessment in the Cameroon savanna using ALOS PALSAR data[END_REF]. Equations 2.6 and 2.7 were also used to convert them to the linear scale. i Q@ = $10$ × jkl +m 9In " C 4 $op! ! (2.4) q Q@ = $ i Q@ 4 10$ × jkl +m 9rst uC (2.5)

i v = $ 10 9w xy +m L C
(2.6)

q v = $ i v × rst u
(2.7)

Where i Q@ and q Q@ stands for gamma and sigma naught in dB, respectively. i j and q v stands for gamma and sigma in linear scale, respectively. The θ is local incidence angle, DN stands for digital number which is pixel intensity value and CF is constant calibration factor equal to -83 dB.

Speckle noise as defined by [START_REF] Gagnon | Speckle Filtering of SAR Images -A comparative study between complex-wavelet-based and standard filters[END_REF], is a common phenomenon in all coherent imaging systems like laser and SAR imagery. The source of this noise is attributed to random interference between the coherent returns issued from the numerous scatters present on a surface, on the scale of a cell resolution. Speckle noise is often an undesirable effect, and so, speckle filtering turns out to be a critical pre-processing step for detection/classification optimization. There are many speckle reduction techniques. In this study "LEE" as a widely used filter and then "multitemporal" filter (four dates, 2007 to 2010, at HH and HV polarizations) were applied on the data to reduce speckle effects [START_REF] Attarchi | Improving the estimation of above ground biomass using dual polarimetric PALSAR and ETM+ data in the hyrcanian mountain forest (Iran)[END_REF][START_REF] Mermoz | Biomass assessment in the Cameroon savanna using ALOS PALSAR data[END_REF]. A 5 by 5 window was used for performing LEE filter on q v and i j in "nest" software. Then multitemporal filter was performed on four images from 2007-2010. Quantitatively the distribution of the backscatter of a homogenous area is reduced significantly after these filters. An illustration of speckle reduction is shown in figure 2.17a,b. The histogram of a forested area before and after filter is seen in figure 2.17c. The standard deviation of the backscatter coefficient (here gamma naught which is shown by "γ˚" hereinafter) is decreased from 1.42 to 0.44 decibel. Image texture which is defined as variation of image tones that are related to the spatial distribution of forest vegetation [START_REF] Roberts | Forest structural assessment using remote sensing technologies: an overview of the current state of the art[END_REF] has proved to be capable of identifying different aspects of forest stand structure. [START_REF] Kayitakire | Retrieving forest structure variables based on image texture analysis and IKONOS-2 imagery[END_REF] showed that textural indices derived from a Grey-level co-occurrence matrices (GLCM) of an IKONOS-2 image are well correlated to the forest structural variables such as age, crown circumference, tree height, stand density and basal area. [START_REF] Nichol | Improved biomass estimation using the texture parameters of two high-resolution optical sensors[END_REF] and [START_REF] Attarchi | Improving the estimation of above ground biomass using dual polarimetric PALSAR and ETM+ data in the hyrcanian mountain forest (Iran)[END_REF] resulted in heterogeneous forests, texture measures are more sensitive to the canopy structure than spectral reflectance, and are more correlated to forest AGB. [START_REF] Trinder | Relating worldview-2 data to pine plantation lidar metrics[END_REF] also indicated that the textural-based models are significantly more efficient than spectral-based models for predicting lidar metrics. The texture measures are able to identify significant differences in image texture independently of image contrast (i.e. backscatter). It also increases the saturation threshold and the biomass range that can be measured (Kuplich et al., 2005;[START_REF] Sarker | Potential of texture measurements of two-date dual polarization PALSAR data for the improvement of forest biomass estimation[END_REF][START_REF] Cutler | Estimating tropical forest biomass with a combination of SAR image texture and Landsat TM data: An assessment of predictions between regions[END_REF]. However texture characterization is influenced by window size, it is not easy to select a superior window because an optimal window depends on different characteristics of forest and GLCM attributes are affected by the window size in different ways [START_REF] Ouma | Optimization of second-order grey-level texture in highresolution imagery for statistical estimation of above-ground biomass[END_REF][START_REF] Trinder | Relating worldview-2 data to pine plantation lidar metrics[END_REF]. Consequently, the mean value of each feature from multitemporal data (from 2007 to 2010) was calculated and used. In order to match PALSAR data, optical images and DEM spatially, the produced maps or indices were resampled to 10 m resolution. Figure 2.18 shows mean values of multitemporal data for eight statistical texture features extracted from HH band using GLCM algorithm. As it is seen, variation of "contrast" and "variance" over study area is low. There is no or very low correlation between canopy height/volume and most extracted features at the location of in situ plots. The "mean" and "correlation" showed highest correlation with canopy height/volume. Despite of low correlation of these features with considered parameters, individually, combination of these measures with other variables (extracted from other sources of data) will be considered in the height estimation. (2.9)

R{| = $ z|R R~• L (2.10)
Several other maps were also produced using multitemporal data of vegetation indices after resampling TM derived indices to 10 m resolution: 1) minimum, maximum and mean values of NDVIs, MVIs and RVIs (So-called min-ndvi, max-ndvi, mean-ndvi, etc.), and also mean-summer and mean-winter values of each index; 2) eight statistical features including "mean", "variance", "homogeneity", "contrast", "dissimilarity", "second moment", "entropy" and "correlation" derived from GLCM texture analysis on mean NDVI map using 3×3, 5×5, 7×7, 9×9 and 11×11 window size. Figures 2.16 shows mean NDVI and eight texture features derived from it using GLCM algorithm, respectively. In contrast to the texture features extracted from PALSAR data, almost all texture features from NDVI have moderate to good correlation with forest canopy height/volume. Correlation between "mean" and Lorey's height, for instance, was 0.68. The "correlation" is the only feature not correlated with H max and H Lorey . As it is seen in the figure 2.19, the variation of "correlation" values over forested area is very low.

. Direct method for estimating H max

Over flat area, estimation of maximum canopy height (H max ) is based on vertical difference between the waveform signal start (Ss) and the ground peak (Gp) [START_REF] Neuenschwander | Characterization of ICESat/GLAS waveforms over terrestrial ecosystems: Implications for vegetation mapping[END_REF]Chen, 2010b), and it is calculated using equation 3.1. Vertical resolution of waveforms is 15 cm for GLAS data (Harding & Carabajal, 2005).

U €!• 9'ƒ$"C = $ 9… †$ : $ ‡tC × 0.15 (3.1)
As it was described in section 2.3.3.2, over non-flat terrain the width of waveform and Gaussian peaks increases by effect of surface roughness. So the distance between signal start and ground peak will not represent the canopy height.

Parametric and non-parametric methods for prediction H max and H Lorey

Over sloped terrain, peaks from ground and surface objects can be broadened and mixed, making identification of ground peak difficult [START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF][START_REF] Pang | Temperate forest height estimation performance using ICESat GLAS data from different observation periods[END_REF]Chen, 2010b). Hence it is necessary to find a way to decrease slope impact on waveform. [START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF] and Chen (2010b) used DEM to include topography effects on height estimations. Some researchers derived terrain information from the waveform itself based on metrics such as leading and trailing edge extent (H lead and H trail , respectively). [START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF][START_REF] Pang | Temperate forest height estimation performance using ICESat GLAS data from different observation periods[END_REF] used multiple transforms of the leading and trailing edge extents to model the H lead and H trail correction factor. In the next step, height correction factor (cf), defined as difference between waveform extent (W ext ) and mean tree height (H mean ), was regressed against H lead and H trail correction factor. Then it was used to estimate mean tree height (H mean = W ext -cf). In present research, Terrain Index (TI) extracted from two sources of DEM (section 2.3.2) was entered as a predictor to consider the topography condition in the models. The effect of using finer resolution DEM (DEM 10 in replace of DEM 90 ) on performance of models was also investigated.

A large number of multiple linear regression (MLR), Random Forest (RF) and Artificial Neural Network (ANN) models were developed in R, employing different combination of metrics extracted from waveforms ("waveform metrics" and "PCs") and also TI.

Random forests, as an ensemble learning method developed by Breiman in 2001, operate by constructing a multitude of regression trees [START_REF] Breiman | Bagging predictors[END_REF]. Each tree in the forest is made of a random subset of observations with replacement and also a random of explanatory variables. So, two important parameters in random forest are the number of trees in the forest and the number of variables in the random subset at each node of tree. Prediction of new set of data for regression application would be the average prediction of all trees [START_REF] Breiman | Random Forests[END_REF][START_REF] Liaw | Classification and regression by randomforest[END_REF]. Therefore, the random forests algorithm (for both classification and regression) contains following steps:

1. Draw n tree bootstrap samples from the original data.

2. For each of the bootstrap samples, grow an unpruned classification or regression tree. So that at each node, the best split is chosen among randomly selected variables.

3. Predict new data by aggregating the predictions of the n tree trees (i.e., majority votes for classification, average for regression).

Fig. 3.3. A schematic of random forest regression

Through random sampling of observations, about one-third of them are not used for any individual tree that is called out of the bag, "OOB", for that tree. The accuracy of a random forest's prediction can be estimated from these OOB data [START_REF] Breiman | Random Forests[END_REF][START_REF] Grömping | Variable importance assessment in regression: linear regression versus random forest[END_REF]. It would be possible to calculate variable importance by determining how much worse would be the OOB predictions, if the data for that variable are randomly permuted [START_REF] Liaw | Classification and regression by randomforest[END_REF][START_REF] Wei | Global patterns and predictions of seafloor biomass using random forests[END_REF]. In fact, it would be possible to find out what would happen with or without the help of that variable. Variable importance measures produced by RF can also sometimes be useful to build simpler model. One warning about RF is that they are dependent on the training set. If the training set is not actually representing the population, it is supposed to obtain inaccurate predictions for values out of training set [START_REF] Horning | Random Forests: An algorithm for image classification and generation of continuous fields data sets[END_REF]. RF regressions were developed using different combination of predictors ("waveform metrics" and "PCs") utilizing "randomForest" package by Liaw & Wiener, 2014 in R. Artificial neural networks are able to perform nonlinear modeling without a prior knowledge about the relationships between input and output variables. It is also a nonparametric and black-box model. Thus they are a more general and flexible modeling tool for forecasting.

A variety of neural network structures have been developed. In this research a Multilayer perceptron (MLP) as the most popular type of neural networks belong to a general class of structures called "feedforward" have been used. An MLP is composed of several layers of neurons (nodes); the first layer "input layer" for distributing the data into network, the last one "output layer" for extraction the result of the network, and remaining layers between input and output are called hidden layers. There are complete connections between neurons in successive layers. Each neuron, except input layer neurons, is obtained by computing weighted sum of previous layer neurons and applying an activation function. Figure 3.4 shows a MLP neural network consisting one hidden layer.

It utilizes a learning technique called backpropagation for training the network. This kind of ANN is based on supervised learning. The idea of backpropagation algorithm is that output of neural network is evaluated against desired output. If results are not satisfactory, weights between layers are modified and process is repeated until error is small enough.

The answer that emerges from a neural network's weights can be difficult to understand and the network's training can take longer than certain other methods of machine learning such as random forests.

Fig. 3.4. A schematic of MLP neural network

Generally, an MLP is characterized by the number of hidden layers, hidden neurons, output neurons and transfer functions. An MLP combined of one input layer, one output layer and one or more hidden layers. In theory a hidden layer with sufficient number of hidden neurons is capable of approximating any continuous function [START_REF] Kaastra | Designing a neural network for forecasting financial and economic time series[END_REF]; output [START_REF] Zhang | Forecasting with artificial neural networks: The state of the art[END_REF]. The number of neurons in input layer equals to the number of variables, and in output layer depends on the application (usually one output neuron for regressions).

There is no formula for setting an optimum number of hidden neurons. [START_REF] Katz | Developing neural network forecasters for trading[END_REF] indicated that optimum number of hidden neurons is between one-third and two or three times the number of input neurons. [START_REF] Bailey | Developing neural network applications[END_REF] suggested that the number of hidden neurons for a three layer ANN should be 75% of the number of input neurons.

The relationship between input and output of a neuron or network is determined by an activation function. There are a number of common activation functions in use with ANN including "Linear", "Gaussian", "Sigmoid", "Hyperbolic tangent" etc. Selection of activation function is arbitrary and is usually determined by response variable [START_REF] Günther | Neuralnet: training of neural networks[END_REF]. In this research the hyperbolic tangent and linear function were used as activation function, respectively in hidden and output layers. A numerous ANN models were developed using different combination of input variables ("waveform metrics" and "PCs"), and different number of hidden layers and neurons. The models were developed in R using "monmlp" package [START_REF] Cannon | Package 'monmlp': Monotone multi-layer perceptron neural network[END_REF].

As mentioned, the idea of using terrain index and edge extents came to remove the broadening effects of sloped terrain. It was questioned if other waveform metrics could improve the result. The selection of predictors (metrics) was initially based on the experience and knowledge about impact of extracted metrics on the forest height (regarding to the literatures). A stepwise regression was also used to select the best combination of predictors. It combines backward elimination and forward selection to reach the best combination of metrics based on AIC criteria. The selected metrics were employed in MLR, RF and ANN models. However numerous models were built using waveform metrics and PC S , a few MLR models developed using waveform metrics were presented in table 3.1 as instance.

Forest volume estimation using ICESat GLAS

Two methods were applied to estimate forest volume. The first method consists of three steps: 1) developing volume-H max and volume-H Lorey relationships. The stronger one was chosen to estimate volume next. The common volume-height relationship (Equation 3.2) used in different literatures (Baghdadi et al., 2014;[START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF][START_REF] Saatchi | Benchmark map of forest carbon stocks in tropical regions across three continents[END_REF][START_REF] Mitchard | Mapping tropical forest biomass with radar and spaceborne lidar in Lop´e National Park, Gabon: overcoming problems of high biomass and persistent cloud[END_REF][START_REF] Healey | A sample design for globally consistent biomass estimation using lidar data from the Geoscience Laser Altimeter System (GLAS)[END_REF], was calibrated based on collected in situ data; 2) estimating height from GLAS data using best model resulted from subsection 3.1.2. It should be mentioned that if volume-H Lorey relationship is chosen at the first step, Lorey's height would be estimated form lidar data; and 3) estimating forest volume (V) using chosen volume-height relationship. This method has been used in several studies (Baghdadi et al., 2014;[START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF][START_REF] Mitchard | Mapping tropical forest biomass with radar and spaceborne lidar in Lop´e National Park, Gabon: overcoming problems of high biomass and persistent cloud[END_REF][START_REF] Healey | A sample design for globally consistent biomass estimation using lidar data from the Geoscience Laser Altimeter System (GLAS)[END_REF]:

{ = ˆ. U ‰ (3.2)
Where V is volume in m 3 /ha and H is H max or H Lorey in m.

The second method estimates forest volume directly from GLAS waveforms [START_REF] Duncanson | Aboveground biomass estimation using spaceborne lidar in managed conifer forests in south central British Columbia[END_REF][START_REF] Zhifeng | Estimating forest aboveground biomass using HJ-1 satellite CCD and ICESat GLAS waveform data[END_REF][START_REF] Nelson | Estimating Siberian timber volume using MODIS and ICESat/GLAS[END_REF][START_REF] Hayashi | Regional forest biomass estimation using ICESat/GLAS spaceborne lidar over Borneo[END_REF]. In fact, a large number of MLR, RF and ANN regressions were developed based on waveform metrics, and PCs to predict forest volume.

Model Validation

A Cross validation allows models to be tested using the full training set by means of repeated resampling; thus, maximizing the total number of points used for testing and potentially, helping to protect against overfitting [START_REF] Rao | On the dangers of cross-validation, an experimental evaluation[END_REF]. In a k-fold cross validation, the dataset is randomly split into k subsets, and each fold uses one of the k subsets as test and the remaining data for training. This process is then repeated k times (the folds), with each of the k subsamples used exactly once as the validation data.

The k results from the folds can then be averaged (or otherwise combined) to produce a single estimation (Fig. 3.5). In this study, based on a 5-fold cross validation, 80 percent of observations were iteratively used for building the models and the remained 20 percent were used for validation. It is worth to note that in RF, as OOB estimate of the error rate is an unbiased estimate of the generalization error, it is not necessary to test the predictive ability of the model using a cross validation procedure [START_REF] Breiman | Random Forests[END_REF]. However, in accordance with MLR and ANN and for a more reliable comparison, a 5-fold cross validation was performed. A number of statistics was calculated between predicted parameter from GLAS data (maximum height, Lorey's height or volume) and correspondent in situ measurements. Adjusted coefficient of determination (R 2 a ) as an indicator of the fit quality [START_REF] Cameron | R-Squared measures for count data regression models with applications to health care utilization[END_REF], Root Mean Square Error (RMSE) as a measure of accuracy [START_REF] Lee | Physically based vertical vegetation structure retrieval from ICESat data: validation using LVIS in White Mountain National Forest, New Hampshire, USA[END_REF], Mean Absolute Error (MAE) as a measure of dispersion [START_REF] Chai | Root mean square error (RMSE) or mean absolute error (MAE)? -Arguments against avoiding RMSE in the literature[END_REF][START_REF] Willmott | Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance[END_REF], Mean Absolute Percentage Error (MAPE) as an expression of accuracy in percentage [START_REF] Makridakis | Evaluating Accuracy (or Error) Measures[END_REF][START_REF] Hyndman | Another look at measures of forecast accuracy[END_REF], and Akaike Information Criterion (AIC) as a means for model selection by trading-off between the goodness of fit of the model and the complexity of the model [START_REF] Burnham | Information and likelihood theory: a basis for model selection and inference, in Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach[END_REF] were used to evaluate the result of predictions. The significance of all models was also tested statistically. All above mentioned statistics were presented in equations 3.3 to 3.7.

(3.3) Š‹OE• = Ž [ •% • : % ' • ' " " •"+ • (3.4) Š K " = Š " : - -× 91 : Š " C • : -: 1 ˜$$$$$$$$$$Š " = 1 : [ •% • : % ' • ' " " •"+ [ 9% • : % ™ C " " •"+ (3.5) ‹š• = [ ›% • : % ' • › " •"+ • (3.6) ‹šoe• = 1 • • ž •% • : % ' • ' % • ž " •"+ × 100 (3.7) šŸo = 2 4 •¡¢ƒ$ 9ŠOEOE£•C¤
Where Š " : determination coefficient, n: number of observations, % • : in situ height for the plot i, % ' • : estimated height for the plot i, % ™ : mean in situ height, : number of predictors in the model, ŠOEOE: residual sum of squares for the fitted model.

Production of canopy height map and its validation

Since GLAS data does not provide continues coverage of the study area, other sources of remotely sensed data were employed to produce a wall to wall height map. Two strategies were used to produce canopy height map. The first one provides a height map based on regression model regardless spatial correlation between the canopy heights. The second one is regression-kriging, a spatial prediction technique that combines a regression of the dependent variable on some predictors with kriging of the regression residuals. The resulted maps were validated using 32 in situ plots dispread over part of study area, accidently.

Canopy height map using regression model

Following steps describe the procedure of canopy height map production using regression model:

1) Applying GLAS height model on all GLAS data: Since canopy height models were built and validated using 60 GLAS waveforms (called GLAS height models hereinafter), the best maximum height and Lorey's height models were employed to estimate these heights at the location of all GLAS footprints over the study area (450 footprints).

2) Developing height model using GLAS heights as reference data and indices extracted from other remote sensing (radar and optical images) and environmental data: As it was mentioned before, to match all indices extracted from PALSAR and optical images (TM and SPOT) and also DEM-extracted variables spatially, they were resampled to 10 m resolution. Since the average size of GLAS footprints is 70 meter in diameter, mean value of all indices in a 7 by 7 window at the center of GLAS footprints were calculated (mode value for categorical variables). Multiple linear regression (MLR) and Random Forest regression (RF) was used to develop canopy height models which is called second height model hereinafter). For developing MLR models, the most correlated indices were entered in a stepwise regression. Selection of indices for RF models was based on both stepwise regression and importance degree of indices. The main advantage of random forest is its incorporation of continuous or qualitative predictors without making assumptions about their statistical distribution or covariance structure [START_REF] Breiman | Random Forests[END_REF].

3) Selection of best second height model based on statistical criteria through cross validation described in section 3.3. 4) Applying the best second height model (For both H max and H Lorey ) on the study area and producing maximum and Lorey's height map.

Canopy height map using regression-kriging method

Kriging is an advanced geostatistical procedure of interpolation that generates an estimated surface from a scattered set of points with z-values. Kriging assumes that the distance or direction between sample points reflects a spatial correlation that can be used to explain variation in the surface. It is a multistep process which includes exploratory statistical analysis of the data, variogram modeling, creating the surface, and (optionally) exploring a variance surface [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF].

Kriging weights the surrounding measured values to derive a prediction for an unmeasured location. The general formula is formed as a weighted sum of the data (Equation 3.8):

¥9OE m C = $ [ ¦ • §9OE • C " •"+ (3.8)
z(S i ): the measured value at the i th location, z(S 0 ): predicted value at prediction location S 0 , n: the number of measured values and λ i : an unknown weight for the measured value at the i th location. In ordinary kriging which is used in this study, λ i depends on a fitted model to the measured points (the fitted semivariogram function), the distance between the measured points and the prediction location and also the overall spatial arrangement of the measured points.

Regression-kriging involves spatially interpolating the residuals from a non-spatial model using kriging, and adding the results to the prediction obtained from the non-spatial model [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF]. As described in section 3.4.1, second height model was built to predict canopy height from remotely sensed and environmental data, and a wall to wall canopy height map was produced by applying this model on the entire study area. But it does not take into account the spatial correlation between the canopy heights. The main objective of this section is to consider spatial correlation between the canopy heights in order to improve height map. To do so, semivariogram analysis was applied to the regression residuals (the difference between the predicted and actual values) to quantify the spatial structure of canopy height. This method has been widely used to analyze spatial structures in ecology [START_REF] Robertson | Geostatistics in ecology: Interpolating with known variance[END_REF][START_REF] Ge | Regression-kriging for characterizing soils with remote sensing data[END_REF][START_REF] Eldeiry | Comparison of regression kriging and cokriging techniques to estimate soil salinity using Landsat images[END_REF]. The semivariogram plots the semivariance as a function of the distance between samples using equation 3.9:

i9¨C = $ + "©9ªC [ ¡ §9OE • C : §9OE • 4 ¨C¤ " ©9ªC •"+ (3.9)
Where$«9¬C is semivariance as a function of distance h, z(S i ) and z(S i +h) are the estimated residuals from the regression models at locations S i and S i +h, a location separated by distance h, N(h) is the total number of pairs of samples separated by distance h. The Where S 2 is the nugget, σ 2 the sill, and a the range of the semivariogram, g(h).

Next, the fitted semivariograms were used in the kriging of the canopy height residuals for maximum and Lorey's height using equation 3.8 and then defining the regression-kriging model. As its name indicates, regression-kriging is consisted of a regression part (± ² 9OE m C) and a kriging part (¥9OE m C) as shown in the equation 3.11.

³9t m C = $ " ² 9OE m C 4 $¥9OE m C = $ " ² 9OE m C 4 $ [ ¦ ´ §9OE ´C • ´=1 (3.11)
Where ³9t m C is canopy height value using regression-kriging method, ± ² 9OE m C is the fitted trend, ¥9OE m C$is the kriged residual, λi are the kriging weights determined by the spatial dependence structure of the residual and z(S i ) is the residual at location S i .

Results and discussion

Estimation of maximum canopy height (H max ) using ICESat GLAS

Direct method for estimation of H max

Maximum canopy height calculated using direct method (vertical difference between signals start and ground peak of GLAS waveform) was compared with in situ H max . The correlation between estimated H max and H max -total was higher rather than the correlation between estimated H max and H max -in or H max -out (refer to section 2.3.1). Figure 4.1 shows estimated H max versus H max -total. The adjusted coefficient of determination (Š K " ) and root mean square error (RMSE) are 0.48 and 9.9, respectively.

Fig. 4.1.

Estimated H max from GLAS data using direct method versus in situ H max

Estimation of H max using regression models (MLR, RF and ANN)

In order to reduce the impact of slope on estimation of canopy height, three statistical methods such as multiple linear regressions (MLR), random forest (RF) and artificial neural network (ANN) were used.

Estimation of H max using MLR

Concerning maximum height, the result of regression models using in situ H max -total was better than H max -in and H max -out. Thus hereinafter to the end of GLAS height models result, H max refers to H max -out. contributed most to the model for this set of independent variables.

As it is seen in table 4.1, the accuracy of the simplest model (RMSE = 6.3 m) based on two metrics (W ext and TI 10 ), was about one meter lower than the accuracy of model 1 (RMSE = 5.0 m). But it should be noticed that the predictor H 50 in model 1 is dependent on ground peak identification. So there is uncertainty in extraction of this metric.

However models containing H trail and H lead produced good result somewhat, but t-statistic showed these predictors are not significant.

As it is observed in figures 4.2 and 4.3, overestimation and underestimation were decreased largely using regression models 1 and 2 in comparison with direct method (refer to Figure 4.1). The use of TI 90 instead of TI 10 in models 1 and 2 led to an R 2 a of 0.83 and 0.72, and RMSE of 5.3 m and 6.9 m, respectively (Equation 4.1 and 4.2). Figure 4.4 shows the estimated H max using equations 4.1 and 4.2 versus in situ H max . As it is seen, using TI 90 instead of TI 10 did not have significant effect on the result of model. Similar outputs were obtained The model offered by stepwise regression (Equation 4.3) produced an R 2 a and RMSE of 0.82 and 5.1 m, respectively (Fig. 4.5). As it is seen, the result of stepwise regression is similar to model 1 from table 4.1. But based on principle of parsimony which tends to prevent overfitting and reduction in prediction ability of the model, a model should be simple as much as possible [START_REF] Vandekerckhove | Model comparison and the principle of parsimony[END_REF]. So however the model resulted from stepwise regression produced good result, it is not preferable because of its complexity.

H max = -0.04336 W ext 2 + 0.41997 W ext . It produced an R 2 a and RMSE of 0.77 and 6.0 m, respectively, the MAE between predicted and observed height was about 4.7 m, and the prediction error was about 22.1% (Fig. 4.6). Based on the t-statistics of regression coefficients, TI 10 1.5 and W ext 1.5 contributed most to the model. In situ H max (m) R 2 = 0.74 RMSE= 6.8 m 1:1 Table 4.4 shows statistics of three RF models based on PCs. The best result concerning RF regressions based on PCs was generated using the same metrics as MLR (three first PCs, W ext and TI 10 ) with an R 2 a and RMSE of 0.66 and 8.0 m, respectively (Fig. 4.8). As it is seen, however PC 1 consists the largest variance of data among all PCs, it has less importance degree rather than PC 2 and PC 3 in the model. This confirms that the informative part of waveform is not always in the first PC.

Estimation of H max using ANN

In order to access the optimal structure of ANN, numerous networks with different hidden layers and neurons and also various iteration rates were assessed. As it is seen in table 4.5, three-layer networks (one hidden layer) performed well in predicting forest height. Based on the results, a three layers network with only two metrics "W ext " and "TI 10 " is able to predict maximum height with an RMSE and R ! " of 5.7 m and 0.80, respectively (Fig. 4.9).

Adding other metrics did not improve the result considerably. As it is seen, an ANN with three metrics "W ext ", "TI 10 " and "H 50 " produced an RMSE and R ! " of 5.4 m and 0.82, respectively (Fig. 4.10).

Discussion on H max estimated using GLAS data

Maximum canopy height retrieved from direct method did not show a good result. Canopy height has been overestimated where there are short trees and mostly under-estimated in tall trees locations (Fig. 4.1). Overestimation is expected especially where short trees are located over a sloped terrain. In these conditions, the elevation of the highest object within a footprint is not necessarily at the top of the tallest tree, and could be a shorter tree located in higher elevation or even terrain instead of any vegetation which could occur for sparse canopy over steep terrain (Chen, 2010b). Deep investigation in our field data confirms footprints possessing short trees are located over a sloped terrain (the range of terrain slope for these footprints except one (20%) is between 40-55%) with low forest volume as a proxy of forest density (plots 22-28 and 30 in figure 2.10b). As it was demonstrated in figure 2.10a, number of trees (by ha) is approximately like most plots (except for one plot), but correspondent basal area are very low rather than other plots. The slope problem has been solved greatly using an MLR model combining terrain information with GLAS's waveform metrics (W ext 2.5 , W ext 1.5 , TI 10 1.5 , Ln(H 50 ) (Fig. 4.2) and also an ANN model employing three waveform metrics (W ext , TI 10 , H 50 ) (Fig. 4.10). The over and underestimation of height has been also decreased considerable. Generally, based on our outcomes, predicting height of short trees using GLAS data is difficult. This was also reported by [START_REF] Nelson | Model effects on GLAS-based regional estimates of forest biomass and carbon[END_REF]. He showed lack of efficiency of GLAS data to accurately measure forest structure in short-tree sparse forests. In total, all three regression methods (MLR, RF, ANN) based on waveform metrics produced greater accuracy in comparison with models based on PCs. The performance of the best MLR, RF and ANN models based on waveform metrics was compared in figure 4.12. As it is observed, RF had the weakest performance specially where there are short trees (plots 22-28 and 30). In terms of using TI 90 instead of TI 10 , models with TI 10 produced just slightly better results. This is contrary to our expectations for producing much more accurate result using local DEM generated from topographic map rather than SRTM DEM. One reason could be that conventional DEMs produced from photogrammetric techniques might not adequately characterize topography over forest areas (NOAA Coastal Services Center, 2012). Conclusively, the SRTM DEM could be an acceptable source of information about terrain variability especially in large extent areas with presence of forest cover. Recent availability to the SRTM DEM 30 for whole world (with more details rather than SRTM DEM 90 ) strengthens this deduction. 4.7 represents some regression models predicting H Lorey . The first model including ln(W ext ) and TI 10 produced the lowest AIC (288.3) with a prediction error of about 24.0% and RMSE of 5.1 m (Fig. 4.13). As mentioned before, significance of coefficients was considered by calculation of t-statistics. However models containing H 50 , H trail and H lead (models 4 and 5) produced good result somewhat. But based on t-statistic criterion, coefficients of these metrics are not significant. To estimate Lorey's height based on PCs, like H max , all PCs or PCs from stepwise regression (20 PCs) were used for prediction. The resulted models produced high error (the RMSE greater than 25 m). MLR models based on only three first PCs, containing 75% of data variance, produced better result (Table 4.8). However based on t-statistics, some coefficients in the model 2 and 3 were not significant.

Estimation of H Lorey using RF

The five best RF models developed based on waveform metrics for estimation of H Lorey were presented in table 4.9. Generally, the four first models produced approximately the same result. Model 1 with four variables (metrics) including Ln(H 50 ), TI 10 1.5 , W ext 2.5 and Ln(W ext ) produced an RMSE of 5.4 m. Statistic of mean absolute percentage error (MAPE)

shows 26.9% of predictions of this model were off. As it is seen, all models presented in this table includes TI 10 or TI 10 1.5 , which indicate the importance of this variable in estimation of Lorey's height over sloped area. The prediction error of model 5 built using metrics W ext , TI, H trail and H lead is greater than the first four models. In fact models containing H trail or H lead showed less performance in comparison with the others. It could be because of uncertainties in extraction of H trail and H lead form waveforms broadened by terrain slope [START_REF] Lefsky | Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveform[END_REF]. Figure 4.14 shows estimated H Lorey using model 1 and model 5 versus in situ H Lorey . As it is seen overestimation and underestimation in model 1 is lower than model 5. which is similar to the result of RF models based on waveform metrics (Table 4.9). In these models, PC 3 has more contribution to the model than PC 2 and PC 1 , however it contains less variance of the data. The optimal structure of ANN was obtained by trial and error. The results showed a three layer network produces suitable outcome in this research. Generally, H Lorey was predicted with an accuracy of about 5 meters using many types of ANN. 4.12). Adding W ext as input variable improved the result significantly. In other word, an ANN with four inputs including three first PCs and W ext (model 2, Table 4.12) estimated H Lorey with higher accuracy (RMSE = 3.4 m and 2 a R = 0.87) (Fig. 4.17). 12.3% of predictions by this model are off from true measurements. Adding TI 10 did not improve the result comparatively to the model 2.

Discussion on H Lorey estimated using GLAS data

MLR, RF and ANN had approximately similar performance in terms of employing waveform metrics as predictors. In other words, all three methods based on waveform metrics are able to predict Lorey's height with an accuracy of about 5 meters. ANN outperformed two other methods when PCs were used as predictors. The superiority of ANN is considerable in case of using PCs as input variables. An ANN model with four input neurons including three first PCs of PCA and W est predicted H Lorey with an accuracy of 3.4 m. The interesting points concerning this model are: firstly, using PCs as input variables which are lack of uncertainties unlike some waveform metrics, especially over sloped terrain; secondly, achieving higher accuracy rather than other models without entering any ancillary data (DEM). In other words, all best models resulted from three statistical methods based on waveform metrics and PCs except ANN based on PCs contained TI extracted from DEM as representative of topography status. Thirdly, in contrast to other models, this model was able to estimate H Lorey properly even in short sparse and tall dense stands. Figure 4.18 demonstrates performance of the best MLR, RF and ANN models based on PCs.

Estimation of forest volume using ICESat GLAS

As it was explained in section 3.2, two methods were applied to estimate forest volume (V).The first method consists of three steps: 1) developing volume-H max and volume-H Lorey relationships ({ = ˆ. U ‰ ); 2) estimating height from GLAS data using best model resulted from statistical methods explained in section 3.1.2.; and 3) estimating forest volume (V) using chosen volume-height relationship. The second method was estimation of forest volume directly from GLAS waveforms. Regards to the first method, it was necessary to find out correlation between forest volume and height. Since two heights; maximum height (H max ) and mean Lorey's height (H Lorey ), were measured in this research, volume-H max and volume-H Lorey relationships were developed based on in situ measurements. As it is observed in figure 4.19, there is stronger correlation between volume and H Lorey rather than H max . And the accuracy resulted from volume-H Lorey (# = 2.6507$% &'()* +.,-/-) is higher than volume-H max (# = 2.65µ$% SKT +.-" ¶ ).

Discussion on volume estimated using GLAS data

As observed, forest volume was estimated with an accuracy of 116.3 m 3 /ha (RMSE%= 25.4) based on volume-H Lorey relationship, and 119.9 m 3 /ha (RMSE%= 26.2) using a neural network model. Concerning volume-H Lorey method, there were several sources of error that resulted propagation of error and low accuracy of volume estimation. Two main sources were: 1) Height-DBH relationships. As known, Lorey's height is not a directly measured parameter on the ground, but is a weighted mean of height. So, to calculate this parameter, height of all trees in each plot were estimated using developed height-DBH relationships.

2) It is known that forest volume is a function of both height and diameter as two essential quantitative factors. But in this research, only third dimension of objects is retrievable from GLAS data. So we built # = 2.6507$% &'()* +.,-/-with an RMSE and R 2 a of 106.6 m 3 /ha and 0.80, respectively. In other words, even with precise estimation of H Lorey from GLAS data, we would expect error about ±100 m 3 /ha in prediction of volume. Deep investigation in field inventory data shows the possibility of having same Lorey's height for completely different forest structure which leads to a different forest volume. To better understanding, three couple of plots were compared in terms of Lorey's height (m), number of trees (n/ha) and volume (m 3 /ha) (Table 4.13). As shown, plots with approximately the same Lorey's height have different volumes. It confirms that estimating forest volume only relying on an average height could cause a high discrepancy with reality especially in uneven aged forests. Concerning RF, the most important variables describing maximum canopy height with lowest error were terrain index (TI), vegetation indices extracted from Landsat-TM and SPOT5 data (min-rvi, mean-summer-ndvi, max-mvi, mvi from dates October2009, June2010 and April2015), mean index from texture analysis of optical images and also texture indices from HH and HV polarization of PALSAR data including mean, correlation and variance (hh-variance, hv-correlation, hv-mean). The RMSE and R ! " were 7.4 m and 0.50, respectively. Figure 4.26 shows estimated maximum height using the best random forest regression based on above variables versus reference maximum height derived from GLAS (GLAS H max ). This model overestimated maximum heights less than 25 m that could be as result of several limitations in height estimation. It should be noticed that he GLAS-based heights were obtained using local GLAS height models developed for a small part of the study site which will lead to height discrepancy especially in heterogeneous forests. The attempt of slope correction in this study was parameterized using only 60 field plots that do not represent all slopes conditions properly, and it could cause the broadening effect of slope on the GLAS waveforms unsolved which will lead to final estimation of height incorrectly.

Canopy height map using regression-kriging

Regression-kriging method was also used to produce canopy height map with considering spatial correlation between canopy heights. The semivariogram of height residuals (maximum height and Lorey's height) and exponential model fitted on them is observed in figure 4.32. The Nugget, Psill and range obtained from semivariogram of H max residuals were 38.35641m 2 , 14.84780m 2 and 2424.932m, respectively. These coefficients for H Lorey were 17.107786m 2 , 6.350782m 2 and 1664.844m. As known, the nugget effect can be attributed to measurement errors or spatial sources of variation at microscales smaller than the sampling interval (or both). In the presented variograms, the nugget effects look too high that could be as a consequence of sparse data. However the residuals of RF models did not exhibit a strong correlation structure (refer to the figure 4.32), regression-kriging method was under consideration to investigate the probability of improving height maps produced using RF models. Thus, kriging of height residuals was performed using kriging weights calculated based on information derived from height residual's semivariogram. This layer was added to the height map produced using RF regression. The resulted height maps were validated using in situ heights on 32 plots in part of study area. The RMSE and R 2 were 5.9 m and 0.72, respectively, for kriged H max and 4.3 m and 0.54, respectively, for kriged H Lorey , Figures 4.33 feature from GLCM analysis on mean-summer-NDVI that was highly correlated with forest height individually, showed also high importance degree in the maximum height RF model. Among PALSAR texture features, "mean", "correlation" and "variance" contributed in the height regression model. However the "variance" was not correlated with forest height individually but it was important for maximum height estimation in combination with other variables. Among environmental data, only TI (terrain index) as a proxy of topography condition was contributed into the model. Other variables including aspect, elevation and geological units did not affect the efficiency of the model noticeably. It is likely to achieve better result if the above environmental data is used in the GLAS-height model (first part of this study). For mean Lorey's height, PALSAR extracted texture features did not show contribution in the best RF model. Among vegetation indices derived from optical data, only NDVI and RVI, and among texture features derived from mean-summer-NDVI, "mean" and "homogeneity" showed high importance degree. Generally TI and indices extracted from optical images had most contribution in the estimation of forest mean and maximum heights. However the RMSE in estimation of H Lorey (5.5 m) is less than H max (7.4 m), both H Lorey and H max are predictable with similar relative error in the range of 22.8 to 23.5% using random forest regression. These models were used to produce a wall to wall height map and estimate H max and H Lorey at locations with no GLAS coverage. In total, based on general knowledge about our forest site and visual interpretation, the resulted maps seem logical and reliable at least for large scale studies. Comparison of predicted heights extracted from H max and H Lorey maps with in situ measurements in small part of study area showed improvement in RMSE for about 2.1 m and 1.2 m for H max and H Lorey , respectively.

In order to improve H max and H Lorey maps, spatial correlation of heights were considered using depiction of height residual semivariogram. Refer to figure 4.32, there was no strong spatial correlation between height residuals. The model fitted on the semivariogram was quite flat that could be a consequence of high nugget effect which itself would be resulted from error in observation data (here GLAS-based heights) or low density of data. So, the height maps resulted from regression-kriging did not improve height accuracy.

Conclusion and perspectives

Conclusion

Measuring biophysical parameters of forest and providing accurate information and knowledge in different scales are vital for forestry plans and ecosystem management.

Remote sensing techniques can provide a less expensive and relatively precise outputs in comparison with field measurements especially for large heterogeneous forests and inaccessible area. This research includes two main parts and aimed to: 1) investigate the capability of GLAS data in estimating forest biophysical parameters including maximum canopy height (H max ), Lorey's height (H Lorey ) and forest volume (V), 2) producing forest height/volume map using integration of ICESat/GLAS, ALOS/PALSAR and optical images and environmental data.

Concerning the first part, numerous MLR, RF and ANN regressions were developed using different sets of metrics including waveform metrics and PCs to estimate each parameter.

In situ measurements were carried out to build and validate regression models. In order to overcome slope effect on GLAS waveforms, terrain index (TI) was derived from digital elevation models provided from topographic maps (DEM 10 ) and SRTM (DEM 90 ) and used as a predictive in regression models.

Regards to the second part, GLAS height models were used to derive H max and H Lorey at the location of all GLAS data in the study area. Then several MLR and RF regression models were developed between GLAS-based heights as reference and indices extracted from PALSAR, Landsat-TM, SPOT5 and DEM to provide second height models. The resulted models were employed for production of height maps. Consequently regression-kriging procedure was implemented to consider the possibility of height values improvement.

Prediction of maximum height (H max ) using GLAS

As expected, H max estimated from direct method (vertical distance between the signal start and the ground peak) did not match highly the real heights (RMSE=9.9m, RMSE% =27.5). It could be a cause of sloped terrain and misidentification of ground peak witch has been also stated by several researches [START_REF] Lefsky | Estimates of forest canopy height and aboveground biomass using lCESat[END_REF][START_REF] Pang | Temperate forest height estimation performance using ICESat GLAS data from different observation periods[END_REF]Chen, 2010b;[START_REF] Xing | An improved method for estimating forest canopy height using ICESat-GLAS full waveform data over sloping terrain: A case study in Changbai mountains, China[END_REF].

In order to decrease effects of terrain slope on estimation of canopy height, different regression methods were employed. An MLR model combining terrain information with GLAS's waveform metrics (W ext 2.5 , W ext 1.5 , TI 10 1.5 , Ln(H 50 ) produced H max with an accuracy of 5.0 m (RMSE% = 13.8

). An ANN model employing three metrics (W ext , TI 10 , H 50 ) estimated H max with an accuracy of 5.4 m (RMSE% = 15). The point regards to these two models is employing H 50 witch is more exposed to uncertainty rather than two other metrics. It is because the extraction of this metric depends on recognition of ground peak witch is difficult on steep area and could contain error. An ANN model using only W ext and TI 10 led to an accuracy of 5.7 m (RMSE% = 16). As observed, these models outperformed direct method (with 13 percent reduction in error of estimation), and over/underestimation of height has been decreased considerably. Although overestimation is still observed in short-trees sparse forest stands ( < 10 m). [START_REF] Nelson | Model effects on GLAS-based regional estimates of forest biomass and carbon[END_REF] has also showed lack of efficiency of GLAS data to accurately measure forest structure in such forests.

Generally models contained TI showed better performance in estimation of forest canopy height. It indicates that TI derived from DEM neutralizes greatly the negative effect of terrain slope on waveform's characteristics. In a study by Chen (2010a) which has been done in three different sites with terrain slope of 20 degrees on average, a simple linear model contains predictors W ext and TI 10 outperformed two other linear and non-linear models containing W ext , H lead and H trail .

PCs-based models (MLR, RF and ANN) did not perform as well as models based on waveform metrics. In overall, they produced better result when models include W ext and TI in addition to PCs.

Comparison of three statistical methods in estimation of maximum canopy height based on waveform metrics indicated that MLR and RF represented respectively, the best and the worst performance. When the regressions were developed based on PCs, the ANN produced, slightly, better result rather than MLR, and RF did net show good result. These results were in contrast with [START_REF] Fayad | Canopy height estimation in French Guiana with lidar ICESat/GLAS data using principal component analysis and random forest regressions[END_REF] that observed approximately the same accuracies in predicting canopy height using MLR or RF models, also waveform metrics or PCs based models. This confirms the local applicability of fitted regressions. It is worth to notice that terrain topography in this research (Planted eucalyptus forests in French Guiana) is mostly flat.

Prediction of Lorey's height (H Lorey ) using GLAS

A simple three layer ANN model using three first PCs of PCA and W est predicted H Lorey with an accuracy of 3.4 m (RMSE% = 12.4). While TI was an important variable in most GLAS height models, this model was able to estimate H Lorey with high accuracy without necessity of participation the terrain information. Generally, ANN showed better performance in comparison with MLR and RF when PCs were used as input variables and RF produced less accurate result comparing MLR. All three methods (MLR, RF and ANN) had approximately similar performance in terms of employing waveform metrics as predictors (RMSE about 5 m).

Effect of terrain index on estimation of height using GLAS

As discussed before, terrain index (TI) was extracted from two sources of DEM: DEM provided from 10m-topographic maps (TI 10 ), and SRTM90 DEM (TI 90 ). Generally models containing TI produced higher accuracy. This metric showed highest correlation with canopy height after W ext , and was important in reduction of broadening effect of terrain slope on waveform metrics.

The result of regression models showed that models containing TI 10 performed slightly better than those including TI 90 . This is contrary to our expectations for producing much more accurate result using local DEM generated from topographic map rather than SRTM DEM. One reason could be that conventional DEMs produced from photogrammetric techniques might not adequately characterize topography over forest areas (NOAA, 2012).

Conclusively, the SRTM DEM may be an acceptable source of information about terrain variability especially in large extent areas with presence of forest cover. Recent availability to the SRTM DEM30 for whole world (with more details rather than SRTM DEM90) strengthens this deduction. However, it is expected to reach higher accuracy using DEM derived from airborne lidar data which has been confirmed by Chen (2010).

Prediction of forest volume (V) using GLAS

Concerning volume predictions, two approaches was employed. The first, estimation of volume using volume-height relationship and the second, volume estimation using regressions developed between in situ volume and lidar based metrics. The result of volume-H Lorey (116.3 m 3 /ha) was slightly better than PCs-based ANN model (119.9 m 3 /ha) but ANN model performed better in very low ( <10 m 3 /ha) and very high ( > 800 m 3 /ha) volume stands. In total, the relative error of forest volume estimated using GLAS data was about 26%. The result of this part is better than findings of [START_REF] Nelson | Estimating Siberian timber volume using MODIS and ICESat/GLAS[END_REF], the only study on estimating forest volume using GLAS data. They predicted timber volume in central Siberia dominated by coniferous with an R 2 a of 0.75 and RMSE of 87 m 3 /ha using a neural network model employing six metrics extracted from GLAS waveform (n=51): h ̄med : a median height which below that cumulative canopy height profile (CHP) is 50% at maximum; h 2-sun : a corrected maximum height; h g1-sun : height of waveform peak with the maximum amplitude above ground peak, f: the slope of the line formed by connecting the signal start point with the peak of the uppermost Gaussian return, r g3 : the waveform area under the 3 rd Gaussian peak, and n g : the number of Gaussian peaks in the waveform. Volume average calculated using GLAS/MODIS considering GLAS shots on all slopes (n=66119) was about 172 m 3 /ha which is about 2.5 times less than mean volume in our study site (mean of 450 m 3 /ha in 60 plots). This led to higher relative error, "percentage of RMSE divided by mean volume", in comparison with our result. In total, however the developed ANN model improved the accuracy of volume estimations in the extreme low and high volume stands, the residual error is still high in such area. It is worth to notice that our findings are based on only 60 plots mostly (about 75%) over stands with volume ranging from 100 to700 m 3 /ha (figure 2.11, Section 2.3.1). Thus, because of low number of plots in very low and high density stands, the network may not be well trained to learn to differentiate. It is needed to investigate scrupulously while increasing the number of observations, applying other statistical methods or participating ancillary data to enhance the accuracy of forest volume estimation in the future studies.

Following section discusses on possible sources of uncertainties in estimation of forest height and volume.

Uncertainties in prediction of height and volume using GLAS

Two main sources of uncertainties in estimation of forest height and volume using ICESat/GLAS were identified. The first one is related to the field data collection, and the other one concerning errors in extraction of waveform metrics.

The general sources of uncertainty in field data could be:

-Time interval between lidar data acquisition and field measurements which has been ignored in the present study because of having deciduous species with low growth rate in their climax age.

-In situ measurement uncertainty which mainly concerns the uncertainty on trees height due to different factors such as measuring tool, measurement procedure, skill of operator and site geography. Aside from the obvious errors associated with wrong measurements of distances or misreading the angles of top and base of trees with the clinometer, there are several less apparent sources of error that can compromise the accuracy of the tree height calculations. An error occurs where 1) the treetop is offset from the base of the tree, or 2) the top of the tree has been misidentified. [START_REF] Larjavaara | Measuring tree height: a quantitative comparison of two common field methods in a moist tropical forest[END_REF] compared tree height measured using tangent method (clinometer) with actual heights using towers

Production of canopy height map

As explained, second height regressions were developed using GLAS-based heights as reference and other remotely sensed and environmental data. In total, RF models outperformed the MLR models in predicting maximum (RMSE = 7.4 m, Š K " $= 0.52) and mean Lorey's height (RMSE = 5.5 m, Š K " = 0.60) in this stage. Concerning H max map produced from RF model, however the accuracy is relatively suitable (RMSE = 5.3 m), comparison of heights extracted from the map with field measurements showed overestimation in samples with maximum height ranging from 23 to 35 m. It is worth to note that the total number of plots included in the assessment is 32 with H max ranging from 23 to 48 m (except one plot = 15 m) and lack of samples especially in short sparse stands. Furthermore, the total area represented by the in situ plots is very small as compared to the total area mapped. In total only 3.2 ha of the 15000 ha of forested lands were sampled. Regards to mean Lorey's height, because of lower range of values, the resulted H Lorey map showed higher accuracy equal to 4.3 m rather than H max . But the relative error for H max and H Lorey is similar equal to 14.8% and 14.4%, respectively. An advantage of production of H Lorey map rather H max map is inclusion of less data sources. As observed before, the best regression model describing maximum canopy height includes indices extracted from DEM, PALAR and optical images (Landsat-TM and SPOT5). It is in accordance with research done in east part of Hyrcanian forests of Iran by [START_REF] Attarchi | Improving the estimation of above ground biomass using dual polarimetric PALSAR and ETM+ data in the hyrcanian mountain forest (Iran)[END_REF]. They found that the joint use of optical and SAR data increases the reliability of the biomass model, significantly. But for Lorey's height, PALSAR extracted indices did not have contribution in the best regression.

The attempt for improving the precision of canopy height map using regression-kriging was unsuccessful in contrast with the result achieved by [START_REF] Fayad | Regional scale rain-forest height mapping using regression-kriging of spaceborne and airborne lidar data: application on French Guiana[END_REF] that reported an improvement of about two meters in terms of RMSE for forest canopy height map. In our study, the exponential model fitted on the height residual semivariogram did not show strong spatial correlation which could be as a result of the heterogeneity of the study area in case of forest structure and topography. While [START_REF] Fayad | Regional scale rain-forest height mapping using regression-kriging of spaceborne and airborne lidar data: application on French Guiana[END_REF] worked on relatively homogeneous flat forests.

In total, there are several limitations to production of height map for the study area. The slope correction attempted here was parameterized using only 60 field plots that do not represent all slopes conditions properly, especially in case of steepest slopes (greeter than 60%). The GLAS-based heights were obtained using local GLAS height models developed for a small part of study site which will lead to height discrepancy especially in heterogeneous forests. To ensure that as many GLAS footprints as possible were included in the analysis, the data from all GLAS laser campaigns from September 2003 on were included and processed using the same algorithms. These data were collected over a more than 5-year period, and therefore do not reflect a static moment in time. Lastly, the validation of the height maps is incomplete and is limited by a lack of field observations for many of the forested lands within study site.

Several other studies have used GLAS data to derive canopy height over large regions, but have typically combined them with 250 m resolution MODIS data rather than 10 to 30 m resolution spot5 and Landsat-TM data. These studies typically developed products that are global assessments of canopy height at a coarse scale [START_REF] Lefsky | Global forest canopy height map from the moderate resolution imaging spectroradiometer and the geoscience laser altimeter system[END_REF][START_REF] Simard | Mapping forest canopy height globally with spaceborne lidar[END_REF][START_REF] Fayad | Regional scale rain-forest height mapping using regression-kriging of spaceborne and airborne lidar data: application on French Guiana[END_REF]. In this study, we demonstrated GLAS data are also useful in mapping at finer resolution, although subject to the limitations identified herein. However further work needs to be done to thoroughly understand and quantify the various sources of error underlying the lack of correspondence between the field observations and the mapped canopy height values, this map provides a good understanding of the distribution of forest canopy height across the study area in a short time and at the lowest cost.

General conclusion

In this research, capability of ICESat GLAS was investigated for estimation of forest maximum canopy height, mean Lorey's height and forest volume in part of Hyrcanian forests of Iran. It was also subjected to provide forest canopy height and volume map if it is possible.

Based on the result, GLAS was able to estimate maximum canopy height and Lorey's height with an accuracy of about 5m and 3.4m, respectively. In contrary to the suitable results related to forest heights (maximum and mean), GLAS data did not meet the desirable achievements in estimation of forest volume. These results address the first main question of this study about ability of GLAS data in estimation of forest biophysical parameters in mountainous heterogeneous forests of north of Iran.

As it was observed in the result, two metrics of waveform extent (W ext ) and terrain index (TI) had key role in estimation of GLAS-based heights. W ext as vertical distance between signal start and end of a waveform is directly related to the canopy height over flat area, but it is extended by increasing of the terrain slope. It challenges the height derivation over severe topography. In this research, using regression models and digital elevation model helped to overcome the impact of terrain slope on waveform characteristics. Generally, all regression models containing TI extracted from DEM outperformed models not-including TI. It confirms the importance of this metric in estimation of height using GLAS data which was under question in our study area. Although, further research needs to be done to address the impacts of slope on height recovery in Hyrcanian forests of Iran using GLAS data.

Parametric and non-parametric statistical methods of multiple linear regression, random forest and artificial neural network were employed for GLAS-based height/volume estimation. Concerning the H max , the best result was obtained using MLR based on waveform metrics. An ANN model based on waveform metrics produced also relatively similar result. Regards to H Lorey , PCs based neural network models outperformed two other statistical methods. A neural network model based on PCs performed also better than MLR and RF regressions in estimation of forest volume. Thus random forest showed the weakest performance in GLAS based forest parameter derivation in our study area.

As a consequence of desirable results concerning height estimation (both H max and H Lorey ), production of height maps was under consideration using regression models and regression-kriging method. Thus a synergy of GLAS, PALSAR, Landsat-TM, SPOT5 and environmental data were used in the analysis. Comparison of resulted maps from regression models with field observations (in a small part of the study area near to GLAS footprints locations) showed promising outcomes. However, further work needs to be done to thoroughly understand and quantify the various sources of error underlying the lack of correspondence between the field observations and the mapped canopy height values. The regression-kriging method did not improve the former height maps that could be as a consequence of heterogeneity of the study area. By this end, the second important question of this research related to the possibility of forest height/volume map production using combination of GLAS and other remotely sensed data was answered.

However terrain slope is very important in estimation of forest parameters using GLAS data, consideration of forest characteristics such as forest type and horizontal and vertical structure of forest may influence the quality of estimations. For instance, it is expected to achieve better results in pure even-aged forest that have simpler vertical structure rather than mixed uneven-aged forest.

With the end of the GLAS data collection in 2009 no new spaceborne lidar data are currently available for updating results. The ICESat-2 mission, which will provide new spaceborne lidar data using the Advanced Topographic Laser Altimeter System, is scheduled for launch in 2017.

Perspectives

-The Hyrcanian forest is a unique natural heritage of global importance. The Caspian region harbors the world's last remaining primary forests of the temperate deciduous forest formation. This research was done on a small area with about 15000 ha forest cover of 1.8 million hectare of Hyrcanian forest. Given the complexity of this forest in terms of vertical and horizontal structures, forest types, forest density, topography etc., it is suggested to do supplementary researches using spaceborn lidar in other parts of mixed broadleaved Hyrcanian forests.

-This research took into account the terrain slope in the estimation of canopy height using GLAS data. Since many other forest characteristics (i.e. forest type, mixture percent and forest age) may affect the quality of estimations, it is interesting to consider these factors in the future studies. However, based on the attained experiences, such a detailed research needs denser coverage of lidar data which is addressed in the upcoming ICESat-2 technology.

-GLAS height models were developed only using leaf-on season lidar data and was applied to all time lidar data. This was a limitation of our research because of low density of GLAS data over study area. As stated by [START_REF] Pang | Temperate forest height estimation performance using ICESat GLAS data from different observation periods[END_REF], the summer period GLAS waveforms capture the returns from forest canopy. The data from early stage of autumn period still contain enough returns from forest canopy, even with lower intensity. The spring period and late autumn period data contain less signals from forest canopy and difficult to estimate forest height. Therefore, it is expected to improve estimations by performing the analysis separately using leaf-on and leaf-off season data in the deciduous broadleaved forests.

-For the second part of this research (providing forest height map from synergy of GLAS, PALSAR, optical and environmental data) some environmental data such as geological units, aspect and elevation maps were used in addition to TI and slope as predictors in building second height models. Among these variables, only TI showed good contribution in the models. This predictor was also used in the first part of this study which was detecting capability of GLAS data for retrieving mean and maximum height. It is probable that considering other environmental data rather than TI in the first step (GLAS height models), lead to good contribution of them in the second step (second height model and finally height map). However as mentioned before, such investigations require dense coverage of lidar data so that there be enough lidar samples representing different conditions or classes of an environmental data.

-In order to provide forest canopy height map using regression model (non-spatial method), MLR and RF models were developed between GLAS-based heights and pre-mentioned remotely sensed data. The selected models and consequently the resulted maps showed overestimation especially for maximum height. It could be as a result of inability of GLAS in estimating short trees height ( <10m) which leads to error in reference heights (less than 10m) estimated from GLAS and consequently the fitted RF model based on optical, radar and environmental data overestimates also tree's height shorter than 20 m. But it is also suggested to test other statistical regression methods.

-A limitation of this research was lack of field observations describing all conditions of the study area. Indeed, a larger database (ranges of slopes, height, etc.) would be very useful to better understand the limitations of the proposed methods. Rather than laborious, cost and time consuming field work, some points in a mountainous forest are inaccessible. Therefore, it would be interesting to employ airborne lidar to collect information in such condition. By the way, a comprehensive database of field observations needs also a denser coverage of lidar data which is addressed in the upcoming ICESat-2 technology.

Fig. 1

 1 Fig. 1.1. Major carbon pools and fluxes of the global carbon balance (FAO, 2007) ............. 1 Fig. 2.1. Location of study area in Iran (top right map) and over Landsat image ...................Fig. 2.2. Schematic illustration of location of trees for height measurement in a field plot .............. Fig. 2.3. a) Schematic illustration of height measurement based on trigonometry; b) Calculation of horizontal distance from tree over sloped terrain............................................. Fig. 2.4. illustration of three phases of field work including navigating and locating the plot center, DBH and height measurement ..................................................................................

Fig. 2

 2 Fig. 2.5. Schematic of the relationship between geoidal heights and ellipsoidal heights ..............

Fig. 2

 2 Fig. 2.6. Geological map provided by geological survey and mineral exploitation of Iran .......... Fig. 2.7. Schematic illustration of GLAS instrument operating from ICESat while orbiting the Earth .................................................................................................................................. Fig. 2.8. Color composite of SPOT-5 images over study area ..................................................

Fig. 2 . 9 .

 29 Fig. 2.9. Tree's Height vs. DBH and the best fitted line ......................................................... Fig. 2.10. a) Distribution of number of trees per hectare (n/ha), b) forest volume (m 3 /ha) in 60 reference plots ........................................................................................................................... Fig. 2.11. Histogram of forest volume (m 3 /ha) for 60 reference plots ...................................... Fig. 2.12. Digital elevation model (DEM 10 ) and four extracted maps .................................... Fig. 2.13. Schematic of laser altimeter pulse spreading over a vegetated area, and a returned waveform ....................................................................................................................... Fig. 2.14. GLAS waveform and some metrics over a terrain of 25% slope in the study area. 1ns corresponds to 15cm sampling distance in waveform ............................................. Fig. 2.15. a) A schematic illustration of lidar derived vegetation height (H), H 100 , and waveform extent (W ext ) for one waveform over flat terrain (solid line) and another over sloped terrain (dash line) (Lee et al., 2011). b) Impact of slope on lidar height retrieval (Lee et al., 2011). c,d) A schematic of footprint located over flat (left) and sloped (right) vegetated terrain, respectively ................................................................................................. Fig. 2.16. Information explained by ten first PCs and other PCs ........................................... Fig. 2.17. Gamma naught (γ˚=γ dB ) before (a) and after (b) LEE and multitemporal filters, and the histogram of backscatters over a forested area (green rectangle) (c) before (blue) and after (red) noise reduction. The standard deviation of γ˚ decreases from 1.42 to 0.44 db ................................................................................................................Fig. 2.18. Eight GLCM texture features extracted from PALSAR-HH ..................................Fig. 2.19. NDVI and eight GLCM texture features extracted from it ..................................... Fig. 3.1. Overview of forest canopy height and volume estimation using GLAS data ...........

Fig. 3

 3 Fig. 3.2. Overview of forest canopy height map using combination of GLAS, PALSAR, optical images and environmental data ............................................................................................... Fig. 3.3. A schematic of random forest regression ................................................................. Fig. 3.4. A schematic of MLP neural network ........................................................................ Fig. 3.5. Schematic diagram of a 5-fold cross validation .......................................................... Fig. 3.6. Illustration of nugget, range and sill components of a semivariogram ...................... Fig. 4.1. Estimated H max from GLAS data using direct method versus in situ H max ...............

Fig. 4 .

 4 Fig. 4.18. Comparison of residual errors produced by MLR, RF and ANN for estimation of H Lorey based on PCs ............................................................................................................. Fig. 4.19. Correlation between volume and a) H max , b) H Lorey .................................................

Fig. 4 .

 4 Fig. 4.30. Lorey's height map produced using RF model ....................................................... Fig. 4.31. Comparison of H Lorey map produced from RF model with in situ H Lorey measures at the location of 32 plots ........................................................................................ Fig. 4.32. Examples of fitted semivariograms of a) maximum height and b) Lorey's height residuals ................................................................................................................................... Fig. 4.33. Maximum height map produced using regression-kriging ..................................... Fig. 4.34. Comparison of H max map produced from regression-kriging with in situ H max measures at the location of 32 plots ........................................................................................

Fig. 4 .

 4 Fig. 4.35. Lorey's height map produced using regression-kriging ......................................... Fig. 4.36. Comparison of H Lorey map produced from regression-kriging with in situ H Lorey measures at the location of 32 plots ..............................................................................

  ) and forest volume (V) for 60 in situ plots ...................................................................................................

  based on waveform metrics ..................................................................................................................................... Table 4.4. Statistics of three RF models for estimation of H max based on PCs ........................ Table 4.5. Properties of two ANN models for estimation of H max based on waveform metrics and the resulted statistics ................................................................................................ Table 4.6. Properties of ANN models for estimation of H max based on PCs and the resulted statistics ...................................................................................................................................

  based on waveform metrics ..................................................................................................................................... Table 4.12. Statistics of three ANN models for estimation of H Lorey based on PCs .................... Table 4.13. Comparison of Lorey's height, number of trees and volume in three couple of plots .....................................................................................................................................

Forests

  are of extreme importance to humans in many ways. They are watersheds, and have economic, environmental and climate control benefits. Forests and natural areas play a very important role in maintaining natural processes. Forests are one of the biggest reservoirs of carbon. They contain up to 80 percent of the aboveground carbon in the terrestrial communities and around 33 percent of the belowground carbon. So they help to keep the carbon cycle and other natural processes working and help reduce climate change. The relationship between forests and climate change is complex. On one hand forests can mitigate climate change by absorbing carbon, while on the other they can contribute to climate change if they are degraded or destroyed. In turn climatic changes may lead to forest degradation or loss -which intensifies climate change further. Figure 1.1 presents the different carbon pools and fluxes of the global carbon balance.

Fig. 1 . 1 .

 11 Fig. 1.1. Major carbon pools and fluxes of the global carbon balance[START_REF] Fao | Carbon sequestration in drylands, report on global soil resources 102[END_REF]

Fig. 2 . 1 .

 21 Fig. 2.1. Location of study area in Iran (top right map) and over Landsat image (bottom right map: the pink polygon shows the border of Nowshahr and the blue frame is the border of study area. The left map shows lidar footprints over hillshade of the study area.

Fig. 2 . 2 .

 22 Fig. 2.2. Schematic illustration of location of trees for height measurement in a field plot; 5 trees within a 35 m diameter plot and 5 in a co-center 70 m diameter plot (outer margin of smaller circle)

Fig. 2

 2 Fig. 2.4. illustration of three phases of field work including navigating and locating the plot center, DBH and height measurement

  Fig. 2.6. Geological map provided by geological survey and mineral exploitation of Iran (The white color corresponds to urban and non-forest area)

Fig. 2 . 7 .

 27 Fig. 2.7. Schematic illustration of GLAS instrument operating from ICESat while orbiting the Earth

  5dB(TBR) * Meets under the following off-nadir angle. For under other angle, it is requested to keep similar specification as far as possible. High Resolution Mode: Off-nadir angle 34.3 deg. (TBR) Direct Downlink Mode: Off-nadir angle 34.3 deg. (TBR) ScanSAR Mode: 4th scan (off-nadir 34.1 deg. (TBR)) Polarimetric Mode: Off-nadir angle 21.5 deg. (TBR) Note: Above descriptions are specifications over the equator.

Fig. 2 . 8 .

 28 Fig. 2.8. Color composite of SPOT-5 images over study area.

Fig. 2 . 9 .

 29 Fig. 2.9. Tree's Height vs. DBH and the best fitted line given in table 4-2; a) Fagus orientalis, b) Carpinus betulus, c) Quercus castanifolia, d) Alnus subcordata, e) Group1, and f) Group2

Fig. 2 .

 2 Fig. 2.10. a) Distribution of number of trees per hectare (n/ha), b) forest volume (m 3 /ha) in 60 reference plots

Fig. 2 .

 2 Fig. 2.12. Digital elevation model (DEM 10 ) and four extracted maps (elevation and aspect are classified maps)

Fig. 2 .

 2 Fig. 2.13. Schematic of laser altimeter pulse spreading over a vegetated area, and a returned waveform (NASA, 2015)

Fig. 2 .

 2 Fig. 2.15. a) A schematic illustration of lidar derived vegetation height (H), H 100 , and waveform extent (W ext ) for one waveform over flat terrain (solid line) and another over sloped terrain (dash line) (Lee et al., 2011). b) Impact of slope on lidar height retrieval (Lee et al., 2011). c,d) A schematic of footprint located over flat (left) and sloped (right) vegetated terrain, respectively.

  the number of observations (

Fig. 2 .

 2 Fig. 2.17. Gamma naught (γ˚=γ dB ) before (a) and after (b) LEE and multitemporal filters, and the histogram of backscatters over a forested area (green rectangle) (c) before (blue) and after (red) noise reduction. The standard deviation of γ˚ decreases from 1.42 to 0.44 db.

  this study GLCM measures were derived from both HH and HV polarization. The GLCM characterizes the texture of an image by calculating how often pairs of pixel with specific values and in a specified spatial relationship occur in an image and then statistical measures are extracted from this matrix. Eight GLCM indices such as mean, variance, homogeneity, contrast, dissimilarity, second moment, entropy and correlation were calculated from this matrix. A window size of 3×3, 5×5, 7×7, 9×9 and 11×11 pixels with horizontal and vertical offset of one was used for extraction of texture indices.

Fig. 2 .

 2 Fig. 2.18. Eight GLCM texture features extracted from PALSAR-HH

Fig

  Fig. 2.19. NDVI and eight GLCM texture features extracted from it

Figure 3 .

 3 Figure 3.2 shows the process of producing forest canopy height map using combination of GLAS, PALSAR, optical images and environmental data. In this flowchart also, gray boxes show input data, simple white boxes present data preparation processes and dot boxes indicate outputs. More details are described in the following sections of this chapter.

Fig. 3 . 1 .

 31 Fig. 3.1. Overview of forest canopy height and volume estimation using GLAS data

Fig. 3 . 2 .

 32 Fig. 3.2. Overview of forest canopy height map using combination of GLAS, PALSAR, optical images and environmental data

Figure 3 .

 3 Figure 3.3 illustrates a schematic figure explaining how random forest works.

Fig. 3 . 5 .

 35 Fig. 3.5. Schematic diagram of a 5-fold cross validation

Fig. 3 . 6 .

 36 Fig. 3.6. Illustration of nugget, range and sill components of a semivariogram

  W ext ) and Terrain index calculated from DEM 10 (TI 10 ). AIC (Akaike Information Criterion) was calculated as a way of selecting the best regression model. The model with lowest AIC score represents the best model. An MLR model combined W ext 2.5 , W ext 1.5 , ln(H 50 ) and TI 10 1.5 (model 1) produced the lowest AIC (296.3) and highest accuracy (5.0 m). Based on the MAPE (Mean Absolute Percentage Error), 16.4% of predictions of this model were off (Fig. 4.2). The t-statistics of regression coefficients shows the relative importance of each metric in the model. Based on this statistics,

  Fig. 4.4. Estimated H max using a) model 4.1, b) model 4.2 versus in situ H max

Fig. 4 . 8 .

 48 Fig. 4.8. Estimated H max using RF regression based on PCs (model 3, table 4.4) versus in situ H max

Fig. 4 .

 4 Fig. 4.12. Comparison of residual errors produced by MLR, RF and ANN for estimation of H max based on waveform metrics

Fig. 4 .

 4 Fig. 4.14. Estimated H Lorey using MLR based on PCs; a) model 1, b) model 5 from table 4.9 versus in situ H Lorey

  on PCs performed slightly better than those based on waveform metrics. A simple neural network employing only three first PCs of PCA produced an RMSE and 2 a R of 4.7 m and 0.76, respectively (model 1, Table

Fayad

  et al. (2014) estimated canopy height in forest sites of French Guiana that terrain topography is mostly flat. They obtained approximately the same accuracies using MLR and RF models based on waveform metrics and PCs.

Fig. 4 .

 4 Fig. 4.18. Comparison of residual errors produced by MLR, RF and ANN for estimation of H Lorey based on PCs

Fig. 4 .b

 4 Fig. 4.19. Correlation between volume and a) H max , b) H Lorey

Fig. 4 .

 4 Fig. 4.21. Estimated volume using MLR based on waveform metrics versus in situ volume

Fig. 4 .

 4 Fig. 4.24. Trend of forest volume bias (observed value-estimated value) by a) in situ forest volume, and b) terrain slope (The lowest to highest value on both volume and slope was coded by 1 to 58,respectively. In other words, plot 1 in graph a, as an example, is not necessarily the same plot 1 in graph b.)

Fig. 4 .

 4 Fig. 4.26. Estimated H max using RF regression based on TI and extracted indices from optical and radar data versus reference H max (GLAS H max )

Fig. 4 .Fig. 4 .

 44 Fig. 4.27. Maximum height map produced using RF model

  Fig. 4.32. Examples of fitted semivariograms of a) maximum height and b) Lorey's height residuals

  Fig. 4.33. Maximum height map produced using regression-kriging

  1. Estimated H max from GLAS data using direct method versus in situ H max ............... Fig. 4.2. Estimated H max using MLR based on waveform metrics (model 1, table 4.1) versus in situ H max .................................................................................................................... Fig. 4.3. Estimated H max using MLR based on waveform metrics (model 2, table 4.1) versus in situ H max ....................................................................................................................Fig. 4.14. Estimated H Lorey using MLR based on PCs; a) model 1, b) model 5 from table 4.9 versus in situ H Lorey ................................................................................................................... Fig. 4.15. Estimated H Lorey using RF based on PCs versus in situ H Lorey (a, b, c: model 1 to 3, table 4.10) ....................................................................................................................... Fig. 4.16. Estimated H Lorey using ANN based on waveform metrics (model 1, table 4.11) versus in situ H Lorey .................................................................................................................. Fig. 4.17. Estimated H Lorey using ANN based on waveform metrics (model 1, Table 4.12) versus in situ H Lorey ........................................................................................................

Fig. 4.4. Estimated H max using a) model 4.1, b) model 4.2 versus in situ H max ....................... Fig. 4.5. Estimated H max using stepwise regression of waveform metrics (Eq. 4.3) versus in situ H max ............................................................................................................................... Fig. 4.6. Estimated H max using MLR based on PCs (model 3, table 3.5) versus in situ H max ......................................................................................................................................... Fig. 4.7. Estimated H max using RF regression based on waveform metrics (model 1, table 4.3) versus in situ H max ................................................................................................... Fig. 4.8. Estimated H max using RF regression based on PCs (model 3, table 4.4) versus in situ H max ................................................................................................................................................. Fig. 4.9. Estimated H max using ANN based on waveform metrics (model 1, table 4.5) versus in situ H max ................................................................................................................... Fig. 4.10. Estimated H max using ANN based on waveform metrics (model 2, table 4.5) versus in situ H max ...............................................................................................................................

Fig. 4.11. Estimated H max using ANN based on PCs (model 3, table 4.6) versus in situ H max ............... Fig. 4.12. Comparison of residual errors produced by MLR, RF and ANN for estimation of H max based on waveform metrics ........................................................................................ Fig. 4.13. Estimated H Lorey using MLR based on waveform metrics (model 1, Table 4.7) versus in situ H Lorey ..................................................................................................................

Table 2 .

 2 1. Inventory form for DBH measurement ..................................................................Table 2.2. Inventory form for height measurement ................................................................ Table 2.3. ICESat operational periods .................................................................................... Table 2.4. GLAS data products ............................................................................................... Table 2.5. Main characteristics of ALOS-1 PALSAR (Japan Space Center, 2012) ............... Table 2.6. Landsat-TM and SPOT 5-HRG spectral bands and their resolutions .................... Table 2.7. Distribution of in situ plots in elevation and slope classes ....................................

Table 2 .

 2 8. Non-linear height-diameter functions ....................................................................

Table 2 .

 2 9. Six selected non-linear height-DBH models and their statistical performance ............

Table 2 .

 2 10. Six volume models and their parameters (produced by FRWO) .........................

Table 2 .

 2 11. Statistical summary of maximum height (H max ), Lorey height (H Lorey

Table 2 .

 2 12. The main information derived from GLA01 and GLA14 .........................................

Table 2 .

 2 13. Definition of metrics extracted from GLAS waveforms .....................................

Table 3 .

 3 1. Some MLR models fitted based on waveform metrics ..........................................

Table 4 .

 4 1. Statistics for five MLR to estimate H max based on waveform metrics ...................

Table 4 .

 4 2. Statistics for three MLR models to estimate H max based on PCs ...........................

Table 4 .

 4 3. Statistics of five RF models for estimation of H max

Table 4 .

 4 7. Statistics for five MLR for estimation of H Lorey based on waveform metrics ..................

Table 4 .

 4 8. Statistics of three MLR models for estimation of H Lorey based on PCs ....................

Table 4 .

 4 9. Statistics of five RF models for estimation of H Lorey based on waveform metrics ............

Table 4 .

 4 10. Statistics of three RF models for estimation of H Lorey based on PCs ........................

Table 4 .

 4 11. Statistics of three ANN models for estimation of H Lorey

Table 2 .1. Inventory form for DBH measurement Plot Number: 75 Slope: 60% Aspect: North Elevation: 1950m Date: 12/09/2013 Tree number

 2

		Plot (r = 18m) Small (r = 35m) Large	Species	DBH (cm)	Considerations
	1	ü	Carpinus betulus	24.5
	2	ü	Fagus orientalis	50	It is one of ten top trees in the plot.
	3	ü	Fraxinus excelsior	34
	4	ü	Fagus orientalis	35
	5	ü	Fagus orientalis	33
	6	ü	Fagus orientalis	32
	. . .
	184	ü

Acer campestre 16

The last measured tree in this plot

Table 2

 2

	Tree number	Species	DBH (cm)	Tan of top of tree angle (%)	Tan of base of tree angle (%)	Slope angle (degree)	Distance to the tree (m)	Plot (r = 18m) Small (r = 35m) Large
	١	Fagus orientalis	42	+8	-62	29	27.90	ü
	٢	Fraxinus excelsior	32	+36	-22	10	29.10	ü
	٣	Quercus castanifolia	57	+6	-58	27	24.90	ü
	٤	Fagus orientalis	34	+105	+44	26	38.80	ü
	٥	Fagus orientalis	50	+19	-45	23	26.70	ü
	٦	Fagus orientalis	54	+7	-60	29	34.40	ü
	٧	Quercus castanifolia	65	+23	-34	16	28	ü
	٨	Quercus castanifolia	83	+12	-50	24	28.90	ü
	٩	Fagus orientalis	44	+32	-26	11	32.80	ü
	١٠	Fagus orientalis	42	+25	-556	27	31.60	ü
	١١	Fagus orientalis	31	+105	+53	26	37

.2. Inventory form for height measurement Plot Number: 75 Slope: 60% Aspect: North Elevation: 1950m Date: 12/09/2013 ü Considerations: Tree number 11 is the closest tree to the center of plot

Table 2 .3. ICESat operational periods Laser Identifier Start Date End Date Days in Operation Duration of repeat orbit cycle (days)

 2

	L1A	2003-02-20	2003-03-21	29	8
	L1B	2003-03-21	2003-03-29	9	8
	L2A	2003-09-25	2003-10-04	10	8
	L2A	2003-10-04	2003-11-19	45	91
	L2B	2004-02-17	2004-03-21	34	91
	L2C	2004-05-18	2004-06-21	35	91
	L3A	2004-10-03	2004-11-08	37	91
	L3B	2005-02-17	2005-03-24	36	91
	L3C	2005-05-20	2005-06-23	35	91
	L3D	2005-10-21	2005-11-24	35	91
	L3E	2006-02-22	2006-03-28	34	91
	L3F	2006-05-24	2006-06-26	33	91
	L3G	2006-10-25	2006-11-27	34	91
	L3H	2007-03-12	2007-04-14	34	91
	L3I	2007-10-02	2007-11-05	37	91
	L3J	2008-02-17	2008-03-21	34	91
	L3K	2008-10-04	2008-10-19	16	91
	L2D	2008-11-25	2008-12-17	23	91
	L2E	2009-03-09	2009-04-11	34	91
	L2F	2009-09-30	2009-10-11	12	91

Table 2 .

 2

		4. GLAS data products
	Short name	Long name
	GLA01	GLAS/ICESat L1A Global Altimetry Data
	GLA02	GLAS/ICESat L1A Global Atmosphere Data
	GLA03	GLAS/ICESat L1A Global Engineering Data
	GLA04	GLAS/ICESat L1A Global Laser Pointing Data
	GLA05	GLAS/ICESat L1B Global Waveform-based Range Corrections Data
	GLA06	GLAS/ICESat L1B Global Elevation Data
	GLA07	GLAS/ICESat L1B Global Backscatter Data
	GLA08	GLAS/ICESat L2 Global Planetary Boundary Layer and Elevated Aerosol
	GLA09	GLAS/ICESat L2 Global Cloud Heights for Multi-layer Clouds
	GLA10	GLAS/ICESat L2 Global Aerosol Vertical Structure Data
	GLA11	GLAS/ICESat L2 Global Thin Cloud/Aerosol Optical Depths Data
	GLA12	GLAS/ICESat L2 Antarctic and Greenland Ice Sheet Altimetry Data
	GLA13	GLAS/ICESat L2 Sea Ice Sheet Altimetry Data
	GLA14	GLAS/ICESat L2 Global Land Surface Altimetry Data
	GLA15	GLAS/ICESat L2 Ocean Altimetry Data

Table 2 . 5 .

 25 Main characteristics of ALOS-1 PALSAR (Japan SpaceCenter, 2012)

	Observation Mode	Fine Resolution	Direct Downlink	ScanSAR	Polarimetric	Remarks
					L band
	Polarization	HH or VV	HH + HV or VV + VH	HH or VV	HH or VV	HH + HV+ VV + VH
	Incidence Angle	8˚ -60˚	8˚ -60˚	8˚ -60˚	18˚ -36˚ (3 scan) 18˚ -40˚ (4 scan) 18˚ -43˚ (5 scan)	8˚ -30˚	Off-nadir Angle: 9.7˚ -50.8˚
		Range	10 m *	20 m *	20 m *	100 m *	30 m *	Number of looks of
	Resolution	Azimuth	10 m (2looks) 20 m (4looks)	10 m (2looks) 20 m (4looks)	10 m (2looks) 20 m (4looks)	100 m	10 m (2looks) 20 m (4looks)	the ScanSAR mode is 8or moreby both range and azimuth.
						70 km (3 scan)
	Swath width	70 km *	70 km *	70 km *	300 km (4 scan)	30 km *
						350 km (5 scan)
	Data Rate	240 Mbps	120 Mbps	120/240 Mbps	240 Mbps
	Radiometric Accuracy		Relative accuracy within 1 scene: < 1dB Relative accuracy within 1 orbiting: < 1.

Table 2 .

 2 6. Landsat-TM and SPOT 5-HRG spectral bands and their resolutions

	Sensor	Spectral bands	Wavelength (micrometers)	Resolution (meters)	Acquisition dates
		Band 1 (Blue)	0.45-0.52	30
		Band 2 (Green)	0.52-0.60	30
		Band 3 (Red)	0.63-0.69	30	04/10/2008
	Landsat 5-TM	Band 4 (Near-Infrared)	0.76-0.90	30	08/11/2009 04/06/2010
		Band 5 (Near-Infrared)	1.55-1.75	30	29/12/2010
		Band 6 (Thermal)	10.40-12.50	120* (30)
		Band 7 (Mid-Infrared)	2.08-2.35	30
	SPOT-5 (take 5)	Band 1 (Green) Band 2 (Red)	0.50-0.59 0.61-0.68	10 10	15/04/2015 20/04/2015
	HRG	Band 3 (Near-Infrared) Band (Shortwave Infrared)	0.78-0.89 1.58-1.75	10 10	19/06/2015 24/06/2015
	* TM Band 6 was acquired at 120-meter resolution, but products are resampled to 30-meter
	pixels.

Table 2 . 7 .

 27 Distribution of in situ plots in elevation and slope classes Maximum canopy height was computed using equation: h= a[tan α-tan β]. As it was described in section 2.2.1, five top tall trees were measured in a small circle plot of 36 m diameter and five top tall trees in a co-center 70 m diameter circle plot but outer margin of small plot. So the tallest tree inside and outside of the small plot was called H max -in and H max -out, respectively. The tallest one among ten top trees was also selected and called H max -total.

			Elevation (m)				Slope (%)
		< 200 1000-1500 1500-2000 <10 10-20 20-40 40-50 > 50
	Plots (phase 1)	3	19	38	6	11	19	12	12
	Plots (phase 2)	0	18	14	3	5	16	3	5

To calculate the height of all trees in each plot (in addition to 11 trees that were measured), a variety of non-linear models relating DBH to height, recommended in different studies were selected and tested (Table

2

.8). These relationships were considered for four species as 1) Fagus orientalis, 2) Carpinus betulus, 3) Quercus castanifolia, 4) Alnus subcordata, and two groups of species (similar in shape and height) as Group1 included Tilia begonifolia, Acer velutinum, Acer cappadocicum, Sorbus torminalis and Fraxinus excelsior, and Group2 included Quercus macranthera, Carpinus orientalis, Parotia persica and Diospyros lotus. These six categories have been chosen based on six forest volume tables produced by Forests, Range & Watershed Management Organization (FRWO) for northern forests of Iran. To select the best regression model among a number of models, several most commonly used criteria such as adjusted coefficient of determination (R 2 a ), Root Mean Square Error (RMSE) and Akaike Information Criterion (AIC) were evaluated

Table 2 . 9 .

 29 Six selected non-linear height-DBH models and their statistical performance. H max and DBH stand for maximum height and diameter at breast height, respectively.

	Species	Model	Parameters a b	RMSE (m)	R 2 a
	Fagus orientalis	% SKT = 1.3 4 891 : ; <>.?@A C	41.794	0.025	5.38	0.65
	Carpinus betulus	% SKT = 1.3 4 8. IJ% H 4 IJ% L	33.039	14.772	3.90	0.48
	Quercus castanifolia % SKT = 1.3 4 891 : ; <>.?@A C	39.574	0.035	5.74	0.30
	Alnus subcordata	% SKT = 1.3 4 891 : ; <>.?@A C	39.698	0.038	3.05	0.47
	Group1	% SKT = 1.3 4 IJ% " 98 4 H. IJ%C " L	2.053	0.143	5.24	0.63
	Group2	% SKT = 1.3 4 891 : ; <>.?@A C	13.263	0.025	2.79	0.44

Group1: Tilia begonifolia, Acer velutinum, Acer cappadocicum, Sorbus torminalis and Fraxinus excelsior Group2: Quercus macranthera, Carpinus orientalis, Parotia persica and Diospyros lotus

Table 2

 2

	Species	Model	A	Parameters b	c
	Fagus orientalis Carpinus betulus Quercus castanifolia # = 89IJ% " × %C > # = 98. IJ% " × %C 4 9H. IJ% " C 4 O 0.000026364 0.000342059 0.0405 0.000023 1.0432 # = 89IJ% " × %C > 0.000058 0.9544 Alnus subcordata 0.000051 0.9595 # = 89IJ% " × %C > Group1 0.000056 0.9539 # = 89IJ% " × %C > Group2 # = 89IJ% " × %C > 0.000035 1.0058

.10. Six volume models and their parameters (produced by FRWO)

Table 2 .

 2 11 summarizes statistics for 60 measured plots that were used for developing and validating the GLAS height and volume models. Moreover, distribution of number of trees (n/ha) and basal area (m 2 /ha) is seen in figure2.10. The frequency distribution of the forest volume for 60 GLAS footprints is shown in figure2.11.

Table 2 .

 2 11. Statistical summary of maximum height (H max ), Lorey height (H Lorey ) and forest volume (V) for 60 in situ plots; Min., Max. and Avg. stands for minimum, maximum and average value, respectively.

	Number		H max (m)			H Lorey (m)		V (m 3 /ha)
	of plots	Min. Max. Avg. Min. Max. Avg.	Min. Max. Avg.
	60	5.3	52.3	36	5.1	36.2	27.4	0.69 996.56 457.44

Table 2

 2

			.12. The main information derived from GLA01 and GLA14
	#	Flag	product Description
	1	i_rec_ndx	GLA01 GLA14	GLAS Record Index: GLAS is recording 40 waveforms per second; each 40-waveform bunch has a record index
	2	i_shot_ctr	GLA01 Shot Counter
	3	i_rng_wf	GLA01 The 1064 nm echo waveform digitizer sample output (0-255)

Table 2 .

 2 13. Definition of metrics extracted from GLAS waveforms

			Metrics	Definition
	W ext n , ln(W ext ), exp (W ext )	Waveform extent
	H lead	n , ln(H lead ), exp (H lead )	Height of lead edge extent
	H trail	n , ln(H trail ), exp (H trail )	Height of trail edge extent
	H 25 n , ln(H 25 ), exp (H 25 )	Height at which 25% of the returned energy occurs
	H 50 n , ln(H 50 ), exp (H 50 )	Height at which 50% of the returned energy occurs
	H 75 n , ln(H 75 ), exp (H 75 )	Height at which 75% of the returned energy occurs
	H 100 n , ln(H 100 ), exp (H 100 )	Height at which 100% of the returned energy occurs
	ln: natural logarithm (the logarithm to the base e=2.718), power n=0, 0.5, …, 3

  index) extracted from Landsat7-ETM+ images in Atlantic rainforest fragments, in southeastern Brazil. They resulted that MVI outperformed in dense humid forests, whereas NDVI is a good indicator of green biomass in deciduous and dry forests. They observed a weaker saturation effect and a higher sensitivity to MVI rather than NDVI over dense canopies in the Atlantic rainforest.[START_REF] Pascual | Relationship between lidar-derived forest canopy height and Landsat images[END_REF] have reported high correlation between NDVI and MVI extracted from Landsat-ETM+ and mean and median lidar derived heights (R > 0.6) in pine forests of the Fuenfria Valley in central Spain.[START_REF] Nichol | Improved biomass estimation using the texture parameters of two high-resolution optical sensors[END_REF] estimated forest biomass based on simple ratio vegetation index (RVI) derived from AVNIR-2 and SPOT-5 with R 2 of 0.59 and 0.39, respectively. They observed a significant improvement in biomass estimation with an R 2 of 0.739 obtained from the combined use of RVI of both sensors.

	In this study three vegetation indices including NDVI, MVI and RVI were extracted from
	Landsat-TM and SPOT5 multispectral bands (Equations 2.8 to 2.10).
	zh{| = $ 9z|R : RC 9z|R 4 RC L }{| = $ 9z|R : }|RC 9z|R 4 }|RC L	(2.8)

Table 3 . 1 .

 31 Some MLR models fitted based on waveform metrics

	#		model	reference
	1	H E = aW ext -bTI	Lefsky et al., 2005
	2	H E = aW ext -bTI -c
	3	H E = aLn(W ext ) -bTI -c	Xing et al., 2010
	4	H E = aW ext -b(H lead + H trail )	Chen, 2010b
	5	H E = aW ext -bH lead -cH trail -d	Baghdadi et al., 2014
	6	H E = aW ext -bH trail	Baghdadi et al., 2014
	7	H E = aW ext -bH trail -c	Baghdadi et al., 2014
	8	H E = aW ext -bH lead
	9	H E = aW ext -bH lead -c
	10	H E = aW ext -bTI + cH lead -d
	11	H E = aW ext -bTI -cH trail -d
	12	H E = aW ext -bTI -cH trail	Baghdadi et al., 2014
	13	H 50 + f
	16	H E = aW ext + bTI + cH 50
	17	H E = aW ext + bTI + cH 50 + d
	18	H E = aW ext -bTI + cH 50 -dH 75 + eH 100 + f
	19	H E = a.W ext	(2.5) + b.W ext (1.5) + cTI (1.5) + d.Ln(H 50 ) + e
	20	H E = b.W ext + h.H 75 (1.5) + i.H 75 (2) +b.W ext (1.5) + cTI + dH trail + eH 50 + fH 75 + gH 100 (2.5) + j
	H

E = aW ext -bTI + cH lead -dH trail 14 H E = aW ext -bTI + cH lead -dH trail -e 15 H E = aW ext -bTI + cH lead -dH trail + eH E : Estimated height (maximum or Lorey's height); a,b,c,d,e,f,g,h and i are coefficients; TI: terrain index. Other metrics were presented in table 2,8.

Table 4

 4

.1 represents five models developed based on

Table 4 . 1 .

 41 Statistics for five MLR to estimate H max based on waveform metrics

	#	Model	Coefficients	RMSE (m)	RMSE (%)	R 2	a	MAE MAPE (%)	AIC	P-value
			a= -0.0042
	1	Hmax = a.Wext 2.5 + b.Wext 1.5 + c.ln(H 50 ) + d.TI10 1.5 + e	b= 0.386 d= -0.052 c= 3.549	5.0	13.8	0.85	4.0	16.4	296.3	3.61e-26
			e= 21.22
	2 Hmax = a.Wext + b.TI10	a= 1.1041 b= -0.4910	6.3	17.5	0.76	5.2	23.0	314.2	2.08e-20
	3 Hmax = a.Wext + b.H 50 + c. TI10	a= 0.93538 c= -0.42024 b= 0.36011	5.8	16.1	0.79	4.7	20.3	307.6	2.10e-22
			a= 1.0845
			b= -0.3822
	4	Hmax = a.Wext + b.TI10 + c. H 50 + d.H 75 + e.H 100 + f	c= 0.6642 d= -0.1114	6.1	16.9	0.77	8.3	35.3	312.3	4.00e-21
			e= -0.2605
			f= -0.992
			a= 1.1092
	5 Hmax = a.Wext + b.TI10 + c. Hlead + d.Htrail	b= -0.4948 c= 0.1067	6.4	17.8	0.75	5.3	22.7	316.9	5.67e-20
			d= -0.1319

Fig. 4.2. Estimated H max using MLR based on waveform metrics (model 1, table 4.1) versus in situ H max

  Estimated H max using stepwise regression of waveform metrics (Eq. 4.3) versus in situ H max MLR regressions using all PCs or PCs from stepwise regression (26 PCs) did not produce good result. Three first PCs, explaining 77.5% of data variance, had the highest performance in our models. The statistical result of developing MLR models using three first PCs were presented in table 4.2. The model consisting only three first PCs did not performed well. Adding W ext and TI 10 improved the result considerably. The smallest AIC (301.1) belongs to model combining three first PCs, W ext and TI 10

	+ 0.02297 H 75	1.5 -2.5984 H 75	1.5 -0.48393 TI 90 + 1.06764 H 50 + 10.8052 H 75 + 0.6756 H 100 1.5 + 0.4155 H trail -53.3838 (4.3)
	Fig. 4.5.

Table 4 . 2 .

 42 Statistics for three MLR models to estimate H max based on PCs

	#	Model	coefficient	RMSE (m)	RMSE (%)	R 2 a	MAE (m)	MAPE (%)	AIC	P-value
			a= -5.257
	1	H max = aPC 1 + bPC 2 + cPC 3 + d	b= 16.065 c= 23.479	10.8	30	0.29	9.0	46.4 369.0 3.32e-5
			d= 35.928
			a= 8.1492
	2	H max = aPC 1 + bPC 2 + cPC 3 + dW ext + e	b= 7.7948 d= 1.3427 c= -2.0450	7.7	21.4	0.63	6.4	28.6 337.3 5.57e-11
			e= -21.5505
			a= 4.8863
			b= 5.6782
	3	H max = aPC 1 + bPC 2 + cPC 3 + dW ext + eTI 10 + f	c= -7.1939 d= 1.3460	6.0	16.6	0.77	4.7	22.1 301.1 2.2e-16
			e= -0.4791
			f= -10.5722

10 have the highest and lowest importance in this model. Models 4 and 5 containing edge extent metrics (H lead , H trail ) did not performed as well as other presented models.

Table 4 . 3 .

 43 Statistics of five RF models for estimation of H max based on waveform metrics

					Importance degree
	#	W ext	TI 10	H trail	H lead	H 50	Ln(H 50 )	TI 10 1.5	W ext 2.5	W ext 1.5	RMSE (m)	RMSE (%)	R	2 a	MAE (m)	MAPE (%)
	1	-	-	-	-	-	1973.6	1348.6	2345.2	2402.9	6.8	18.8	0.72	5.4	28.0
	2	5372.6	2596.0	-	-	-	-	-	-	-	7.2	20	0.70	5.7	30.7
	3	3366.4	1716.0	-	-	2928.0	-	-	-	-	7.3	20.3	0.69	5.6	30.5
	4	2336.6	1270.8	885.7	1308.9	2222.2	-	-	-	-	7.7	21.4	0.70	6.3	35.2
	5	3270.0	1830.0	1069.5	1650.6	-	-	-	-	-	8.5	23.6	0.62	6.7	38.8

Table 4 .

 4 4. Statistics of three RF models for estimation of H max based on PCs

			Importance degree
	#	W ext	TI 10	PC 1	PC 2	PC 3	RMSE (m)	RMSE (%)	R 2 a	MAE (m)	MAPE (%)
	1	-	-	1807.3	1914.4	3839.5	9.3	25.8	0.50	7.4	40.6
	2	2920.2	-	1081.0	1338.3	2427.0	8.0	22.2	0.63	6.1	33.1
	3	2359.4	1113.7	928.5	1131.8	2153.9	8.0	22.2	0.66	6.1	36.0

Table 4 .

 4 5. Properties of two ANN models for estimation of H max based on waveform metrics and the resulted statistics

			Properties of network
	#	Input	Number of	Number of	Iteration	RMSE (m)	RMSE (%)	R 2 a	MAE (m)	MAPE (%)
			hidden	hidden	rate
			layers	neurons
	1	W ext , TI 10	1	2	30	5.7	15.8	0.80	4.6	20.1
	2	W ext , TI 10 , H 50	1	2	20	5.4	15	0.82	4.3	17.1

Table 4 .6.

 4 Properties of ANN models for estimation of H max based on PCs and the resulted statistics Fig. 4.11. Estimated H max using ANN based on PCs (model 3, table 4.6) versus in situ H max

			Properties of network
	#	Input	Number of hidden layers	Number of neurons hidden	Iteration rate	RMSE (m)	RMSE (%)	R 2 a	MAE (m)	MAPE (%)
	1 PC 1 , PC 2 , PC 3	1	3	20	8.8	24.4	0.53	6.6	31.3
	2	PC 1 , PC 2 , PC 3 , W ext	1	4	15	6.6	18.3	0.73	5.0	20.3
	3	PC 1 , PC 2 , PC 3 , W ext , TI 10	1	4	15	5.6	15.5	0.81	4.3	17.4

Table 4 .

 4 7. Statistics for five MLR for estimation of H Lorey based on waveform metrics

	#	Model	Coefficients RMSE (m)	RMSE (%)	R 2 a	MAE MAPE (%)	AIC	p-value
	1	H Lorey = a.ln(W ext ) + b.TI 10 +c	a= 27.6671 c= -67.8802 b= -0.3454	5.1	18.6	0.71	3.9	24.0	288.3 8.18e-18
	2 H Lorey = a.W ext + b.TI 10	a= 0.8079 b= -0.3252	5.4	19.7	0.70	4.1	23.0	293.1 2.19e-16
	3	H Lorey =a.W ext b.Ln(W ext ) + c.TI 10 2.5 +	a= 0.0007 c= -0.3651 b= 7.1533	5.7	20.8	0.67	4.4	28.7	300.8	2.2e-16
	4	H Lorey = a.W ext + b.H 50 + c.TI 10	a= 0.7671 c= -0.305 b= 0.0876	5.4	19.7	0.67	4.3	23.5	294.4 1.85e-16
			a= 0.7475
	5	H Lorey = a.W ext + b.TI 10 + c.H lead + d.H trail + f	b= -0.3408 d= 0.0369 c= -0.0503	5.6	20.4	0.63	4.5	25.9	297.7 2.61e-15
			f= 3.3432

Fig. 4.13. Estimated H Lorey using MLR based on waveform metrics (model 1,

Table 4.7) versus in situ H Lorey

Table 4 .8.

 4 Statistics of three MLR models for estimation of H Lorey based on PCs

	#	Model	Coefficient	RMSE (m)	RMSE (%)	R 2 a	MAE (m)	MAPE (%)	AIC	P-value
			a= 3.856
	1	H Lorey = aPC 1 + bPC 2 + cPC 3 + d	b= 12.192 c= 20.590	7.8	28.4	0.30	6.4	43.9	340.3 1.68e-5
			d= 26.939
			a= -7.7682
	2	H Lorey = aPC 1 + bPC 2 + cPC 3 + dW ext + e	b= 2.2925 d= 0.8996 c= 7.6185	6.6	24.1	0.51	5.2	31.6	320.0 4.19e-9
			e= -11.5824
			a= -4.5647
	3	H Lorey = aPC 1 + bPC 2 + cPC 3 + dW ext + eTI 10 + f	b= -0.2831 e= -0.3080 c= 1.5421 d= 0.9006	5.4	19.7	0.66	4.0	24.1	304.9 1.42e-11
			f= -4.4943

Table 4 . 9 .

 49 Statistics of five RF models for estimation of H Lorey based on waveform metrics

					Importance degree
	#	W ext	TI 10	H trail	H lead	H 50	Ln(H 50 )	TI 10 1.5	W ext 2.5	W ext 1.5	Ln(W ext )	RMSE (m)	RMSE (%)	R	2 a	MAE (m)	MAPE (%)
	1	-	-	-	-	-	740.2	490.0	1018.5	1056.8	-	5.4	19.7	0.68	3.9	26.9
	2	-	599.8	-	-	-	-	-	1308.3	-	1378.7	5.5	20	0.67	3.9	26.1
	3	1527.4	604.2	-	-	1088.4	-	-	-	-	-	5.6	20.4	0.66	4.2	29.8
	4	2317.0	861.1	-	-	-	-	-	-	-	-	5.7	20.8	0.66	4.2	28.5
	5	1352.7	626.9	437.7	731.8	-	-	-	-	-	-	6.6	24.1	0.56	5.1	37.9

Table 4 .

 4 10. Statistics of three RF models for estimation of H Lorey based on PCs

			Importance degree
	#	W ext	TI	PC 1	PC 2	PC 3	RMSE (m)	RMSE (%)	2 a R	MAE (m)	MAPE (%)
	1	-	-	1003.7	1120.0	2178.2	6.4	23.3	0.59	4.8	35.3
	2	1448.9	-	681.4	766.0	1366.2	5.9	21.5	0.65	4.3	31.2
	3	1276.0	686.7	563.3	615.0	1133.5	5.6	20.4	0.72	4.3	31.3

Fig. 4.15. Estimated H Lorey using RF based on PCs versus in situ H Lorey (a, b, c: model 1 to 3, table 4.10)

  Three simple ANN models were presented in table 4.11. H Lorey resulted from model 1 with two metrics (W ext and TI 10 ) was compared with in situ height in figure 4.16. Based on MAPE statistic, 23.2 percentage of this model predictions are off.

		Table 4.11. Statistics of three ANN models for estimation of H Lorey based on waveform metrics
			Properties of network
	#	Input	Number of hidden layers	Number of hidden neurons	Iteration rate	RMSE (m)	RMSE (%)	2 a R	MAE (m)	MAPE (%)
	1	W ext , TI 10	1	2	10	5.1	18.6	0.72	3.7	23.2
	2 W ext 2.5 , Ln(W ext ), TI 10	1	3	10	5.0	18.2	0.72	3.5	20.9
	3	W ext , TI, H 50	1	3	10	5.2	18.9	0.69	4.0	23.5
		Fig. 4.16. Estimated H Lorey using ANN based on waveform metrics (model 1, table 4.11) versus in
				situ H Lorey

Table 4 .

 4 12. Statistics of three ANN models for estimation of H Lorey based on PCs

			Properties of network
	#	Input	Number of hidden layers	Number of hidden neurons	Iteration rate	RMSE (m)	RMSE (%)	2 a R	MAE (m)	MAPE (%)
	1	PC 1 , PC 2 , PC 3	1	3	35	4.7	17	0.76	3.5	17.4
	2	PC 1 , PC 2 , PC 3 , W ext	1	4	50	3.4	12.4	0.87	2.5	12.3
	3	PC 1 , PC 2 , PC 3 , W ext , TI 10	1	4	30	3.6	13	0.85	2.6	15.5

Table 4 .

 4 13. Comparison of Lorey's height, number of trees and volume in three couple of plots Concerning statistical methods, the result of an ANN based on PCs is approximately the same as result obtained by MLR model based on waveform metrics. But positive points about ANN model are: 1) It does not need an ancillary data (DEM); 2) It relies on PCs and W ext which are not prone to uncertainties in contrast with some waveform metrics which have been used in MLR model.Generally, almost the same accuracy as volume-H Lorey method was obtained using the ANN model. ANN performed slightly better where there exist very low ( <10 m 3 /ha) and very high ( > 800 m 3 /ha) volume. In most other points ANN produced more error rather As it is observed, there is no obvious trend in volume estimation by increasing in forest density or terrain slope. This confirms that the heterogeneity of forest reduces the ability of lidar data to estimate forest volume. A collection of forest properties including forest type, horizontal and vertical structure of forest, and topographical properties may impress predictions. Consideration of reference plots conditions demonstrated dependency of volume on diverse factors. For instance, it happened to have high number of trees per hectare but low volume and reverse, and also equal number of trees per hectare or equal mean height but different volume(Table 4.13). It is expected to have higher accuracies in homogenous forests.

	Plot	Lorey's height (m)	Number of trees (n/ha)	Volume (m 3 /ha)
	1	16.7	39	5.2
	3	17.7	697	160.5
	13	32.9	484	996.6
	17	32.8	322	499.2
	23	6.5	8	0.7
	25	7.5	237	20.7

 [START_REF] Yang | The potential of Weibull-type functions as flexible growth curves[END_REF][START_REF] Zhang | Cross-validation of non-linear growth functions for modelling tree height-diameter relationships[END_REF]Peng et al., 2001% = 1.3 4 891 : ; <>.?@A C D [START_REF] Batista | Performance of height-diameter relationship models: analysis in three forest types[END_REF][START_REF] Zhang | Cross-validation of non-linear growth functions for modelling tree height-diameter relationships[END_REF] % = 1.3 4 8$. ; E > ?@AFD G [START_REF] Huang | Comparison of nonlinear height-diameter functions for major Alberta tree species[END_REF][START_REF] Peng | Nonlinear height-diameter models for nine boreal forest tree species in Ontario. Ministry of natural resources[END_REF][START_REF] Ahmadi | Non-linear height-diameter models for oriental beech (Fagus orientalis Lipsky) in the Hyrcanian forests, Iran[END_REF] et al., 1992;[START_REF] Peng | Nonlinear height-diameter models for nine boreal forest tree species in Ontario. Ministry of natural resources[END_REF] adjacent to the trees in a moist tropical forest. They measured trees with five technicians and obtained an RMSE on the tree heights of minimum 2.88 m. They concluded that these methods produced unbiased height estimates but also high level of random error. In our research, we minimized these sources of errors by walking around the tree and viewing it from different angles to distinguish the actual top from other branches. For some trees, measurements were done from different angles by two measurers to compare and justify the measurements. After all, it is expected to have an RMSE of about 2 to 3m on in situ tree height. The DBH was also measured using a caliper with 0.5cm precision.

-Concerning Lorey's height and volume, it was needed to obtain the height of all trees located in each plot. As explained in section 2.2.1, all trees were measured in case of DBH, but only eleven trees in case of height. So, height-DBH relationships were used to address this requirement (Table 2.9, section 2.3.1) which will cause some bias from the true heights. Field volume was also calculated using local volume allometric models provided by FRWO (Table 2.10, section 2.3.1) which may contain uncertainties. Following paragraphs give an overview on the uncertainties associated to the field volume calculation:

As shown in Table 2.9, the accuracy of estimating tree height from DBH is between 3 and 5.7 m. Since the RMSE of volume allometric models has not been reported by FRWO, the precision on the estimation of tree volume (V) using allometric model of Carpinus betulus, for instance, is calculated:

Where, V is in m 3 , H, and DBH are in m and cm, respectively. The relationship between the precision on the estimation of volume and the precision on height and DBH can be written as equation 5.1:

Where ∆V/V is the relative precision on the estimation of volume, ∆DBH/DBH and ∆H/H are the relative error on DBH and height, respectively. The coefficient b is equal to 1.0432 (Table 2.10). For a DBH of 7.5 and 124cm, as minimum and maximum DBH of Carpinus betulus, respectively, the estimated height would be 12.4m and 30.8m (refer to height relationship for Carpinus betulus, table 2.9). According to a DBH accuracy (ΔDBH) of 0.5cm and height accuracy (ΔH) of 3.9m (Table . 2.9), the maximum and minimum relative error on estimation volume would be 46.7% (for low H and DBH) and 14% (for high H and DBH), respectively.