• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Using GEDI lidar data and airborne laser scanning to assess height growth dynamics in fast-growing species:a showcase in Spain

    2021-04-14 06:59:12JuanGuerraHernndezandAdriPascual
    Forest Ecosystems 2021年1期

    Juan Guerra-Hernández and Adrián Pascual

    Abstract

    Keywords: Forest inventory,Productivity, Multi-temporal LiDAR, Full-waveform (FW), Satellite mapping

    Introduction

    Forest inventories employing airborne laser scanning(ALS) data have become common in many countries(Wehr and Lohr 1999; Lefsky et al. 2002; N?sset et al.2004;White et al.2016).The ALS-based forest inventory methods are typically categorized into two groups: the area-based approach (ABA) (N?sset 2002; Guerra-Hernández et al. 2016) and individual tree detection(ITD) (Hyypp? et al. 2001). So far, most operational forest inventories framed under ALS data utilization have been designed under the ABA method (Maltamo et al.2014). Under the ABA method, stand attribute models are fitted with sample plots using metrics calculated from ALS data and then these models are used to predict stand attributes wall-to-wall for whole inventory area, typically using square grid cells whose area exactly matches the size of the sampling plots as inventory units(N?sset 2002). As a result, the elevation of tree vegetation or the predicted forest inventory variables can be expressed as a spatially continuous raster map. The science and the implementation of ALS-based forest inventory is solid, evolving fast and mainstream worldwide(Maltamo et al. 2014).

    During the last decades, nationwide laser scanning has helped forest managers to integrate LiDAR in forest inventory (Nilsson et al. 2017). The impact of ALS data is strong on forest decision-making when data from National Forest Inventories is estimated with remote sensing sources under good co-registration goodness(Pascual et al. 2019). The number of countries that have been scanned at least once relying on ALS technology has increased and in many areas forest dynamics and natural disturbances can be estimated using multitemporal LiDAR (Hopkinson et al. 2008). The time-gap between ALS surveys, the growth rates of the species and forest management actions challenge the management of fast-growing plantations (McRoberts et al.2018). As a showcase, the rotation age for Eucalyptus stands in the Northwest of Spain ranges from 13 to 20 years while the time gap between the last two nationwide ALS surveys for the area was half the rotation (7 years). Therefore, capturing the fast growth rates under commercial silviculture can benefit from independent and flexible sources of 3D data operated on a shorter interval basis compared to ALS nationwide surveys(Silva et al.2017).The launch of the NASA’s Global Ecosystem Dynamics Investigation (GEDI) satellite mission in 2018 has paved the way to the integration of global, multitemporal and full-waveform laser data into operational forestry (Duncanson et al. 2020). The estimation of forest structure using GEDI laser statistics expands the frontiers of forest monitoring and support a change of paradigm on the spatio-temporal evaluation of forest productivity (Bontemps and Bouriaud 2013; Coops et al.2018).

    The GEDI missions is continuously scanning over the earth’s land surfaces since 2019 using optical Light Detection and Ranging (LiDAR) to monitor tropical and temperate forests ecosystems (Dubayah et al. 2020a).The billions of GEDI laser shots are expanding the frontiers of forest monitoring worldwide especially in remote forest areas (Silva et al. 2017; Duncanson et al. 2020).The GEDI laser instrument generates total of 8 beam ground tracks that are spaced approximately 600 m apart on cross-track direction. Each beam track consists of ~25 m footprint samples spaced every 60 m along the ground track. The GEDI data science products include arrays of laser waveform metrics for each footprint sample, including canopy height and profile metrics on above-ground height or ground elevation (Hancock et al.2019;Dubayah et al.2020a).The understating and development of solutions based on GEDI data is a promising research topic (Duncanson et al. 2020).One of the applications of the GEDI technology is the enhancement of large-scale forests maps generated with ground-based National Forest Inventories (NFI).

    The aim of the investigation was to test the capability of GEDI to support ALS scanning when it comes to describe dynamics in fast-growing forest ecosystems. The 2.5-year time gap between ALS surveys and the first GEDI orbits in the Northwest of Spain reflects the scenario in which GEDI laser shots can be used to update forest structure and detect spatio-temporal changes.This investigation is one of the first attempts to validate the capability of GEDI to detect spatial and temporal changes using nationwide ALS surveys as baseline. We tested the capability of GEDI data science products under challenging fast-growing commercial forest ecosystems and the small-scale forestry carried out in Galicia, a mosaic of private and public forest patches (2 ha on average).

    The objectives of this study were:i)to use GEDI to detect disturbances at regional level, ii) to investigate,evaluate and explain the source of uncertainty when interpreting GEDI footprint ground shots and iii) to validate the height growth dynamics detected from ALS to GEDI compared to ground NFI data. The paper emphasizes on the filtering mechanisms and geospatial operations required to combined ALS data and GEDI data.The GEDI laser data will become a reference to monitor forest dynamics worldwide once researchers and practitioners fully understand the technology and its integration under operational forest monitoring.

    Material and methods

    Training area

    Fig.1 Location of the training area in Spain in the region of Galicia.The extent of the three forest ecosystems assessed in the province of Lugo is displayed using the Forest Map of Spain(FMS)

    The training area is the province of Lugo, located within the region of Galicia in the Northwest of Spain(Fig. 1). The landscapes are characterized by rugged orography and mild temperatures with low thermal oscillation between winter and summer and frequent rainfalls. The area is a rural region and forestry is one of the principal economic activities. The conditions of the region boost the high productivity of the main tree commercial species in Galicia and even Spain: Eucalyptus globulus Labill., Pinus pinaster Ait.and Pinus radiata D. Don. Galicia is one of the most important regions in Spain in terms of forestry production. To overview the economics, around 6.5 million m3are harvested every year and the turnover is steadily growing due to the intensive use of forest biomass as supply of the pulp industry, a major driver of the economy in Galicia. Around 30% of the 1.4 million hectares of forest land comprises pure and even-aged pine stands, mainly maritime pine (P. pinaster) (271,281 ha) and radiata pine (P. radiata) (96,(916,71777 ha) (MAGRAMA 2018). The high productivity of P. pinaster and Pinus radiata forests, for instance,allows profitable timber production with short rotation lengths. The integration of laser scanning technology in the region is becoming mainstream in operational forest inventory (Gon?alves-Seco et al.2011; González-Ferreiro et al. 2012). However, the short rotation and the fast dynamics under commercial forestry have constrained the impact of nationwide ALS surveys on the decision-making arena. The use of 3D data collected on shorter schedules can support the industry while giving forest owners of small-scale forest properties the chance to monitor remotely located assets. The fragmentation of the landscape, the fast growth of forest species and the silviculture turn Galicia into a laboratory to test the measuring capability of GEDI laser technology.

    The Forest map of Spain

    The last version of the Forest Map of Spain (FMS) at 25-m resolution was used to select the forest areas corresponding to the most important commercial species in Galicia (MAPA 2018). The FMS contains wall-to-wall information at polygon level extracted from the National Forest Inventory (NFI) including species occupation, forest stage of development or forest cover, which was assessed by photointerpretation of aerial RGB imagery from the PNOA project (Plan Nacional de Ortofotografía Aérea) and ground NFI data. The integration of data science products from GEDI and ALS surveys can improve the description of FMS units (polygons). The scope of the study was narrowed to a bounding box of 73.2 km width and 140.6 km long (Easting 588,902.2-661,214.47, Northing 4,695,137.4-4,835,977.0 UTM)covering around a million hectares, which is 2% of the country and almost the full extent of the Lugo Province.The research evaluated the three main fast-growing species in Spain: P. pinaster, P. radiata and Eucalyptus spp.The final scale of the assessment after filtering data reached a total of 211,346 ha in three groups: 47,820 ha(P. pinaster), 72,620 ha (P. radiata) and 90,906 ha (Eucalyptus spp.).

    Airborne laser scanning data

    The ALS surveys are part of the PNOA project (MATA 2019), which allows the use of nationwide laser data for free. The nominal pulse density of the 2015-2017 ALS surveys, was 1.0 pts.·m?2on average with a measuring accuracy of 20 cm in the vertical profile. The scheduled calendar for the flightpaths was available. The scale of the LiDAR processing was regional (Galicia) and the FMS polygons were used to account for the spatial distribution of the three species across the region. The raw LiDAR data files retrieved from the PNOA server comprised 282 GB of data. The licensed software Lastools(Isenburg 2020) was used to process the point clouds and compute 3D laser statistics using Intel?Core i9-7900X workstation with 64 GB of RAM and 10 cores. As a guide for the reader, we recommend the papers from Wehr and Lohr (1999) and Lefsky et al. (2002) to understand the principles of ALS and our previous study on ALS processing steps using large-scale 3D point clouds(Pascual et al. 2020). Briefly, lasheight were used to normalize the classified point cloud of ALS echoes. The lascanopy command was used to extract the 99th height percentile and CC (canopy cover) metrics from the normalized point cloud using a circle of 25 m of diameter for the center of each GEDI footprint ground shot, while lasclip was used to extract the normalized point cloud for the extent of FMS polygons. Then, lascanopy was also used to generate a 25-m raster grid of the 99th height percentile for the training area. The height_cutoff and cover_cutoff parameters were set to disregard ALS echoes below 2 m when computing ALS statistics. The 25-m grid raster for the 99th height percentile (ALSH99)was generated for the three species within the extent of the bounding box created to retrieve GEDI data. The grid scale of 25 m matched the ground footprint of GEDI laser shots (Dubayah et al. 2020b) to avoid the bias resolution when comparing laser statistics.

    Filtering operations were very important along the research project as they were applied to process the ALS surveys, to process the GEDI laser shots and to combine both while acknowledging the edge effect of the irregular boundaries of FMS polygons and the resolution dependence of the computed statistics. The notation for filters along the manuscript follows their chronological implementation in the research workflow.Shots GEDI not completely contained within FMS polygons in the training area were removed to avoid edge effects(Filter#1).Statistics from the 3D point clouds were computed for each grid cell. Growth dynamics under multi-temporal laser scanning are better detected when using the highest height percentiles on ALS echoes distribution(Hopkinson et al. 2008; Socha et al. 2020). The proper benchmarking between GEDI and ALS statistics should be carried out using ALS point cloud versus GEDI FW to avoid resolution dependance and the averaging effect on ALS statistics when using raster pixels (Packalen et al. 2019). From the center points of the valid GEDI laser shots considering the 25-m footprint, we clip the 3D point cloud using only the ALS echoes within the ground GEDI footprints.Hereby,the 99th height percentile(ALSH99)from the ALS data was computed for each GEDI shot.

    The GEDI satellite LiDAR

    The bounding box was used to select the GEDI beams within the 1-million-ha training area. The rGEDI package(Silva et al. 2020) was used to retrieve and process the GEDI data under the 3.6 version of the R statistical software (R Core Team 2020). There were 31 GEDI orbit tracks in h5 format available to date March 15th, 2020(Fig. 1). The GEDI laser shots within the bounding box considered in the research were collected in 2019 between April 21st and August 3rd. The GEDI laser equipment uses 3 lasers to produce 8 laser beams creating laser shots every 60 m along the longitudinal track in the ground.The distance in the cross-track direction for two consecutive longitudinal tracks is 600 m for a given orbit. When combining GEDI shots from different acquisition dates and orbits, the cross-track separation between tracks decreased. For instance, the study area contained distances in the cross-track direction down to 60 m when we merged the multi-temporal GEDI shots collected in ten different days along the specified time interval. To date,the stage of development of the GEDI mission restricts the possibility of wall-to-wall mapping across our training area (Hofton and Blair 2020) for which the level of landscape fragmentation is high(Fig.2).

    The use of GEDI Level 3 to produce wall-to-wall estimates should be carefully assessed when evaluating small-scale properties as occurs in Galicia. Hereby, the authors framed the research using the elevation statistics in Level 2A data science product (Dubayah et al. 2020b).The full-waveform (FW) for each GEDI laser shot was processed to compute the following statistics for each of the 25-m footprint ground shot. The focus of using Level 2A was to compute the same height percentiles as with the ALS data. The relative height was computed along 100 unitary intervals from the minimum to the maximum elevation (RH100), and RH99was the selected statistic on which to benchmark ALS data.

    Geospatial operations and statistics with GEDI, ALS and the FMS

    Fig.2 Partial overview of the P. pinaster forests in the training area showing the aerial satellite image and canopy height model derived with airborne laser data in the background.The ground beams of the GEDI shots(GEDI_level2A)intersecting forest areas are presented

    Fig.3 Spatial layout of Filters#2-8 use to disregard GEDI shots displayed sequentially. a Filter #2,b Filter#3, c Filter#4-6 and d Filter #7-8.The maps correspond to forest areas classified as P.pinaster. The raster map for ALSH99 (99th height percentile,LP99 in the legends)derived from ALS surveying is shown over the 2017 orthophoto

    The benchmark between ALS and GEDI was carried out within the boundaries of FMS polygons intersecting GEDI ground tracks (Fig. 3a). For a given GEDI laser shot, we derived the RH99from the 25-m footprint ground shots (GEDI99), the 99th percentile of the ALS echoes height distribution within the extent of the 25-m diameter samples (ALS99). From the initial set of 288,529 laser shots within the bounding, we applied the following filters:

    1) Selection of GEDI shots completely contained within the species-specific FMS polygons(Filter #2).

    2) Selection of GEDI shots located further than 30 from the edge of the ALS-based raster map boundaries (Filter #3).

    3) Selection of GEDI shots whose GEDI99was above 70 m (Filter #4).

    4) The column quality flag was used to disregard all GEDI shots classified as 0 meaning (Filter#5) the technical and quality attributes for the given shot number were not according to standards. See detail of interpretation of L2A Quality Flag in Dubayah et al. (2020b). A quality flag value of 1 indicates the laser shot meets criteria based on energy,sensitivity,amplitude, and real-time surface tracking quality.

    5) All pairwise values of GEDI99and ALSH99ranged between 2 and 70 m to remove outliers (Filter #6).The upper bound was set based on the ground NFI tree height measurements across Galicia for the three species.

    6) We disregard GEDI shots for which the canopy cover (CC)associated ALSH99was below 70% to reduce the impact multi-layered and sparse forest conditions on the GEDI laser reflectance homogeneity (Filter#7).

    7) The average height increment (ALSH99-GEDI99)and its standard deviation was computed for FMS polygons comprising 3 or more GEDI shots within their boundaries. Preliminary tests showed the importance of this filter to remove outliers for very large positive increments (Filter #8).

    The final selection of GEDI shots were used to compute height differences from GEDI to ALS statistics corresponding to the ground extent of GEDI shots. As a showcase, we display the geospatial operations from Filter #2 to Filter #6 for a sample of P. pinaster stands in the training area (Fig. 3).

    The estimation of height growth dynamics required from the exact date for which ALS data was collected for all GEDI points. For the set of valid GEDI shots, we computed the difference in days and in years between ALS survey acquisition retrieved from PNOA server(2015-2017) and time of the corresponding GEDI FW was registered (Dubayah et al. 2020b). The difference between GEDI99and ALS99_3Dwas calculated for each GEDI shot available after Filter #8. The assessed mean and standard deviation values for height increments at FMS polygon level were converted to annual height increments for better interpretation. The information for the time intervals is presented in Table 1.

    The sign of the difference was used to classified GEDI points as candidate metrics for harvests (negative) and natural dynamics (positive). The mean value and the standard deviation of height increments computed from ALS to GEDI was computed at FMS polygon level. The average values were compared to ground validation data and the proposed thresholds to cut off very large positive increments. The tendency and dispersion of the height increments was assessed in preliminary tests toward the effect of slope, the area of the FMS polygons and goodness of the DEM (Digital Elevation Model) derived from both laser sources considering ALS as the reference value. The results were visually displayed to confirm the capability of GEDI lasers to detect harvests (i.e., negative values higher than 5 m per year), fast growing dynamics(positive increments) while the transition region denoted as “uncertainty” was carefully inspected using satellite imagery to confirm the occurrence of cuttings and addressing the effect of intra-polygon variability in the FMS. The dataset was enlarged by adding polygon area,elevation, slope in degrees and the near distance from each GEDI point to the boundary of the FMS polygon(greater than 25 m after applying Filter #3). The geospatial operations were iteratively implemented in the R statistical software (R Core Team 2020). The captured laser-based increments were validated using an independent set of high-quality ground NFI observations.

    Validation data for ground elevation

    For the final DEM surface generation, the ALS ground points were converted to a TIN surface raster using 2-m raster DEM. A memory-efficient-streaming technology under three parallel processes was computed using the blast2dem command available in Lastools (Isenburg 2020). For ground ALS, the average of the 2-m raster grid cells was computed at 25-m diameter footprint level, while the elev_lowestmode attribute in GEDI Level 2A provided the ground elevation of GEDI laser shots.The vertical accuracy of ground elevation was evaluated using the differences between ground ALS and ground elevation GEDI (elev_lowestmode) within the 25-mdiameter footprint. The following statistics was calculated: correlation coefficient of Pearson (r), the Root Mean Squared Error (RMSEz, Eq. 1), mean of the differences referred to as biasz (Eq. 2), and relative bias(rbiasz) (Eq. 3).

    Table 1 Summary information for the three species on time intervals between ALS adquisition date and GEDI laser shooting computed for the set GEDI shots analysed

    Table 2 Summary of the 256 samples from National Forest Inventory of Spain (NFI) located in Lugo Province and used for validation

    where n is the number of GEDI shoots, yiis the mean elevation of ALS raster surface model for the 25-m diameter footprint GEDI shot i and ^yiis the ground elevation estimation from the 25-m footprint GEDI defined by the elev_lowestmode attribute in the GEDI Level 2A.

    The accuracy of ground elevation at footprint GEDI level was also tested via linear regression model using DEMALSas variable dependent as ground truth value(Eq. 1). To determine the accuracy of ground elevation at footprint GEDI the following statistics were computed from the model: model efficiency (adjR2, Eq. 2), the overall root mean square error (RMSE, Eq. 3), the relative root mean square error (rRMSE, Eq. 4).

    where n is the number of GEDI shoots as the total number of observations used to fit the linear regression,yiis the mean elevation of ALS raster surface model from the 25-m footprint GEDI shot i and ^yiis the estimated ground elevation from the fitted linear regression model using elev_lowestmode the GEDI Level 2A as dependent variable.

    Validation data for height growth dynamics

    The rotation of the NFI in the Northern and most Atlantic provinces of Spain is set to 5 years, half or a third of the normal rotation for a province (10-15 years). The last two rounds of the NFI in Galicia are denoted as 4th and 5th NFI, respectively. We used the permanent NFI samples located in the province of Lugo for the three species to compute the annual height increment using ground tree measurements. The Spanish NFI protocol for the tree measurements can be found in álvarez-González et al. (2014) and in previous research by the authors (Pascual et al. 2021; Guerra-Hernández et al.2021). The ground validation of heights increments for Eucalyptus globulus can be substantially biased with NFI data so we restricted the validation to P. pinaster and P.radiata using 198 and 60 plots, respectively. Data for the 4th NFI data for P. pinaster and P. radiata were measured between March and September 2009. The measurement of the 5th NFI was between May and October 2018. The difference in dominant height between the inventories ranged between 0.01 and 1.471 m per year in NFI samples for P. radiata, while the upper growth bound for P. pinaster samples ranged from 0.01 to 1.1 m. The maximum height increments were used to filter GEDI ground shots whose height dynamics (i.e.,ALSH99subtracted from GEDI99) was strictly positive.The maximum values for the ground data were validated using studies in Galicia for the species (Diéguez-Aranda et al. 2009, 2012). The validation data is presented in Table 2 and the flowchart of the proposed methodology is shown in Fig.4.

    Fig.4 Flowchart of the research showing the collection and integration of laser data from both airborne and GEDI satellite laser into the operational forest inventory to estimate height growth dynamics

    Results

    Comparison of ALS and FW lidar-derived ground elevation

    Fig.5 a Scatterplots of LF and ALS-derived DEM elevation at level of footprint GEDI from Filter #5(quality flag value=1,n=9875,Filter #5),b Boxplot of elevation differences between GEDI (elev_lowestmode)and DEMALS_2m considering terrain slope classes at GEDI footprint level.The upper and lower whiskers extend from to the highest and lowest value respectively within the 1.5 times the interquartile range

    The FW ground elevation and ALS lidar-derived DEM(2 m resolution) data was strongly correlated (r=0.99;n=9875, Fig. 5a). In terms of RMSEz, the vertical accuracy across study area was 4.48 m (rRMSEz=0.79%)using ALS data as ground truth value. Mean differences between FW-derived elevation (elev_lowestmode) and ALS-derived elevation (denoted as DEMALS_2m) was 0.18 m (biasz=0.18 m) at footprint GEDI level. FW ground elevation estimations were positively biased meaning FW predominantly underestimated ground elevation when compared with ALS-based DEM. Linear regression improved the accuracy of FW ground height estimations from ALS data in terms of RMSE (RMSE=0.23 m and rRMSE=0.041%). Although slope effects also increased the dispersion values of differences between FW and observed ALS ground elevation (Fig. 5b), we did not observe a clear patter and substantial improvements in our results at FMS polygon level adding a filter based on steep conditions (Fig. 5b).

    Filtering and detection of fast-growing areas and harvests

    The use of filters significantly constrained the pool of valid GEDI shots: after Filter #2 the total set comprised 56,235 shots:18,269 for Eucalyptus spp.,13,527 for P.radiata and 24,439 for P.pinaster.After Filter#8,the final pool of valid GEDI shots decreased by 92.3%, reaching a total of 3566 observations:1218 for Eucalyptus spp.(16,313 ha),1726 for P. radiata (16,710 ha) and 622 for P. pinaster (7445 ha).The array of FMS polygons reached 384 observations.

    The negative height increments computed from the 2015-2017 ALS surveys to the 2019 GEDI scanning showed similar patterns for the three cases (Fig. 6). Negative increments classified as “harvest” corresponded to tree elevations above 15 m in ALS and those cases corresponded to low height values(i.e.,early stages of stand development)in 2019 when computing GEDI shot information. The negative annual height threshold of 5 m was sufficient to detect homogenous FMS polygons entirely harvested(Fig. 6). The mean annual height increment calculated for FMS polygons was displayed using the individual GEDI ground beams to confirm the results in those cases where the harvest coincided with the 2017 orthophotomap available for Galicia in the PNOA repository.The spatial layout of harvested polygons was displayed for the three species(Fig.7).

    The positive height dynamics computed from ALS surveys to GEDI scanning at FMS polygon level were, on average, higher comparing the mean values but more similar when comparing median values, especially for P.radiata. The variability of annual growth estimates with the ALS-GEDI method was higher than NFI ground observations (Fig. 8). The median and maximum annual height increments in NFI plots were 0.56 (SD=0.37)and 1.47 m for P. radiata, while for P. pinaster the values were 0.44 (SD=0.24) and 1.11 m, respectively.Under the ALS-GEDI approach the values were 0.62(SD=0.42) and 2.93 m for P. radiata, and 0.61 (SD=0.63) and 2.55 m for P. pinaster, respectively. The values for Eucalyptus spp. were 0.69 and 2.55 m, respectively.The proportion of FMS polygons classified as “growth”decreased by a third when adding Filters #7 and #8.Same as for harvests, we present examples of FM polygons classified as “growth” (Fig. 9).

    Uncertainty and scaling when interpreting height growth dynamics

    Fig.6 Mean height difference between 99th height percentiles compute from GEDI shots and airborne laser scanning data. The mean values were computed using polygons in the Forest Map of Spain(FMS).The inside-polygon GEDI shots were classified considering the above-ground height derived from GEDI(left),from ALS data(right)and the sign of the height difference.a,b: Eucalyptus; c, d: P.pinaster; e,f: P.radiata

    The “uncertain” region as presented in Fig. 6 comprised the greater number of non-positive height increments. We provide a graphical description of the most common sources of error when calculating mean annual height increment. The most recent aerial image from PNOA supported the verification process (Fig. 10) regarding the following identified sources:

    1) Variability of forest structure withing FMS polygons, especially for Eucalyptus and P. radiata from intensive thinning,partial clear-cuts and irregular planting.

    2) No-adjustment of FMS polygons to harvesting operations and real forest edges.

    3) Erratic behavior of variable GEDI99due to the positioning error of shots in sparse forest areas showing presence of gaps and discontinuities between forest patches, which was often the case for P.pinaster stands and P.radiata.

    Discussion

    The combined use of ALS-based forest maps and the timely GEDI laser data to estimate growth dynamics certainly depended on the spatial co-registration between laser sources. The research focused on the ability to aggregate, scale and re-use existing ALS-based raster data for the assessment of fast height growth dynamics and intensive management at a scale never collected before.The multi-temporal assessment of height dynamics in the ALS-GEDI approach arose the strengths and limitations of data fusioning. It is well accepted that ALS-based forest inventory is solid (Maltamo et al.2014). Therefore, the challenge comes from the integration of GEDI laser data to support forest mapping missions. There are important caveats in the interpretation of this work regarding data filtering, spatial co-registration or the singularities of commercial forestry and forest ownership in the area. In the paper,we highlighted the geospatial operations and filters required to combine GEDI data and large-scale ALSbased mapping.

    Fig.7 Showcases of FMS polygons classified as harvests from the height dynamics between ALS(2015-2017)and GEDI scanning(2019).a,b Clear-cuts in Eucalyptus stands;c Unevenly distributed planting for P.pinaster; d Clear-cut for P.pinaster; e Intra-polygon variation on management decisions; f P.radiata thinned stand.The 2017 orthophotomap is the background

    Calibration of ground level

    The nominal accuracy of horizontal co-registration for GEDI observations is between 8 and 10 m(Luthcke et al. 2020) while in the vertical the accuracy is close to 50 cm (section 4.31 in Dubayah et al.2020b). The vertical accuracy from GEDI sensor (filter criteria quality flag=1) was 4.45 m in terms of RMSEz with a positive bias of 0.18 m (biasz)considering as ground truth the DEM derived from the ALS surveys. Our values were similar to what Hancock et al. (2019) reported (bias of less 0.22 m and RMSE of 5.7 m using ALS as baseline) and higher than values Silva et al. (2018) showed (RMSE=1.64 m, bias=0.29, n=2987, spatial resolution=25 m). The consistent pattern in the GEDI FW sensor for ground elevation was under-estimation, which likely impacted our calculations by means of annual non-positive increments as documented by Hancock et al. (2019):underestimation of upper RH metrics at low canopy cover (see Fig. 7d in Hancock et al. 2019). The integration of GEDI data will benefit from improved horizontal positioning accuracy and locally calibrated algorithms (Luthcke et al. 2000, 2005).

    Fig.8 Average annual height growth for Eucalyptus spp. (Ec),P.pinaster(Pp)and P.radiata(Pr)estimated for the ALS-GEDI approach. Mean(dot)and median(bar line)values were computed using GEDI shots positioning within the polygons of the Forest Map of Spain.The results for NFI data were computed using the dominant height values at plot-level in Lugo for the species in the 4th and 5th NFI

    The processing of ALS data is always oriented to rescale the elevation of many laser echoes to ground level.The way ground bare land was calculated differed from the ALS surveys to GEDI data. In the first case, the multiple multi-layered ALS echoes per laser pulse allowed the estimation of the DEM using multiple ground points,while for the case of GEDI the estimation of all laser statistics relied on one FW laser pulse per beam. The estimation of ground level from satellite LiDAR is challenging (Duncanson et al. 2020) and the process becomes more complex in sparse forest ecosystems as we observed for the case of sparse P. pinaster stands and low-canopy cover conditions. We used the description of the GEDI geolocation procedure (see Section 6 Luthcke et al. 2000) to interpret our analyses. As previous studies reported (Hyde et al. 2005; Silva et al. 2018),negative increments were attributable to the spatial configuration of canopy elements in footprints located across the transition areas between non-forest and forest areas and multi-layered forests (Fig. 11). The accuracy of ground finder is a significant source of error in FW laser scanning to estimate bare-ground elevation before calculating RH metrics (i.e., above-ground heights). Despite disregarding the addition of a slope filter based on the DTM assessment in this study, slope can affect ground and above-ground height estimates in other ecosystems covering wide altitudinal gradientes (Hofton et al. 2002;Silva et al. 2018).

    Height growth dynamics and the training area

    The use of Galicia as training area favored the detection of fast-growing dynamics and the occurrence of scattered disturbances. The time gap of 2-3 years between the ALS surveys in Galicia and first GEDI shots was enough to ensure height increments beyond the precision of the laser sensors (Luthcke et al. 2020). Harvests were detected and their evaluation could have been simplified by accounting for the absolute values computed for each GEDI shot. The distinction between harvests and GEDI shots classified as uncertain was according the annual threshold of 5 m. The forest area affected by cuttings was greater but heights increments were computed at FMS polygon-level, the intra-polygon variation of forest attributes and the large size of FMS minimized the detected total harvested area. A representative FMS polygon was a mosaic of small forest patches for which management, ownership and species composition could be diverse. Consequently, many of the FMS polygons allocated in the uncertain region comprised harvested area. One way to detect those conflictive cases is to evaluate the homogeneity of height annual increments within the boundaries.

    Fig.9 Areas classified as fast-growing when accounting for height dynamics using laser data.a,b Even aged and homogenous polygons are presented for Eucalyptus spp.;c, d P.pinaster stands;e mature P.radiata plantation;f early-stage P.radiata plantation

    The comparison of FMS polygons attending to the variability of the mean annual increments was a good proxy to detect variability across the landscape and detect those management units above a certain variability threshold on which harvests must have occurred above a certain height. The results showed, for instance, how large Eucalyptus plantations can be as homogeneous as small pine stands, and it makes sense if we consider a large plantation as the sum of small homogeneous forest units. The recognition of landscape fragmentation in the FMS was not as accurate as desirable, and that is the main reason for the extent of the “uncertain” level in our results(Fig. 12). The reasoning is supported by previous studies that documented the issue under multitemporal observations based on ALS surveys (Hopkinson et al. 2008).

    Fig.10 Showcases of non-positive height dynamics from ALS to GEDI statistics in Eucalyptus areas focusing on canopy structure.Low (blue)and dense canopy cover(yellow,CC>70%)cases are displayed with the corresponding both ALS point clouds and simulated GEDI full waveformusing gediWFSimulator command from rGEDI package(Hancock et al.2019)

    The results confirm that the accuracy of FW estimates of forest growth depends on the complexity of the horizontal and vertical structural diversity of the vegetation.From the earned experience and extensive visual inspection attending to the spatial layout of GEDI shots within FMS polygons, the presented approach is robust for homogeneous ecosystems on which treatment prescriptions have not impacted horizontal tree spacing. In this sense, the annual rates of lidar estimated growth from more homogeneous P. radiata plantations were close to(but slightly higher than) the observed field growth estimate for the 2009 to 2018 period of NFI plots. However,due to the P. pinaster forest structural diversity of Galicia (even-aged and uneven-aged forest), annual forest growth estimates were more variable than in homogeneous P. radiata plantations. The observed noise on sparse, naturalized, and more multi-layered P. pinaster forests could be explained by the estimation of bareground elevation from the GEDI full-waveform laser and the associated variability when computing ALS statistics in the presence of vegetation gaps (Hyde et al. 2005;Pascual et al.2019).The growth dynamics captured from GEDI to ALS surveying were slightly higher (median values) to the NFI plots used as the validation data, although the ground-based evidence reported in Galicia showed faster growth rates (Diéguez-Aranda et al. 2012).The impact of laser precision and co-registration requires the determination of some upper bounds for the height increments and, in this research, we benefit from the support of extensive studies on growth and yield modeling in the region to identify that upper bound to filter “growth” from “outliers”.

    Fig.11 Showcases of non-positive height dynamics in P.radiata polygons with GEDI footprints located in the proximity of forest edges,discontinuities and low canopy cover

    Fig.12 Showcases of areas classified as uncertainty when accounting for height dynamics between ALS and GEDI data and the spatial configuration of ground footprints in the transition areas between non-forest and in heterogenous patches in terms of canopy cover

    Co-registration of ALS and GEDI data

    The technical specifications of the GEDI laser system and official recommendations suggest considering between 8 and 10 m as the precision accuracy to determine the easting and northing the GEDI ground beams(Luthcke et al. 2020) (GEDI manual). The high degree of fragmentation in the landscape and the small size of forest assets conflicted with the nominal co-registration goodness of GEDI shots. Validation studies are being encouraged from the development team to improve the co-registration performance of the GEDI laser mission.The 30-m buffering distance towards the edge FMS polygons used in this study could have been higher but the impact on the computations might have been not as relevant as expected due to the heterogeneity of FMS polygons. However, by restricting the analyses to GEDI shots of >70% of forest cover and using the 99th height percentile to calculate height growth differences we certainly controlled the impact of edge effects. On the other hand, a slope filter was not used in our study. We did not observe substantial improvements in our results at Lugo province since the number of GEDI shots with high slopes values was low in our dataset. However, this type of filters could be restrictive in other forest areas topographically more complex with steep slopes. The reduction of 92.3% on the valid GEDI shots was important but necessary to isolate valid GEDI shots from the factors of uncertainty regarding positioning and ALS referencing. The constraint of data use is not relevant as the mission is collecting data on a weekly basis so the almost 300,000 GEDI shots initially evaluated could turn into 2 million by the end of 2020 while reducing the between-track distance. The reason for not using GEDI L3 products (i.e., gridded statistics on heights, forest cover or ecological indexes) was the high between-track distance for the small-scale and dispersed forest properties in the region. We evaluated the capability of GEDI to detect changes, but the accuracy of the measurements could be tested if ALS and GEDI were collected simultaneously. In future studies, Extremadura region (Spain)will be tested including the GEDI Level 3 as: i) the scale of forest management units is greater, ii) growth rates are lower and iii) forest management is not as intense as in commercial forestry in Galicia. The assimilation of GEDI data would benefit from parametrization. We hypothesize with the utility of allowing the user to rescale the estimates for upper RH percentiles based on accurate ALS-based ground information. The presented results and ideas can help GEDI scientists to improve and adjust the data processing mechanism.

    Conclusions

    The use of GEDI laser data provides valuable insights for forest industry operations especially when accounting for rapid spatio-temporal changes in forest structure high-quality ALS laser scanning data as baseline. Harvest and fast-growing forest are detected as well as the variability of forest attributes. The ALS-supported filtering methods described in this paper bring relevant insights to implement GEDI mission data, which is still operative. The correction of edges and the behavior of the laser signal in sparse forest conditions are important factor to consider and mitigate. Biome-specific algorithms to calibrate GEDI data science products with dense ALS data can upgrade the quality of the ground DEM and easy the assimilation of GEDI statistics under the novel and multi-temporal conceptualization of forest productivity.

    Acknowledgements

    We gratefully acknowledge Vicente Sandoval and Elena Robla from National Forest Inventory Department for supplying the inventory databases NFI data and the latest version of FMS. We also thanks to i) Dr. Carlos Alberto Silva(University of Maryland) for the relevant remarks, suggestions and his advice using rGEDI package and ii) Juan Gabriel álvarez González (Department of Agroforestry Engineering, University of Santiago de Compostela) for his support with NFI database.

    Authors’contributions

    The two authors were involved in all research stages from experiment conceptualization to scientific reporting. The author(s) read and approved the final manuscript.

    Funding

    This work was partially supported by‘National Programme for the Promotion of Talent and Its Employability’ of the Ministry of Economy, Industry, and Competitiveness (Torres-Quevedo program) via postdoctoral PTQ2018-010043 to Dr. Juan Guerra Hernández.The authors also thank to Forest Research Centre, a research unit funded by Funda??o para a Ciência e a Tecnologia I.P. (FCT), Portugal (UIDB/00239/2020). Arizona State University,USDA Forest Service and the Asner Lab supported Dr. Adrián Pascual in the final stages of the research.

    Availability of data and materials

    The research data and materials are available from the corresponding author on reasonable request.

    Ethics approval and consent to participate

    Not applicable.

    Consent for publication

    We express our consent to publish the paper in the journal.

    Competing interests

    The authors declare that they have no competing interests.

    Author details

    13edata,Centro de iniciativas empresariais, Fundación CEL. O Palomar s/n,27004 Lugo, Spain.2Forest Research Centre, School of Agriculture, University of Lisbon, Tapada da Ajuda, 1349-017 Lisbon, Portugal.3Center for Global Discovery and Conservation Science, Arizona State University, Hilo, HA 96720,USA.

    Received: 2 July 2020 Accepted: 1 February 2021

    美女国产视频在线观看| 国产成人精品婷婷| 久久韩国三级中文字幕| 又爽又黄无遮挡网站| 十八禁国产超污无遮挡网站| 成人亚洲精品一区在线观看 | 成人鲁丝片一二三区免费| 99热这里只有是精品50| 国产白丝娇喘喷水9色精品| 亚洲真实伦在线观看| 丰满人妻一区二区三区视频av| 久久6这里有精品| 国产乱来视频区| 国产探花在线观看一区二区| 我的老师免费观看完整版| 婷婷六月久久综合丁香| 91精品伊人久久大香线蕉| 日本黄色片子视频| 一级毛片aaaaaa免费看小| 欧美日韩综合久久久久久| 欧美日韩精品成人综合77777| 综合色av麻豆| 成年免费大片在线观看| 91精品一卡2卡3卡4卡| 国内少妇人妻偷人精品xxx网站| 亚洲18禁久久av| 嫩草影院新地址| 人妻系列 视频| av专区在线播放| 人妻夜夜爽99麻豆av| 美女国产视频在线观看| 午夜老司机福利剧场| 欧美 日韩 精品 国产| 日韩成人伦理影院| 亚洲欧美精品专区久久| 成人午夜高清在线视频| av线在线观看网站| 亚洲经典国产精华液单| 久久这里只有精品中国| 一个人观看的视频www高清免费观看| 久久精品久久久久久久性| 欧美变态另类bdsm刘玥| 91午夜精品亚洲一区二区三区| 黑人高潮一二区| 插阴视频在线观看视频| videos熟女内射| av天堂中文字幕网| 久久久久久久久久黄片| 熟女人妻精品中文字幕| 亚洲欧美成人精品一区二区| 久久国产乱子免费精品| av卡一久久| 精品久久久久久成人av| 日韩,欧美,国产一区二区三区| 国产午夜精品久久久久久一区二区三区| 综合色av麻豆| 女的被弄到高潮叫床怎么办| 看十八女毛片水多多多| 亚洲精品乱久久久久久| 久99久视频精品免费| 日韩一区二区视频免费看| 22中文网久久字幕| 91久久精品国产一区二区三区| 水蜜桃什么品种好| 国产探花极品一区二区| 欧美日韩在线观看h| 欧美性感艳星| 国产91av在线免费观看| 黄色日韩在线| 国产亚洲91精品色在线| 国产视频内射| 成年版毛片免费区| 精品午夜福利在线看| 亚洲最大成人手机在线| 日韩电影二区| 午夜福利高清视频| 国产精品精品国产色婷婷| 国产乱来视频区| 高清视频免费观看一区二区 | 免费看日本二区| 久久久久国产网址| 亚洲伊人久久精品综合| 中文乱码字字幕精品一区二区三区 | 国产熟女欧美一区二区| 777米奇影视久久| 亚洲欧美中文字幕日韩二区| 国产成人一区二区在线| 搡女人真爽免费视频火全软件| 毛片一级片免费看久久久久| 中文在线观看免费www的网站| 你懂的网址亚洲精品在线观看| 国产 亚洲一区二区三区 | 天天躁日日操中文字幕| 高清在线视频一区二区三区| 午夜福利高清视频| 一级二级三级毛片免费看| 人妻系列 视频| 国产色婷婷99| h日本视频在线播放| 日韩不卡一区二区三区视频在线| 久久久a久久爽久久v久久| 你懂的网址亚洲精品在线观看| 国产精品一区二区三区四区久久| 国产亚洲一区二区精品| 免费av不卡在线播放| 七月丁香在线播放| 欧美xxxx性猛交bbbb| 小蜜桃在线观看免费完整版高清| 深夜a级毛片| 两个人视频免费观看高清| 亚洲久久久久久中文字幕| 欧美最新免费一区二区三区| 日本猛色少妇xxxxx猛交久久| 日韩 亚洲 欧美在线| 18+在线观看网站| 久久国产乱子免费精品| 日本欧美国产在线视频| 国产伦理片在线播放av一区| 精品一区二区免费观看| 国产黄色视频一区二区在线观看| 免费大片黄手机在线观看| 国产不卡一卡二| 97热精品久久久久久| 精品久久久噜噜| 22中文网久久字幕| 欧美激情在线99| 老师上课跳d突然被开到最大视频| 国产免费福利视频在线观看| 国产黄色小视频在线观看| 精品99又大又爽又粗少妇毛片| av福利片在线观看| 久久精品综合一区二区三区| 精品久久久久久成人av| 男女啪啪激烈高潮av片| 色综合站精品国产| 久99久视频精品免费| 久久精品久久久久久久性| 一级二级三级毛片免费看| videossex国产| 日日啪夜夜爽| 你懂的网址亚洲精品在线观看| 色吧在线观看| 性色avwww在线观看| 又爽又黄a免费视频| 国产白丝娇喘喷水9色精品| 欧美bdsm另类| 亚洲aⅴ乱码一区二区在线播放| 99久国产av精品| 欧美不卡视频在线免费观看| 成人国产麻豆网| 一区二区三区乱码不卡18| 国产一区有黄有色的免费视频 | 高清日韩中文字幕在线| 能在线免费看毛片的网站| 国产亚洲av片在线观看秒播厂 | 热99在线观看视频| 欧美性感艳星| 日韩亚洲欧美综合| 成人二区视频| kizo精华| 国内少妇人妻偷人精品xxx网站| 丰满人妻一区二区三区视频av| 国产精品久久视频播放| 亚洲精品国产av蜜桃| 特级一级黄色大片| 午夜日本视频在线| 九色成人免费人妻av| 高清欧美精品videossex| 精品亚洲乱码少妇综合久久| 国内精品一区二区在线观看| av在线观看视频网站免费| 搡老乐熟女国产| 亚洲精品第二区| 久久久久国产网址| 久久久精品94久久精品| 精品人妻一区二区三区麻豆| 中文欧美无线码| 精品久久久精品久久久| 伦精品一区二区三区| 菩萨蛮人人尽说江南好唐韦庄| 在线免费十八禁| 欧美+日韩+精品| 国产精品国产三级国产av玫瑰| 少妇人妻精品综合一区二区| 九九在线视频观看精品| 日本爱情动作片www.在线观看| 国产色爽女视频免费观看| 国产极品天堂在线| 亚洲最大成人手机在线| 国产白丝娇喘喷水9色精品| 久久久国产一区二区| 蜜桃亚洲精品一区二区三区| 91午夜精品亚洲一区二区三区| 免费观看无遮挡的男女| 18+在线观看网站| 九九久久精品国产亚洲av麻豆| 国产淫片久久久久久久久| 内射极品少妇av片p| 亚洲一级一片aⅴ在线观看| 国产乱来视频区| 色网站视频免费| 国产亚洲最大av| 精品亚洲乱码少妇综合久久| 国产精品无大码| 免费观看的影片在线观看| 日韩成人伦理影院| 亚洲内射少妇av| 蜜桃亚洲精品一区二区三区| 亚洲精品视频女| 国产精品久久久久久av不卡| 欧美激情国产日韩精品一区| 色综合站精品国产| 成人美女网站在线观看视频| 精华霜和精华液先用哪个| 男人爽女人下面视频在线观看| 十八禁网站网址无遮挡 | 精品久久久久久久人妻蜜臀av| 日韩精品有码人妻一区| 亚洲欧美一区二区三区国产| 日本-黄色视频高清免费观看| 日本免费a在线| 欧美一级a爱片免费观看看| 91aial.com中文字幕在线观看| 亚洲va在线va天堂va国产| 国产又色又爽无遮挡免| 亚洲18禁久久av| 国产精品蜜桃在线观看| 日韩人妻高清精品专区| 欧美日韩国产mv在线观看视频 | 成人欧美大片| 免费电影在线观看免费观看| 欧美 日韩 精品 国产| 午夜福利网站1000一区二区三区| 一区二区三区乱码不卡18| 尾随美女入室| 国国产精品蜜臀av免费| 亚洲久久久久久中文字幕| 日韩av在线免费看完整版不卡| 丝袜美腿在线中文| 国产伦理片在线播放av一区| 国产成人免费观看mmmm| 亚洲av免费高清在线观看| 欧美+日韩+精品| 成人亚洲精品av一区二区| 国产爱豆传媒在线观看| 午夜精品在线福利| 精品国产露脸久久av麻豆 | 精品久久久久久久人妻蜜臀av| 国产精品久久久久久精品电影| 成年免费大片在线观看| 三级国产精品欧美在线观看| 一级毛片我不卡| 久久精品国产亚洲av涩爱| 国产av不卡久久| 久久精品国产亚洲av天美| 婷婷色麻豆天堂久久| 亚洲欧美清纯卡通| 国产高清不卡午夜福利| 91精品伊人久久大香线蕉| 国产亚洲精品av在线| av又黄又爽大尺度在线免费看| 女的被弄到高潮叫床怎么办| 中文欧美无线码| 亚洲欧美中文字幕日韩二区| 亚洲欧美一区二区三区黑人 | 日本熟妇午夜| 中文天堂在线官网| 国产欧美另类精品又又久久亚洲欧美| 一本一本综合久久| 少妇的逼水好多| 91av网一区二区| 在线观看免费高清a一片| 国产大屁股一区二区在线视频| av福利片在线观看| 国产一区亚洲一区在线观看| 久99久视频精品免费| 尾随美女入室| 在线观看一区二区三区| 大片免费播放器 马上看| av黄色大香蕉| 成人毛片a级毛片在线播放| 国产精品精品国产色婷婷| 国产一区有黄有色的免费视频 | 午夜激情福利司机影院| 一区二区三区高清视频在线| 亚洲最大成人手机在线| 肉色欧美久久久久久久蜜桃 | 直男gayav资源| 日韩欧美精品v在线| 最近2019中文字幕mv第一页| 亚洲精品aⅴ在线观看| 亚洲va在线va天堂va国产| 国产精品一区二区三区四区免费观看| 美女xxoo啪啪120秒动态图| 美女被艹到高潮喷水动态| 亚洲精品乱久久久久久| 天堂影院成人在线观看| 日韩制服骚丝袜av| 直男gayav资源| 亚洲精品成人久久久久久| 午夜激情福利司机影院| 亚洲精品亚洲一区二区| 国产精品一区www在线观看| 人人妻人人澡人人爽人人夜夜 | 午夜精品在线福利| 一级毛片aaaaaa免费看小| 亚洲电影在线观看av| 久久久久久国产a免费观看| 久99久视频精品免费| 18+在线观看网站| 午夜日本视频在线| 欧美成人一区二区免费高清观看| 伊人久久国产一区二区| 欧美三级亚洲精品| 久久久久久伊人网av| 久久精品夜色国产| 欧美激情久久久久久爽电影| 国产一区亚洲一区在线观看| 日本黄色片子视频| 国产高清不卡午夜福利| 床上黄色一级片| 中文字幕免费在线视频6| 身体一侧抽搐| 黑人高潮一二区| 一个人观看的视频www高清免费观看| 国产亚洲精品av在线| 菩萨蛮人人尽说江南好唐韦庄| 午夜激情欧美在线| av福利片在线观看| 午夜亚洲福利在线播放| 国产伦精品一区二区三区四那| 亚洲成人精品中文字幕电影| 91精品伊人久久大香线蕉| 免费黄网站久久成人精品| 婷婷色麻豆天堂久久| 婷婷色综合www| 国产精品一区二区三区四区免费观看| 丰满乱子伦码专区| 18禁动态无遮挡网站| 大陆偷拍与自拍| 国产黄a三级三级三级人| 女人被狂操c到高潮| 美女脱内裤让男人舔精品视频| 久久久精品94久久精品| 高清毛片免费看| 男女啪啪激烈高潮av片| 性插视频无遮挡在线免费观看| 日韩一区二区视频免费看| 少妇丰满av| 成人性生交大片免费视频hd| 免费看光身美女| 亚洲精品色激情综合| 国产亚洲一区二区精品| 精品国产露脸久久av麻豆 | 国产69精品久久久久777片| 汤姆久久久久久久影院中文字幕 | 国产免费福利视频在线观看| 国产淫片久久久久久久久| 寂寞人妻少妇视频99o| 免费观看的影片在线观看| 一级毛片黄色毛片免费观看视频| 99久国产av精品国产电影| 日本-黄色视频高清免费观看| 午夜免费激情av| 一级毛片黄色毛片免费观看视频| 国产精品熟女久久久久浪| 夜夜看夜夜爽夜夜摸| 婷婷色综合www| av.在线天堂| 春色校园在线视频观看| 国内精品美女久久久久久| 禁无遮挡网站| 亚洲精品aⅴ在线观看| 亚洲国产色片| freevideosex欧美| 欧美3d第一页| 欧美日韩精品成人综合77777| 久久久久久久大尺度免费视频| 美女主播在线视频| 99热这里只有精品一区| 久久久久久伊人网av| videos熟女内射| 亚洲一级一片aⅴ在线观看| 国产高清三级在线| 床上黄色一级片| 成人午夜精彩视频在线观看| 我要看日韩黄色一级片| 免费看光身美女| av在线播放精品| 少妇熟女aⅴ在线视频| 女人十人毛片免费观看3o分钟| 九草在线视频观看| 欧美日韩国产mv在线观看视频 | 成人亚洲精品一区在线观看 | 亚洲自偷自拍三级| 99久久精品一区二区三区| 成人特级av手机在线观看| 久久久欧美国产精品| 欧美97在线视频| 国产一区二区在线观看日韩| 直男gayav资源| 亚洲av.av天堂| 国产 一区 欧美 日韩| 日韩三级伦理在线观看| 久久久精品免费免费高清| 综合色av麻豆| 不卡视频在线观看欧美| 卡戴珊不雅视频在线播放| 黄色日韩在线| 久久久色成人| 亚洲精品乱码久久久久久按摩| 99热全是精品| 少妇被粗大猛烈的视频| 老师上课跳d突然被开到最大视频| 黑人高潮一二区| 日韩电影二区| 美女xxoo啪啪120秒动态图| 一个人看视频在线观看www免费| 久久97久久精品| 国产精品1区2区在线观看.| 性色avwww在线观看| 亚洲欧美精品自产自拍| 97精品久久久久久久久久精品| 91午夜精品亚洲一区二区三区| 亚洲精品第二区| 亚洲欧美精品专区久久| 亚洲精品乱码久久久v下载方式| 七月丁香在线播放| 久久这里有精品视频免费| 听说在线观看完整版免费高清| 伊人久久精品亚洲午夜| 亚洲高清免费不卡视频| 我的老师免费观看完整版| 久久精品夜色国产| av国产免费在线观看| 99久久人妻综合| 黑人高潮一二区| a级一级毛片免费在线观看| 欧美性感艳星| 国产精品一区二区在线观看99 | 黄色日韩在线| 成年人午夜在线观看视频 | 亚洲国产最新在线播放| 国产免费一级a男人的天堂| 非洲黑人性xxxx精品又粗又长| 深夜a级毛片| av线在线观看网站| 国产av不卡久久| 亚洲精品一二三| 搡老妇女老女人老熟妇| 成人亚洲欧美一区二区av| 午夜精品在线福利| 国产av在哪里看| 亚洲精品456在线播放app| 日本wwww免费看| 免费少妇av软件| 男女国产视频网站| 国产精品一区二区三区四区免费观看| 日本欧美国产在线视频| 永久免费av网站大全| 亚洲四区av| 欧美日韩综合久久久久久| 亚洲成色77777| 亚洲自偷自拍三级| 高清在线视频一区二区三区| 国内揄拍国产精品人妻在线| 久热久热在线精品观看| 国产黄频视频在线观看| 建设人人有责人人尽责人人享有的 | 国产永久视频网站| 能在线免费看毛片的网站| 人妻夜夜爽99麻豆av| a级毛色黄片| 亚洲美女视频黄频| 可以在线观看毛片的网站| 欧美一区二区亚洲| 国产老妇女一区| 亚洲国产色片| 免费观看a级毛片全部| 日韩欧美国产在线观看| 亚洲精品国产成人久久av| 免费大片黄手机在线观看| 高清欧美精品videossex| 亚洲精品日韩在线中文字幕| 欧美高清性xxxxhd video| 国产精品一区二区三区四区久久| 久久亚洲国产成人精品v| 91精品伊人久久大香线蕉| 麻豆久久精品国产亚洲av| 三级国产精品片| 男女下面进入的视频免费午夜| 国产亚洲av片在线观看秒播厂 | 精品久久久噜噜| 美女xxoo啪啪120秒动态图| 波野结衣二区三区在线| 国产精品美女特级片免费视频播放器| 美女大奶头视频| 免费黄频网站在线观看国产| 亚洲av电影在线观看一区二区三区 | 午夜福利在线观看免费完整高清在| 久久精品夜色国产| 日韩不卡一区二区三区视频在线| 黄色欧美视频在线观看| 搡老妇女老女人老熟妇| 寂寞人妻少妇视频99o| 久久精品久久久久久久性| 美女国产视频在线观看| 日韩在线高清观看一区二区三区| 亚洲精品影视一区二区三区av| 人妻夜夜爽99麻豆av| 亚洲成人一二三区av| 亚洲一级一片aⅴ在线观看| 日韩在线高清观看一区二区三区| 欧美xxxx黑人xx丫x性爽| 中文字幕av成人在线电影| 狠狠精品人妻久久久久久综合| 国产成人精品婷婷| 少妇熟女aⅴ在线视频| 亚洲av免费高清在线观看| av播播在线观看一区| 又爽又黄a免费视频| 国产精品爽爽va在线观看网站| 亚洲av二区三区四区| 日韩 亚洲 欧美在线| 免费播放大片免费观看视频在线观看| 成人欧美大片| 91精品伊人久久大香线蕉| 男的添女的下面高潮视频| 午夜福利成人在线免费观看| 午夜激情久久久久久久| 久久精品国产亚洲av涩爱| 亚洲四区av| 91av网一区二区| 黄色欧美视频在线观看| 一区二区三区四区激情视频| 纵有疾风起免费观看全集完整版 | 亚洲精品,欧美精品| 午夜久久久久精精品| 亚洲av男天堂| 黄色欧美视频在线观看| 一级a做视频免费观看| 精品人妻一区二区三区麻豆| 人人妻人人澡人人爽人人夜夜 | 亚洲美女搞黄在线观看| 国产 亚洲一区二区三区 | 国产日韩欧美在线精品| 色综合站精品国产| 国产成人freesex在线| 一区二区三区免费毛片| 国产一级毛片七仙女欲春2| 毛片女人毛片| 久久久久久九九精品二区国产| 三级国产精品片| 免费观看a级毛片全部| 精品久久久久久成人av| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 免费看日本二区| 精品国产三级普通话版| 国产一区有黄有色的免费视频 | 美女高潮的动态| 最新中文字幕久久久久| 久久久久久久久久久丰满| 在线免费观看的www视频| 免费少妇av软件| 久久久久久国产a免费观看| 黄色一级大片看看| 国产美女午夜福利| 亚洲av男天堂| 熟妇人妻久久中文字幕3abv| 麻豆成人午夜福利视频| 成人国产麻豆网| 欧美高清成人免费视频www| 亚洲va在线va天堂va国产| 精品人妻视频免费看| 精品一区二区三卡| 91久久精品国产一区二区成人| 亚洲激情五月婷婷啪啪| 中文字幕人妻熟人妻熟丝袜美| 国产精品女同一区二区软件| 麻豆久久精品国产亚洲av| 熟妇人妻不卡中文字幕| 国产高清国产精品国产三级 | 乱系列少妇在线播放| 亚洲怡红院男人天堂| 亚洲国产日韩欧美精品在线观看| 亚洲精品国产av成人精品| 97人妻精品一区二区三区麻豆| 成年女人在线观看亚洲视频 | 国产 一区 欧美 日韩| 69人妻影院| 男女边摸边吃奶| 真实男女啪啪啪动态图| 嫩草影院新地址| 男女啪啪激烈高潮av片| 丝袜喷水一区| 日日摸夜夜添夜夜爱| 久久久午夜欧美精品| 久久久久久久亚洲中文字幕| 久久综合国产亚洲精品| 中文天堂在线官网| 国产在视频线在精品| 日韩成人av中文字幕在线观看| 国产精品一区二区性色av| 日本三级黄在线观看| 欧美日韩视频高清一区二区三区二| 国产高潮美女av| 精品亚洲乱码少妇综合久久| 一个人免费在线观看电影| 边亲边吃奶的免费视频| 精品亚洲乱码少妇综合久久| 日韩成人av中文字幕在线观看| 在线观看av片永久免费下载| 观看免费一级毛片| 丝袜美腿在线中文| 2021天堂中文幕一二区在线观|