Waveform LiDAR concepts and applications for potential vegetation phenology monitoring and modeling: a comprehensive review

ABSTRACT Researchers continually demonstrated through published literature how LiDAR could create unparalleled measurements of ecosystem structure and forest height. There are a number of studies conducted utilizing waveform LiDAR products for terrestrial monitoring, but those that deal specifically with the assessment of space-borne waveform LiDAR for monitoring and modeling of phenology is very limited. This review highlights the waveform LiDAR system and looks into satellite sensors that could link waveform LiDAR and vegetation phenology, such as the proposed NASA’s Global Ecosystem Dynamics Investigation (GEDI) and the Japanese Experimental Module (JEM)-borne LiDAR sensor named MOLI (Multi-footprint Observation LIDAR and Imager). Further, this work examines the richness and utility of the waveform returns and proposes a spline-function-derived model that could be exploited for estimating the leaf-shooting date. The new approach may be utilized for ecosystem-level phenological studies.


Introduction
Phenology, the study of the timing of recurrent biological events (Lieth 1974) driven by environmental change is perhaps the simplest process in which to track changes in the ecology of species in response to climate change (Intergovernmental Panel on Climate Change (IPCC) 2007). Globally, according to the Intergovernmental Panel on Climate Change (IPCC), most of our knowledge about changes in terrestrial phenology comes from the northern hemisphere, with majority of the studies coming from Europe alone. Although no specific numbers were shown regarding the spread of the individual phenological areas of study, Google scholar has a rough estimation that shows phenological studies conducted on forest biomes and plant species constituting the major bulk. Using various keywords such as plant phenology, forest phenology, animal phenology, marine phenology, and other word combinations such as phenological studiesforest and vegetation phenology always appears at the top pages. Table 1 shows the estimated numbers based on the first three result pages.
The popularity of plant phenology is partly attributed for being the most responsive and easily observable traits in nature, and most importantly for its implications on the ecosystem-atmosphere interactions. Changes of the phenological events, such as the emergence and senescence of leaf, are not only critical to survival and reproduction (Rathcke and Lacey 1985) but also indicate important climatic variations (Schwartz and Hanes 2009;Peñuelas et al. 2013). In the last two decades, the attention on phenology has shifted to its function as a proxy in detecting global climate change (Jackson et al. 2001;White, Brunsell, and Schwartz 2003;Parmesan 2006;Inouye 2008;Miller-Rushing and Primack 2008;Peñuelas, Rutishauser, and Filella 2009;Morin et al. 2010;Cook, Wolkovich, and Parmesan 2012a;Yang et al. 2017) as modeled growing season length affects tree competition and vegetation dynamics (Kramer, Leinonen, and Loustau 2000) and correlates with terrestrial CO 2 uptake (Piao et al. 2017). Further, more and more studies use phenological information to make inferences on the relationships of life cycle events of vegetation to their environment (Reed et al. 1994;Fitzjarrald, Acevedo, and Moore 2001;Schmid et al. 2003;Gordo and Sanz 2010;Wilsey et al. 2018).
Phenological indices that had been assessed in Europe (Menzel and Fabian 1999;Menzel 2000;Chmielewski and Rotzer 2001;Ma et al. 2016) had shown earlier springtime phenological events that could contribute to an increase biomass formation, which is part of a global increase in biospheric activity. In western North America, the same advancement in spring timing was observed (Cayan et al. 2001;Schwartz, Ahas, and Aasa 2006). In order to monitor the implications of the changing plant phenology on biosphere-atmosphere interactions, which is highly variable, phenological models have to be in place. For decades, phenological models have been developed to predict the timing of leaf onset, such as the thermal time model (Cannell and Smith 1983), alternating model (Murray, Cannel, and Smith 1989), sequential model (Hänninen 1990), parallel model (Kramer 1994), and alternating model (Murray, Cannel, and Smith 1989). The start of the leafing marks the start of the photosynthetic season for a deciduous forest (Delpierre et al. 2017) and is a major determinant of forest duration (Donnelly, Yu, and Caffarra 2017) and an indication of the onset of the net carbon uptake (White, Thornton, and Running 1997). Additionally, Richardson et al. (2009) showed that the day of leaf emergence controls the annual carbon sequestration in forests. Carbon uptake refers to the net annual per-hectare removal of carbon from the atmosphere while the land is in a given land-use category. The length of the carbon uptake period across a latitudinal and continental gradient of deciduous forests explained 80% of the spatial variance in annual Net Carbon Exchange Of Ecosystems (NCEE) (Baldocchi et al. 2001). The importance of the timing of the forest leaf onset to the yearly variability of the NCEE had already been described and conveyed both theoretically (White, Running, and Thornton 1999) and experimentally (Black et al. 2000;Schmid et al. 2003).
Within the last couple of decades, the importance of long-term forest phenological data has been re-emphasized by examining existing data sets for evidence of species change and for temperature response (Menzel and Estrella 2001;Melaas, Friedl, and Richardson 2016). The low-tech approach offered by phenology has now being linked with the high-tech approach through the comparison of satellite observations with those on the ground (e.g. Duchemin, Goubier, and Courrier 1998;Guyon and Lagouarde 1999). In fact, more and more studies are now making use of the remote sensing data to look at forest leaf-out, leaf-fall (Zhang et al. 2004;Robin et al. 2007;Liang and Schwartz 2009;White, Pontius, and Schaberg 2014).
Remote sensing satellites could measure forest features such as leaf phenology that are difficult if not impossible to collect with ground crews (Ustin et al. 2004;DeFries 2008). Although ground-collected forest measurements are often more accurate than satellite data at the point and time of collection, satellites gather data at wide extent areas, sample the full range of variation in forest metrics, and capture broad trends and dynamic change in the world's forests (Houghton 2005). Hence, satellite datasets allow for the integration of ground measurements and monitoring of vegetation dynamics at regional to global scales (Myneni et al. 1997).
Several methods exist to estimate the Start Of Season (SOS) (see de Beurs and Henebry 2010 for the list) or "distinct green coloration in the forest" (Thoman and Fathauer 1998) from satellite data. The identification of the SOS in terms of the Day Of Year (DOY) is the key to the seasonal characterization of vegetation. There is, however, no standard approach to determine the SOS (Reed, White, and Brown 2003;White et al. 2009;Fu et al. 2014). One common method is the threshold that uses sets of Normalized Difference Vegetation Index (NDVI) or a value calculated from the minimum and maximum NDVI to know the start, peak, end, and length of the growing season (Lloyd 1990;Suzuki, Nomaki, and Yasunari 2003;Shang et al. 2018). NDVI is a quantity that measures greenness and vigor of vegetation (Tarpley, Schneider, and Money 1984). Related to the threshold method is the compound median filter approach that also uses NDVI for monitoring phenological changes (Kogan 1995). In the inflection point approach, the time of transition from the temporal NDVI profile is detected and metrics are derived with time derivatives or logistic functions (Badhwar 1984;Moulin et al. 1997;Zhang et al. 2003Zhang et al. , 2004. The curve-derived method uses delayed moving average (Reed et al. 1994) and time of largest NDVI increase (Kaduk and Heimann 1996). However, these methods are difficult to apply at global scales, and generally do not account for ecosystems characterized by multiple growth cycles (Zhang et al. 2003;Yu et al. 2017). The modeling approach by de Henebry (2004, 2005) models land surface phenology (LSP) using NDVI as a function of accumulated growing degree day (AGDD). Their phenological model used a simple regression model that assessed whether institutional changes in Kazakhstan had an effect on land surface phenology. The results showed that across different eco-regions, the model explained a significant portion of NDVI variation and could be used to assess significant changes in land surface phenology.
Passive remote sensing has been effective in monitoring cycles of vegetation phenology. For instance, the Advanced Very High-Resolution Radiometer (AVHRR) has been extensively used for monitoring vegetation phenology (Lloyd 1990;Reed et al. 1994;White, Thornton, and Running 1997;Delbart et al. 2006;Zhang, Liu, and Dong 2017). AVHRR offers the only source of global data that could be used to analyze seasonal-to-decadal-scale dynamics in global vegetation (e.g. Zhou et al. 2001). However, the lack of precise calibration, poor geometric registration, and difficulties involved in cloud screening AVHRR data could result in high levels of noise (Goward et al. 1991). Then, there is the Moderate-Resolution Imaging Spectroradiometer (MODIS) instrument that provides daily reflectance data and possesses seven spectral bands that are specifically designed for land applications with spatial resolutions as small as 250 m ( Table 2). The improved spectral, radiometric, and geometric qualities of MODIS data provide an effective means of monitoring global environmental changes. The products generated by MODIS on board Terra and Aqua present an unprecedented chance for scientists to develop long-term records of vegetation phenology. Time-series MODIS data have been used for phenology detection (Zhang et al. 2003;Sakamoto et al. 2005;Peng et al. 2017;Parente and Ferreira 2018). Phenological models have been studied to predict the onset of leaf using MODIS data (Kang et al. 2003;Kim and Wang 2005;Ahl et al. 2006). Zhang et al. (2003) developed a method to estimate phenological events based on the curvature-change rate for MODIS time series. However, confusion on the calculation of specific dates of vegetation phenological transitions still exist due to the mix of methods and definitions being used (White and Nemani 2003). Also, small-scale (<500 m) topographical variability in the order of 50 m could result in a 1 to 2-week difference in SOS (Fisher, Mustard, and Vadeboncoeur 2006). It is also imperative to establish whether temporal NDVI derived from the satellite sensors can uncover the SOS dates since the spatial extent of satellite data permits for a larger area of study than field observations. Table 2 provides a list of satellite sensors that may be available for phenological studies of vegetation. The primary advantage of coarse-resolution satellite datasets is their long-term and frequent revisit time. Medium (e.g., Landsat) to high (e.g., IKONOS) spatial resolution sensors are limited by temporal resolution that could render them insufficient in mapping phenological dynamics in vegetation (Gonzalez-Sanpedro et al. 2008;Verbesselt et al. 2010). In addition, due to cloud cover and archive restrictions, Landsat availability outside the United States is much lower (Wulder et al. 2016). There have been studies that integrated higher temporal-resolution images with Landsat using image fusion techniques (Gao et al. 2006;Roy et al. 2008;Zhu et al. 2010) and showed comparable results with those using MODIS dataset (Bhandari, Phinn, and Gill 2012;Li et al. 2019).
Sentinel-2 has a richer spatial and spectral content compared to Landsat, and with a higher revisit time of 5 days (Table 2). However, Sentinel-2 acquisitions are available only from March 2017 onward. There are technological tradeoffs between the spatial resolution of passive sensors and the dimension of the swath they could capture (Rosenqvist et al. 2003;Andersen 2009). The set of commercial high-resolution satellites like the IKONOS that has a spatial resolution of 4 m and a swath width of 13.8 km is in one extreme. In the other extreme is the set of coarseresolution sensors like AVHRR that has a resolution of 1 km and a swath width of 3000 km. For satellites, this tradeoff in spatial resolution and image size affects the time it takes for a satellite to image the same location on Earth (temporal resolution). Coarse-resolution satellites could sense the entire Earth fast and revisit the same location every few days in contrast with high-resolution satellites, while moderate-resolution satellites could take 10 to 20 days to revisit a location. Temporal resolution is imperative in monitoring the timing of phenology. Detecting the precise date of a phenological event depends upon the sampling frequency or with the frequency of passage and pixel size associated with remote sensing instruments mounted on a satellite or tower.
Active sensors (e.g., light detection and ranging known as LiDAR) only constitute a small proportion of the present satellite fleet, but based on future scheduled launching; it would dramatically increase in number and complexity in the next few years. This would make new types of phenological analysis possible. The first LiDAR satellite, the Ice, Cloud, and land Elevation Satellite (ICESat) that ended in the fall  (Dubayah et al. 2014). The products from these active sensors would revolutionize mapping and the monitoring of forest biomass. In the case of GEDI, it has been designed to create unparalleled measurements of ecosystem structure and forest height using LiDAR waveforms (Qi and Dubayah 2016). These satellites are elaborated further in a later chapter of this paper.
LiDAR remote sensing promises improved accuracy of biophysical measurements and extends the spatial analysis to the third dimension. However, the algorithms, models, feature detection for LiDAR analysis are still to mature as data availability grows (Zhan et al. 2011). Further, there has been limited use of the technology for phenological monitoring and modeling. LiDAR research has focused on the measurement of canopy height and structure (Dubayah et al. 1997;Drake et al. 2002a;Khosravipour et al. 2014;Popescu et al. 2018;Silva et al. 2018) and estimation of aboveground biomass (Drake et al. 2002a; even at a plot level: Zhao, Popescu, and Nelson 2009;Powell et al. 2010;Clark et al. 2011;Rosette et al. 2012;Zolkos, Goetz, and Dubayah 2013;Marvin et al. 2014), with results illustrating unprecedented accuracy and consistency (Hurtt et al. 2004), even in complex canopies (Drake et al. 2002b;Salas and Henebry 2016). With the current technology, it remains a challenge to distinguish primary forests from older secondary forests in remote sensing image. The additional highresolution LiDAR systems available for fusions (e.g. Asner et al. 2008;Swatantran et al. 2011) show that there is a promise to address the challenge. Nevertheless, research seems to be stuck at using the active system on case studies that focus on a static condition at a single point in time. The question still remains whether or not LiDAR data could benefit the research arena that focuses on phenological monitoring and modeling.
The objective of this paper is to (1) present the requirements for phenological monitoring and modeling and (2) look into the possibility of waveform LiDAR, specifically space-borne systems, to address those requirements. Further, this paper is broken down into three main sections: (1) an introductory section that provides the background of phenology (2) a discussion section on LiDAR and its capabilities, including LiDAR instruments such as the proposed GEDI and MOLI (3) the later part covers the ways full-waveform LiDAR could be utilized for phenological monitoring and modeling, which contains an overview of the potentials of LiDAR that may have remained untapped.

Terrestrial applications of LiDAR
LiDAR is active remote sensing with main components such as the laser scanner, GPS receiver, the Inertial Monitoring Unit (IMU) and the control unit (Reutebuch, Andersen, and McGaughey 2005). The concept of the technology is to send laser pulses toward a target to obtain high-resolution measurements of surface altitudes, either from airborne or space platforms (Parker, Harding, and Berger 2004;Vierling et al. 2008). The information collected by the laser scanner is geo-referenced with position and orientation information from the GPS and IMU.
LiDAR system measures the elapsed time and amount of return energy pulses scattered back from the target to the receiving sensor. The time elapse (t) of the pulse traveling a distance (d) at the speed of light (c) is described by equation 1: where the speed of light is assumed constant. The time component is essential for canopy height estimation. Taking the difference of the time recording at the top of the canopy (t 1 ) and the time recording at the underlying ground (t 2 ) may result in an estimate of the tree height (h) as shown in equation 2.
Typical commercial LiDAR instruments operate at around the 1.064 µm wavelength to insure high atmospheric transmittance and minimal loss of signal as the pulse travel through the atmosphere.
LiDAR returns could be categorized as either discrete or waveform. The two variations differ in terms of how each vertically and horizontally sample a threedimensional canopy structure. The vertical sampling is described in relation to the number of range samples recorded for each emitted later pulse, while horizontal sampling is the area of the footprint coverage and the number of footprints per unit area ).

Discrete return LiDAR
Early discrete return LiDAR systems typically recorded only one discrete return, either the first peak or the final peak (Naesset 1997;Magnussen, Eggermont, and LaRiccia 1999) or two, the first peak and the final peak (Lefsky et al. 1999a) in the reflected wave ( Figure 1). The first generation LiDAR was used to obtain surface elevations from airborne systems (Krabill et al. 1984) and usually covered diameter footprints as small as 1 m, 5 cm, to 30 cm (Arp, Griesbach, and Burns 1982). Due to the relatively small geographic area covered by the sensor, challenges in the analysis of the datasets, and the lack of a standardized method for the geolocation have limited the use of conventional LiDAR sensors (Lefsky et al. 1999a). Despite the intense research efforts, practical applications of small-footprint LiDAR have not progressed as far, mainly because of the current cost of LiDAR data (Farid, Goodrich, and Sorooshian 2006).
Mapping large areas using small beam size requires expensive and extensive flying , although small-footprint LiDAR has shown its capability for detailed mapping of the ground and mapping snags and understory shrubs (Martinuzzi et al. 2009). Additionally, tree height estimation using small-footprint and limitedreturn LiDAR systems could pose major questions in the estimation of the "real" canopy heights. Having only one or two reflective returns cannot ascertain a measurement of a shot that penetrates all the way to the ground . This, however, changed when in 2000 commercial systems introduced multiple (3 to 5) returns per pulse ( Figure 1). Still, with discrete return systems, the full vertical structure of the vegetation is not sensed and only a portion of the canopy and ground is recorded . There is also a possibility of missing the tree tops (Zimble et al. 2003) and down-bias the canopy heights. These deficiencies could cause inherent information about the reflecting object and its geometric and reflection characteristics totally neglected. With fewer guarantees that the top of the canopy would have a return record from discrete return systems, it is typical to apply calibration to LiDAR measured canopy heights (ML, DB, and DA 2004). Another approach to ease the limitation is to consider only the upper 10% or so of the canopy height (Popescu and Wynne 2004). The method employed by Harding et al. (2001) used the relative Canopy Height Profile (CHP), which could create plot-level CHPs derived by summing all individual CHP within a study plot.

Waveform LiDAR
The waveform LiDAR system, in contrast, records many returns per emitted pulse within the vertical structure of forest canopies (Figure 1), hence, showing details of the intercepted surfaces or the proportion of the canopy complexity. Dubayah et al. (1997) and Hurtt et al. (2004) demonstrated how the LiDAR data could accurately retrieve canopy heights.  converted the elapsed time difference between the peaks of the two most prominent modes in the amplitude waveform into range to estimate the height. However, the highest peak of the largest mode of the waveform does not necessarily correspond to the highest point of the canopy, but is instead a function of the canopy structure. For a complex forest system, it is likely to have a local variation in tree heights. Salas and Henebry (2016) introduced a novel method by exploring the asymmetrical shape of the waveform and its key profile landmarks (e.g., location and magnitude of peaks) as a canopy height indicator. By using synthetic and real LiDAR datasets, Salas and Henebry (2016) found that the method could work for any type of waveform, especially for a single-peak type where canopy height (canopy peak location minus ground peak location) cannot be estimated due to the absence of one peak. The method proved to be robust for any type of waveform noise (Salas, Amatya, and Henebry 2017).
The waveform LiDAR recorded returns could also give an ample amount of data that could indicate the type of canopy structure and describe the vertical canopy volume distribution (Lefsky et al. 1999a). As stands age or grow, the vertical distribution of canopy elements changes (Lefsky et al. 1999b). The waveform could exhibit these changes relative to younger stands (Salas and Henebry 2016) and successional canopy stages (Gu, Sen, and Sanchez-Azofeifa 2018). For instance, bimodal distributions are associated with the presence of an understory that may occur in more mature stands. Older stands characterized by canopy gaps and trees of multiple ages and sizes exhibit a more even distribution of canopy components.
Full-waveform reflections are offered by large-footprint LiDAR (>10 m in diameter). Using large footprints reduces the cost of mapping large forest regions (Blair, Rabine, and Hofton 1999), rapidly records vertical canopy profiles (Drake et al. 2002a), and provides a more defined vertical arrangement of forest structure from canopy top to the ground surface (Dubayah et al. 1997;. For research, where the vertical profile is a critical measure, it is best to utilize fullwaveform LiDAR when laser energy is densely sampled as opposed to a discrete system . The general diagram in Figure 1 shows the concept of the fullwaveform LiDARfrom how the vertical structure (at the scale of the individual trees) could be an important component to calculating the ground and canopy height.
Waveform LiDAR characterizations of ground and vegetation profiles are consistently accurate over the past years of research showing its huge potentials for broad-scale applications (Lefsky et al. 1999a;Lefsky et al. 2002a;Drake et al. 2002a, Zhao, Popescu, andPopescu et al. 2011;Gwenzi and Lefsky 2014). Plant biomass is approximately 50% carbon (Drake et al. 2003;Gibbs et al. 2007), and estimates of the total aboveground biomass in forest ecosystems are critical for carbon dynamics studies (Asner et al. 2013;Laurin et al. 2014;Réjou-Méchain et al. 2015) at multiple scales and periods (Simonson et al. 2016). A consistent relationship between forest canopy structure from LiDAR and aboveground biomass has been illustrated in three different biomestemperate deciduous, temperate coniferous, and boreal coniferous biomes (Lefsky et al. 2002) using a variety of techniques including height quantile relationships and simple regression (Lim and Treitz 2004;Lefsky et al. 2005;Nie et al. 2017). More and more results from studies could have significant implications on how global observations from spaceborne LiDAR instruments should be used to produce global estimates of terrestrial aboveground biomass (Drake et al. 2003;Chi et al. 2015;Margolis et al. 2015;Hu et al. 2016;Wang, Sun, and Xiao 2018). Results showed that the generality of the relationships between LiDAR metrics and aboveground biomass in closedcanopy Neotropical forests were affected by forest dynamics such as leaf phenology.
Waveform LiDAR dataset could be utilized to estimate the vertical distribution of light transmittance ) and vertical foliar profiles. There was a good agreement (r 2 > 0.70) of the Scanning LiDAR Imager of Canopies by Echo Recovery (SLICER) dataset against the upper and lower distribution of canopy components ). In addition, these types of datasets are important to NASA's Earth Science Enterprise (ESE) Ecology/ Carbon cycle program (Blair, Hofton, and Luthcke 2001). Another waveform LiDAR application is on the estimation of LAI and canopy cover Nie et al. 2016). There have been efforts to measure the fraction of canopy cover using space-borne LiDAR by comparing ground reflectance to canopy reflectance (Lefsky et al. 2005). These efforts require canopy gaps to be large enough to encompass the footprint size. Ground or understory returns indicate canopy openness and LiDAR waveforms could provide this information. However, accurate retrieval of canopy coverage would also require retrieval of the reflectance from other materials involved, as wave returns record not only intercepted leaves but also branches. Leaf inclination angle and topography must be determined precisely as well (Zhao et al. 2013). Although more has to be done in this aspect of LiDAR remote sensing, the present approaches of canopy cover estimation would suffice in the absence of much better methods.
There are developments in inverting models against LiDAR observations (Ranson and Sun 2000). A 3D LiDAR waveform model was inverted to estimate tree height, fractional cover, and overstory LAI of a coniferous forest (Koetz et al. 2006). Using two separate datasets, the approach successfully demonstrated the potential of model inversion to retrieve horizontal and vertical forest structure from LiDAR data. However, processing difficulties limited the accuracy of the obtained results. A follow-up study by Koetz et al. (2007) combined two different sensor types, LiDAR and imaging spectrometers, and successfully derived a comprehensive canopy characterization relevant for the assessment of biomass and productivity of vegetation. However, the method needs to be tested on actual measurements of the respective sensors and validated against a large variability of real field measurements.
Whereas LiDAR provides detailed forest structure information, spectral remote sensing is more sensitive to vegetation composition and phenology. The fusion of the capabilities of different digital image data technologies is a valuable tool in remote sensing analysis (Pohl and Van Genderen 1998) and could stimulate more scientific studies on LiDAR. It is through fusion that canopy cover, LAI, tree density, etc., could be best recovered and provide a much more comprehensive view for understanding the ecosystems. For instance, the use of hyperspectral and LiDAR datain tandemto create maps predicting species abundance patterns (derived primarily from Airborne Visible and Infrared Imaging Spectrometer, AVIRIS, data) augmented with coincident patterns of stem size or height (derived primarily from LVIS data) provided a useful adjunct to traditional canopy inventory approaches (Anderson et al. 2008). Geological outcrop has been analyzed by fusing terrestrial LiDAR and hyperspectral products (Buckley et al. 2013;CChu et al. 2016). Forest biomass estimates have been improved by integrating LiDAR and hyperspectral image data (Lucas, Lee, and Bunting 2008;Swatantran et al. 2011;Manzanera et al. 2016). Also, the system was used together with AVIRIS to map the three-dimensional spectral and structural properties of Hawaiian forests (Asner et al. 2008) with results highlighting the location and fractional abundance of each invasive tree species throughout the rainforest sites. Cao et al. (2018) modeled forest aboveground biomass in arid and semi-arid regions of China using waveform LiDAR and the Chinese satellite ZiYuan-3 (ZY-3). The authors found that through the fusion of two datasets, AGB estimates were optimized.
Aboveground biomass was accurately estimated across the entire forested region of southern Quebec province in Canada by integrating DEM information from the Shuttle Radar Topography Mission (SRTM), spaceborne waveform LiDAR from the ICESat Geoscience Laser Altimeter System (GLAS), airborne profiling LiDAR collected over ground plots, and ground inventory plot data (Boudreau et al. 2008). Also, GLAS data have been combined with LandSat data to create accurate regional biomass and height estimates (Helmer, Lefsky, and Roberts 2009;Wang, Sun, and Xiao 2018) and MODIS to create a national biomass map for China (Chi et al. 2015). However, Milenković et al. (2017) cautioned about the effects that changing footprint size would have on forest biomass estimates, specifically when using spaceborne waveform LiDAR. Anderson et al. (2016) highlighted the benefits of waveform LiDAR for investigating vegetation dynamics across a range of ecological systems.

Waveform LiDAR for phenological investigations
However, the full potential of waveform LiDAR as a tool for phenological studies has yet to be realized and is currently unpopular when based on the number of literature published thus far. Articles that are primarily focused on LiDAR for forest phenology application are rare. Table 3 shows the keywords and the number of results using Google Scholar. While there are limitations of the journals that are included and indexed by Google Scholarfor example, Elsevier journals were not included before the middle of 2007the list of scholarly literature across an array of publishing formats and disciplines is sufficient enough to provide a general overview of the popular and unpopular fields of study.
Considering the first 50 results, only four papers are somehow associated with LiDAR applied to phenological studies. It is even much worse when "waveform LiDAR" and "phenology" are searched together. One related paper looked at differentiating coniferous and deciduous tree species using two different seasonal leaf-off and leaf-on LiDAR datasets (Kim et al. 2009) from the temperate forest zone in southern US. The authors also investigated the distinction between evergreen coniferous and broad-leaved deciduous trees, analyzed LiDAR intensity values from trees with different foliage characteristics, such as the presence or absence of foliage. Collin, Long, and Archambault (2010) highlighted the strong relationship (r 2 = 0.87) between a new index called the Normalized Difference LiDAR Vegetation Index (NDLVI) and the salt-marsh vegetation measured in-situ. NDLVI was derived from a dual-wavelength LiDAR using inverted traditional NDVI equation: where R max is the intensity value of the red or "Raman peak" (Raman 1928) at 645 nm red LiDAR channel. IR max is the intensity value of the NIR peak at 1064 nm LiDAR channel. The underlying principle of the NDLVI was to enhance the greenness of plant communities without predominance of chlorophyll b (the red LiDAR channel corresponds to the second absorption peak of chlorophyll b). Another paper published by Simonson, Allen, and Coomes (2018), studied four LiDAR variables (maximum, mean, standard deviation and skewness of vegetation heights) and how they were influenced by seasonal canopy changes. The specific contribution of the study, however, was in the investigation of subtler effects of tree leafing phenology rather than the leafon/leaf-off dichotomy. The work by Calders et al. (2015) provided insight in the stability of terrestrial LiDAR measurements when monitoring changes in spring phenology. Calders et al. (2015) determined the exact SOS dates for two different species composition of forests. To improve species classification, Hovi et al. (2016) acquired waveform LiDAR data during early summer and showed how intra-species variations in crown and branch morphology and vigor were logically linked with waveform features. Finally, a related study that used simulated waveforms also showed the full-waveform canopy LiDAR capturing the seasonal and vertical variations of NDVI (Morsdorf et al. 2009). The inadequacy in terms of the number of scientific studies conducted using LiDAR regarding monitoring or even modeling of forest phenology comes with a wellfounded rationale. Based on the preceding section of this paper on applications of waveform LiDAR, the system has been efficient at case studies that focus on a static condition (e.g., vegetation structure, height) at a single point in time. Also, there are compelling issues in avoiding LiDAR for phenological monitoring and modeling and just going for passive remote sensing otherwise. First, the waveform LiDAR data files are very large. Since phenological studies require repetitive measurements, LiDAR data could quickly fill up computer servers and drives. Second, the need for effective access to and storage of scan data, coupled with the absence of a universal format standard, has caused developers of LiDAR software to implement their own storage format, which, pay little attention to enabling import and export options. The file format of the LiDAR depends on the form of representation. In its raw form, LiDAR is a series of points stored as x, y, z where x and y could be longitude and latitude; z is the intensity of return. A simple ASCII file where each line has a coordinate (x, y, z) separated by commas (tabs, spaces, etc.) could be used to represent the data and easily accessible with a text editor. Only recently a binary file format (LAS) endorsed by the American Society for Photogrammetry and Remote Sensing (ASPRS) for the interchange of LiDAR data has been gaining popularity and support. The.LAS file format addresses issues on data sharing and flow.
The second issue involves the frequent return time or the acquisition timing that should be anchored in the seasonal progression and phenological stage of vegetation. Take for instance the coarse-resolution passive sensors that allow for near-daily monitoring of phenology; an important consideration for landcover classification and detection of future climate change (Morton et al. 2005;Goetz et al. 2005). The frequent return time of MODIS (Table 2) permits phenology-based mapping of tropical deforestation with 89% accuracy. Repeated passes that capture seasonal leaf dynamics could improve the detection of seasonal tropical forests (Morton et al. 2005). Although flights for airborne laser systems could be planned for daily LiDAR data gathering, space-borne waveform LiDAR systems could be limited to 14 days of orbit repeat period (Dubayah et al. 1997). In the case of NASA's proposed GEDI, it is dependent on the orbit cycle of the International Space Station (ISS).
The ability to fly during all hours of the day (i.e., during the night, as it is independent of solar illumination) makes LiDAR better than optical remote sensing in this regard. However, rainy days should not be considered because infrared light does not penetrate water vapor. LiDAR data acquisitions during leaf-on and leaf-off conditions may be shelved when the acquisition window occurs at a period of optically thick cloud or continuous precipitation (Drake et al. 2003). Thus, as in the case for airborne LiDAR, the best way to optimize the timing for data acquisition to monitor vegetation conditions (tree leaf-on, leafout, senescence, etc.) is to plan ahead of the anticipated flight days and schedule longer time windows.
Amidst the issues of LiDAR application to phenology, more scientists are beginning to utilize LiDAR datasets (Calders et al. 2015). As technological capabilities accelerate, the future of geospatial research is on the integration of different systems and sensors, from the ground to space-borne (Dold and Groopman 2017;Parshin et al. 2018). That being so, the next bulk of research may focus more on linking both the structural (waveform LiDAR) and spectral signatures (optical data) for forest canopy characterization (Fang et al. 2018). Waveform LiDAR in itself may offer novel ways that could help expound vegetation phenological processes, as sensor application and machine computing converge to accurately capture real-world dynamics. The following sections would examine a potential methodology for the phenological study of vegetation using waveform LiDAR.

Potential vegetation phenology analysis using waveform LiDAR
Based on waveform simulations, Koetz et al. (2006) used the updated 3D waveform model to calculate the Canopy Fractional Cover (CFC) or the proportion of the landscape occupied by green or non-green vegetation. The model constructed a 3D-representation of the observed forest stand taking into account the number and position of trees, tree height, crown geometry and shape as well as the exposition of the underlying topography. The crown itself was described as a turbid scattering medium parameterized by its foliage area volume density, the Ross-Nilson G-factor (Nilson 1971) and the foliage reflectance. The updated model also allowed for the input of LAI instead of the foliage area volume density. Results of the study showed that inversion of the LiDAR waveform radiative transfer model provided reliable estimates of model parameters describing the vertical as well as the horizontal canopy structure when linked to a synthetic dataset derived from a forest growth model ZELIG (Urban 1990). The parameters were all retrieved with an overall high correlation coefficient and low RMSE (r 2 = 0.92 and RMSE = 0.08, respectively).
LiDAR models of canopy fractional covers have been tested over several LiDAR datasets collected across a variety of forest ecozones and canopy structural classes (Hopkinson and Chasmer 2009). The models use sensor returns or the amount of backscatter elements that are strong enough to register a distinct amplitude of reflected energy. There are types of returns in terms of peak property (Miura and Jones 2010): type 1 (singular returns), type 2 (first of many returns), type 3 (intermediate returns), and type 4 (last of many returns). Single returns are for those for which there is only one dominant backscattering surface encountered. Multiple returns could give a picture of the dominant canopy and ground elements along the pulse path. LiDAR pulses are intercepted by the top canopy, lower canopy vegetation, low-lying understory and/or the ground surface. In the case of the ground surface, it is inevitable that the returns from the ground have passed through canopy gaps (Hopkinson and Chasmer 2007). Many gaps within the canopy will result into a fractional cover of close to zero, whereas few gaps will give CFC a value close to unity. These assumptions are based from Beer-Lambert Law (equation 4) with gap fraction (P) equivalent to the transmittance (T): where I 0 is open skylight intensity above canopy, I l is the light intensity after traveling a path length through the canopy and k is the extinction coefficient, which can be approximated to a value of 0.5 in a canopy of spherical leaf distribution (Martens, Ustin, and Rousseau 1993) but generally varies between about 0.25 and 0.75 for the natural needle-and broad-leaf canopies (Jarvis and Leverenz 1983). Some authors (e.g. Barilotti 2006) used the total reflected energy from the canopy to ground profile as being some proportion of the total available laser pulse intensity, and the reflected energy from ground level as a similar proportion of the transmitted pulse energy. Adapting equation 4, the return power relationship is given by equation 5: where ΣI b is below canopy power (the sum of all ground return intensity) and ΣI t is the total power (sum of all intensity) for the entire canopy to ground profile. Another model was created to highlight the incident pulse intensity as it enters the canopy. This is one way to estimate canopy FC at near-nadir angles. Adapting again equation 2, CFC can be estimated from the ratio of the sums of canopy to total return intensities (equation 6).
where ΣI canopy is the sum of all canopy level return intensity and ΣI total is the sum of all return intensity for the entire canopy to ground profile. To take into account the two-way energy transmission loss (both into and out of the canopy), the CFC can be computed as the square root of the proportion of the ground level energy reflectance to the sum of all return intensity for the entire canopy to ground profile. The Beer's law modification of equation 6 will become equation 7: where ΣI ground is the sum of all below canopy return intensity. Among the models tested by Hopkinson and Chasmer (2009), the canopy-to-total ratio was the one that demonstrated a strong correlation to the canopy fractional cover from field data. The authors also suggested that a combination of equations 6 and 7 could provide a more robust estimation of the CFC across a range of canopy classes. While it looks like that these equations only highlight the utility and applicability of multiple discrete return intensities as defined by Miura and Jones (2010), they also function for waveform LiDAR return intensities.
Take, for example, Figure 2 shows a single waveform. Instead of only taking specific canopy returns at x 1 , x 2 , x 3 , and x 4 (Figure 2(a)), ΣI canopy could become the sum of all canopy level return intensity at a specified threshold, going from left to right (Figure 2(b)). This method will not only limit the calculations in terms of the dominant backscattering, rather, it will define other elements along the pulse path by taking into account other significant waveform pulses. ΣI total could become the sum of all return intensity for the entire waveform extent. Modifying equation 6 will result to: This proposed optimized technique could extract inherent information of the waveforms by the integration of single reflections. In one paper, waveform decomposition was introduced by Gaulton et al. (2010) to find almost all reflections of the wave, enabling detection and resolution of small dominated and young generation trees. The model decomposes the LiDAR waveform in a region of interest by fitting a series of Gaussian pulses to the waveform. It utilizes coordinate positions of the reflections and optimally the pulse width.
This concept of using pulse width (Gaulton et al. 2010) could be integrated with equation 8 to account for the likelihood of waveform return intensity similarities.
but I canopy 1 = product of dominant intensity and its pulse width or I*w and I canopy 2 is the intensity of all other non-dominant pulses within the specified threshold.
Allowing the inclusion of the pulse width, equation 9 becomes: The derived fractional cover from equation 10 is still in a static condition or this P represents a single LiDAR FC in a specific time. To employ phenology means to look at the seasonal variation of the LiDAR-based FC, which then could be illustrated in Figure 3. Fractional vegetation cover is strongly correlated with NDVI (Carlson, Perry, and Schmugge 1990), even deriving the fraction of green vegetation from NDVI (Gutman and Ignatov 1998).
As mentioned in the previous section of this paper, what may be seen as a deterrent in the use of LiDAR for phenology is the frequent return time or the acquisition timing of the waveform LiDAR sensor. Temporal resolution of the planned space-borne LiDAR system is about 2 weeks. This timing, however, is still effective in estimating the start of season and End Of Season (EOS), which are considered as ecosystem-level, not species-level, phenological metrics (White et al. 2002).
There are several phenological leaf onset models that have been published for tree species (e.g. Cannell and Smith 1983). The existence of a number of models shows that the nature of the different phenological phases and their physiological and morphological status is not yet fully understood. Each of the numerous quantitative methodologies is only suited to its specific research question. While some methods are complex, this paper will attempt a simple approach to model phenology using waveform LiDAR product. With the aid of Figure 3, an algorithm could be formulated to estimate the SOS and the EOS from a given LiDAR CFC time series.
For the approximation of a given LiDAR CFC, the piecewise defined continuous spline function S (t) could be applied (equations 11 to 13). An almost similar approach, a triangular function, was used by Häninnen (1994) and Schaber and Badeck (2003).
S t ð Þ ¼ N 1 otherwise (dormancy phase, DP)(13) The three splines represent the three different phenological phases (Figure 4) using an increasing parabolic function for the leaf-on and a decreasing parabolic function for the senescence phase. The function is dependent on five parameters: t 1 (time 1), t 2 (time 2), t 3 (time 3), N 1 (LiDAR FC 1), and N 2 (LiDAR FC 2). The locations of these parameters and the quantitative amount of LiDAR CFC could be resolved from Figure 4 since the plot of LiDAR CFC vs. DOY is not smooth, and could sharply detects time transitions. This model approach is not complex and also less sensitive against circumstances where there is considerable noise in the LiDAR FC curve, in contrast with methods based on a modified second derivative of the NDVI.
As a warning, this method has not been tested yet and needs to be calibrated with other existing phenological models. Statistical analysis, such as minimizing the absolute error, may be applied to measure the quality of fit.
The same process above could be applied when trying to assess the phenology of the canopy understory or the discrimination between younger and older stands as well as the discrimination between dense and thinned stands. Koetz et al. (2007) simulated LiDAR waveforms as a function of the forest stand structure plus the sensor specifications, and put them against the results of a forest growth model and the reflectance spectra from the imaging spectroscopy. Results showed that different forest growths could lead to differences in the shapes of the LiDAR waveforms. This was not the case, however, for the reflectance where different combinations of forest heights/ structure give rise to similar spectra (Buddenbaum and Seeling 2008).
This characteristic of LiDAR to define the structure could be exploited by considering equation 10 and replacing the canopy intensity parameter by, say, ground/soil parameter or understory parameter. A similar temporal CFC may be plotted and monitored for SOS and EOS phases. This is one way to measure and monitor processes such as forest growth.

Satellite sensors to link LiDAR and phenology
Few LiDAR systems have ever flown in space, owing to limitations involving high power, high cost (Farid, Goodrich, and Sorooshian 2006), and the availability of robust laser sources. There used to be only one satellite-mounted LiDAR sensor, ICESat GLAS (the Geoscience Laser Altimeter System, carried on the Ice, Cloud and Land Elevation Satellite), and although its large footprint (~65 m) was not explicitly designed for vegetation applications, the mission is being used to estimate forest biomass and improve DEMs (Harding and Carabajal 2005; NASA (National Aeronautics and Space Administration) 2009; Nelson et al. 2009). Unfortunately, the ICESat mission is no longer collecting data as the final GLAS laser ceased firing in 11 October 2009. Currently, there are no space-borne LiDAR systems designed specifically to measure vegetation canopies. There had been an attempt like NASA's Vegetation Canopy LiDAR (VCL) mission in the 1990 s. VCL, however, never made it to space. The VCL would have been the first multi-beam waveform-recording LiDAR to fly in spaceholding five lasers, each VCL orbit could sample an area five miles across. Due to serious mission viability concerns like costly trials, and the development dropping way behind schedule, the mission was canceled. Then, there is a NASA-led, proposed Global Ecosystem Dynamics Investigation (GEDI) -3D LiDAR that is scheduled for late 2018 launch. The motivation of the GEDI mission includes the responses of terrestrial biomass, which stores a large pool of carbon, to changing climate and land management. Although a 25-m footprint size is not considered a high resolution compared to a small footprint LiDAR, the scientific objectives of the mission include modeling finer-scale structure (Stysley et al. 2015). Moreover, GEDI is designed to collect LiDAR waveform observations for the world's forests during cloud-free periods (Stysley et al. 2015;Qi and Dubayah 2016) over 14 parallel tracks. However, these observations would only be within the orbit of the ISS (i.e. 50°N to 50°S), thereby excluding much of the northern boreal forest. When launched, GEDI would be the first satellite mission to provide global, highresolution observations of forest vertical structure. Other structure information that could be extracted from the GEDI waveform includes surface topography, canopy height metrics, and canopy cover metrics.
Another proposed vegetation LiDAR sensor is the Japanese International Space Station (ISS)/Japanese Experimental Module (JEM)-borne named MOLI (Multi-footprint Observation LIDAR and Imager) (Murooka et al. 2013). Similar to the GEDI, MOLI has a full-waveform pulse type that could distinguish canopy heights, crown depth, and other forest features. Apart from the LiDAR capability, MOLI has a high-resolution imager that consists of RGB bands, a swath width of 500 m, in a ground resolution of 5 m. The imager is adopted to provide information on tree crown size, height and field data for the conversion to AGB (Murooka et al. 2013). MOLI is planned for 2021 launch just after NASA/GEDI finished a twoyear experiment (Asai et al. 2018).
The U.S. National Research Council (NRC), in a document known as the Decadal Survey (National Research Council 2007), has identified another LiDAR space mission of paramount importance to the U.S. scientific community for monitoring the status and function of the biosphere. Apart from the GEDI, the ICESat-2 was launched in September 2018. Like its predecessor ICESat-1, the main science goals of ICESat-2 are specifically targeted at measuring ice sheet changes and sea ice thickness. It also supports multidisciplinary applications and measures vegetation biomass. ICESat-2 operates with three pairs of beams, each pair separated by about 3 km cross-track, spaced at 90 m (Markus et al. 2017). Each of the beams has a footprint of 17 m in diameter with an along-track sampling interval of 0.7 m. Instead of using an analog full-waveform type, ICESat-2 employs the photon counting Advanced Topographic Laser Altimeter System (ATLAS). In this LiDAR system, the arrival time of each photon is recorded within the vertical distribution of the reflected signal, resembling a full waveform. ICESat-2 data would be time-tagged photon elevations as a function of distance along-track .
ICESat-2 would be capable of producing a global vegetation height surface with 3 m accuracy at 1 km spatial resolution. With measurements that span a minimum of 3 years, ICESat-2 would contribute to large-scale biomass assessment, would help understand the global distribution of carbon on land, and would complement those of the GEDI and MOLI space missions, whose orbit and LiDAR component is more specifically targeted at vegetation and ecosystem science objectives. The next section of the paper would discuss system specifications of the waveform LiDAR products that may have an important bearing on phenological monitoring and modeling.

Laser footprints
The proposed GEDI and MOLI products could provide transects of vegetation vertical canopy profiles overall biomes at 25 m spatial resolution (Table 4). A laser footprint of 25 m has been designed to adequately resolve vegetation height and structure. This is seen as one of the advantages of the large-footprint LiDAR systems over small-footprint LiDAR systems. While small-footprint LiDAR may under-represent the canopy structure, especially the true tree height due to the sensing of portions of forest elements, GEDI and MOLI guarantee a reflection from the top of each canopy within the sampled area as well as intra-tree and inter-tree gaps. For high-density tropical forests where canopy diameters could reach to 10 m to 25 m, footprint size of GEDI and MOLI is suited the best.

Orbit repeat period
The revisit period for space-borne waveform LiDAR helps in understanding better temporal canopy variability. Higher frequency measurements provide time-varying information of leaf-off and leaf-on. When the waveform LiDAR datasets are acquired in leaf-off conditions, the LiDAR intensity information could be used to differentiate between deciduous and conifer trees (Reutebuch, Andersen, and McGaughey 2005). The near-infrared reflectance from branches and stems is much lower than that from live foliage; therefore, the intensity of the LiDAR reflections from leaf-off deciduous tree crowns is significantly lower than that from leaf-on conifer crowns (Andersen 2009 50°north and south latitude, and is in an approximate 4-day revisit cycle (Qi and Dubayah 2016). The limitation of GEDI is that it only gets 2 years on the ISS, after which MOLI continues the observations for another year or two.

Swath width
Different systems utilize different scan angles that could allow the collection of data in a much larger swath. Swath width of GEDI and MOLI could provide unequaled decimeter level vertical absolute accuracies of "true ground" even in highly vegetated regions. Space-borne wide-swath LiDAR with fullwaveform collection capability could provide the following (Blair, Hofton, and Luthcke 2001): landscape scale (10 km swath) coverage, full earth phenological coverage at <10 m pixels within 1 year, near-100% coverage/illumination, change detection measurements at sub-centimeter relative vertical accuracy, subtle topographic change beneath vegetation, and vegetation and land cover changes.

Beyond space-borne LiDAR
Associated with the failed VCL mission is the operational LiDAR Vegetation Imaging System (LVIS) (Blair, Rabine, and Hofton 1999). It is an airborne, wide-swath, large footprint simulator developed by NASA's Goddard Space Flight Center. LVIS collects waveforms using a 25-m footprint over a 2-km swath width. Figure 5 shows examples of the waveforms derived from the LVIS instrument over the La Selva Biological Station in Costa Rica in two different years, 1998 and 2005. The scientific community could look into the location of maximum canopy and soil returns (in terms of time) and their temporal shifting. The weaker ground returns in Figure 5b are caused by high canopy closures and low-lying vegetation. It is worthy to mention that this is one of the advantages of waveform LiDAR over 2D images: the vertical component of 3D LiDAR data allows the researcher to separate ground vs. vegetation information, which is a prerequisite to most LiDAR applications, many of which focus exclusively on one or the other component. Also, in the absence of space-borne waveform LiDAR, data from the LVIS can be exploited to test the phenological concept presented above.

Conclusion
The body of knowledge and literature on useful LiDAR applications has rapidly expanded in recent years. This means that the remote sensing community is gaining interest in the ability of LiDAR to measure topography, capture vegetation height and cover, and record complex profile of canopy height. LiDAR has limitations in the arena of phenology monitoring and modeling. In the case of a waveform system, one limitation is the spatially sparse distribution of the target on the ground that could give an inaccurate reconstruction of the model. This limitation could likely to be resolved by the combination of LiDAR with additional data from spatially extensive sensors such as high-resolution optical images or radar data. Fusion of data from different sources requires the understanding of the different data platforms, dates of acquisition and resolutions. Apart from fusing multiple sensors, the fusion of LiDAR with different modeling algorithms such as the Support Vector Machine (SVM), Random Forest (RF) and Extreme Learning Machine (ELM) could also build more robust empirical relationships that could be extrapolated to areas where datasets are sparse. Another limitation is the lack of precise algorithms and approaches for operational use of the data, especially for phenology. This paper has attempted to illustrate the prospect of using waveform LiDAR for estimating the onset of leaf, including the few LiDAR issues that may deter its full applications. It has been presented in here through a vast list of literature that the waveform LiDAR product may not only show potentials in estimating and mapping a wide range of information such as canopy height and vertical distribution profiles, but may also be exploited using its distinct characteristics to monitor and model leaf phenology. The new model uses a spline function; however, a combination of linear or even triangular functions may result in more defined estimates.
Although the model utilizes CFC in its estimation procedure, there are also options for translating it to LAI or other vegetation indices when necessary in the analysis. The LiDAR research approach could put importance on phenological assessments to measure and monitor processes such as forest growth and eventually carbon sequestration. One of its strengths is the ability to look at leaf onset from the canopy and the understory, even to the non-green elements. There is also a possibility of looking at the spatial variation in the timing of phenology, which occurs within the woodland. Further, the algorithm uses return intensity, an attribute that describes the strength of the beam backscattering. It is dependent on the reflectance characteristics of the target; therefore, it could potentially be used in target discrimination. However, it is unclear how much the variation in accurately predicting the CFC is due to differences, for instance, in LiDAR sensor types or variability among biomes. This approach has to be thoroughly tested for all types of waveform shapesvery complex, slightly distorted, asymmetric, flattened and peaked.
At present, there is no existing space-borne LiDAR system designed specifically for vegetation canopies. The hope is placed on the proposed GEDI, MOLI, and ICESat-2 missions. The greatest limitation of these instruments is the acquisition timing that should be based on the seasonal progression and phenological stage of vegetation. However, at the ecosystem-level, the revisit periods of GEDI and MOLI are still effective in estimating SOS and EOS.
With the wider utilization of LiDAR data usually hampered by lack of familiarity with the waveform measurement approach, more models based on the waveform could be proposed and developed to exploit the rich information content of such a dataset, especially in the connections between the LiDAR signal and foliage parameters. Furthermore, new predictive models for canopy fuel parameters estimation should also be looked into and how the new modeling concept presented in this paper can be integrated.

Notes on contributor
Eric Ariel L Salas received his BS in Civil Engineering from University of San Carlos in Cebu, Philippines, and his MS in Geo-Information from Wageningen University and Research in Wageningen, The Netherlands. He earned his PhD in Geospatial Science and Engineering with specialization in Remote Sensing Geography from South Dakota State University in Brookings, USA.