• <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

    免费av不卡在线播放| 日韩av在线免费看完整版不卡| 久久精品夜夜夜夜夜久久蜜豆| 欧美性猛交╳xxx乱大交人| 免费观看的影片在线观看| 久久精品国产99精品国产亚洲性色| 精品久久久噜噜| 久久久久久久久久成人| av播播在线观看一区| 黑人高潮一二区| 欧美成人a在线观看| 国产亚洲一区二区精品| 久久久欧美国产精品| 亚洲av二区三区四区| 国产 一区精品| 91久久精品国产一区二区三区| 午夜福利视频1000在线观看| 精品久久久久久久久久久久久| 天天一区二区日本电影三级| 色噜噜av男人的天堂激情| 国产亚洲午夜精品一区二区久久 | 国产伦理片在线播放av一区| 精品一区二区三区人妻视频| 又粗又硬又长又爽又黄的视频| 国语对白做爰xxxⅹ性视频网站| 老司机福利观看| 天堂av国产一区二区熟女人妻| 夜夜看夜夜爽夜夜摸| 搡女人真爽免费视频火全软件| 日本av手机在线免费观看| 国产精品一区www在线观看| 久久人妻av系列| 国产精品一区二区在线观看99 | 91aial.com中文字幕在线观看| or卡值多少钱| 久久精品国产99精品国产亚洲性色| 91精品一卡2卡3卡4卡| 一区二区三区高清视频在线| 七月丁香在线播放| 久久久久九九精品影院| videos熟女内射| 欧美最新免费一区二区三区| 免费电影在线观看免费观看| 国产亚洲av嫩草精品影院| 直男gayav资源| 黄色一级大片看看| 九九久久精品国产亚洲av麻豆| 女人十人毛片免费观看3o分钟| 免费看av在线观看网站| 熟妇人妻久久中文字幕3abv| 亚洲欧美日韩无卡精品| 91在线精品国自产拍蜜月| 一级黄片播放器| 麻豆一二三区av精品| 精品一区二区三区人妻视频| 国产一区二区在线av高清观看| 免费观看的影片在线观看| 亚洲国产精品国产精品| 不卡视频在线观看欧美| 自拍偷自拍亚洲精品老妇| 国产视频首页在线观看| 99热这里只有是精品在线观看| 欧美变态另类bdsm刘玥| 日本五十路高清| 综合色av麻豆| 国产乱来视频区| 免费一级毛片在线播放高清视频| 国产免费一级a男人的天堂| 色吧在线观看| 国产精品一区www在线观看| 精品国内亚洲2022精品成人| 欧美潮喷喷水| 又粗又硬又长又爽又黄的视频| 亚洲一级一片aⅴ在线观看| 日日撸夜夜添| 啦啦啦啦在线视频资源| 高清毛片免费看| 日韩欧美三级三区| 九色成人免费人妻av| 国产欧美另类精品又又久久亚洲欧美| 精品免费久久久久久久清纯| 九九热线精品视视频播放| 欧美成人午夜免费资源| 亚洲熟妇中文字幕五十中出| 少妇熟女aⅴ在线视频| 建设人人有责人人尽责人人享有的 | 日韩av在线大香蕉| 亚洲国产成人一精品久久久| 青春草国产在线视频| 男人和女人高潮做爰伦理| 蜜臀久久99精品久久宅男| 国产亚洲5aaaaa淫片| 日韩精品有码人妻一区| 又粗又硬又长又爽又黄的视频| 国产成年人精品一区二区| 国产黄a三级三级三级人| 国产av不卡久久| 久久久精品欧美日韩精品| 亚洲av电影在线观看一区二区三区 | 国产精品女同一区二区软件| 少妇猛男粗大的猛烈进出视频 | 亚洲精品乱码久久久久久按摩| 男人舔奶头视频| 美女大奶头视频| 久久久久久久久大av| av又黄又爽大尺度在线免费看 | 如何舔出高潮| 老女人水多毛片| 精品99又大又爽又粗少妇毛片| 欧美另类亚洲清纯唯美| 又黄又爽又刺激的免费视频.| 老师上课跳d突然被开到最大视频| 亚洲欧美清纯卡通| 91精品伊人久久大香线蕉| 97在线视频观看| 爱豆传媒免费全集在线观看| 亚洲国产精品成人综合色| 久久草成人影院| 两性午夜刺激爽爽歪歪视频在线观看| 黄片wwwwww| 欧美一区二区亚洲| 极品教师在线视频| 欧美日韩综合久久久久久| 国产精品久久久久久精品电影| 91精品伊人久久大香线蕉| 久久精品国产亚洲av天美| 精品酒店卫生间| 亚洲成人中文字幕在线播放| 精品久久久久久久末码| 中文字幕人妻熟人妻熟丝袜美| 亚洲成色77777| 在线观看一区二区三区| 亚洲精品国产av成人精品| av黄色大香蕉| 国内少妇人妻偷人精品xxx网站| 黑人高潮一二区| 最近手机中文字幕大全| 久久精品熟女亚洲av麻豆精品 | 国产三级在线视频| 久久6这里有精品| 亚洲美女搞黄在线观看| 美女脱内裤让男人舔精品视频| 亚洲一级一片aⅴ在线观看| 99久久精品国产国产毛片| 夜夜爽夜夜爽视频| 久久精品人妻少妇| av视频在线观看入口| 日本三级黄在线观看| 美女脱内裤让男人舔精品视频| 国产高清视频在线观看网站| 麻豆精品久久久久久蜜桃| 国产伦在线观看视频一区| 国产真实伦视频高清在线观看| 69人妻影院| 女的被弄到高潮叫床怎么办| 免费观看在线日韩| 听说在线观看完整版免费高清| 美女高潮的动态| 黄色一级大片看看| 国产欧美另类精品又又久久亚洲欧美| 色视频www国产| www.av在线官网国产| 久久久成人免费电影| 亚洲欧美清纯卡通| 国产伦理片在线播放av一区| 国产伦精品一区二区三区视频9| 一区二区三区乱码不卡18| 天堂网av新在线| 国内少妇人妻偷人精品xxx网站| 欧美性感艳星| 十八禁国产超污无遮挡网站| 国产精品国产三级专区第一集| 超碰av人人做人人爽久久| 欧美极品一区二区三区四区| ponron亚洲| av在线天堂中文字幕| 中文字幕亚洲精品专区| 成年女人看的毛片在线观看| 舔av片在线| 特大巨黑吊av在线直播| 在线天堂最新版资源| 国产亚洲91精品色在线| 三级国产精品片| 亚洲最大成人中文| 国产一区亚洲一区在线观看| 麻豆av噜噜一区二区三区| 亚洲经典国产精华液单| 97超视频在线观看视频| 国产精品国产三级国产专区5o | 禁无遮挡网站| 国产免费男女视频| 国产精品,欧美在线| 精品人妻熟女av久视频| 少妇高潮的动态图| 亚洲自偷自拍三级| 美女高潮的动态| 在线观看66精品国产| 国产av不卡久久| 国产视频内射| 在线天堂最新版资源| 成年女人永久免费观看视频| 久久精品国产鲁丝片午夜精品| 国产高潮美女av| 99九九线精品视频在线观看视频| 久久久久久久亚洲中文字幕| 国产精品国产三级专区第一集| 亚洲五月天丁香| 99久久精品一区二区三区| www日本黄色视频网| 最近手机中文字幕大全| 伊人久久精品亚洲午夜| 天堂√8在线中文| 特大巨黑吊av在线直播| 天天躁夜夜躁狠狠久久av| 亚洲aⅴ乱码一区二区在线播放| 久久久久久久亚洲中文字幕| 亚洲欧洲国产日韩| 国产黄色视频一区二区在线观看 | 午夜福利视频1000在线观看| 国产精品不卡视频一区二区| 九色成人免费人妻av| 国产精品久久久久久久久免| 男人舔奶头视频| 欧美成人一区二区免费高清观看| 国产免费男女视频| 变态另类丝袜制服| 久久久久久久久久久丰满| 日韩欧美 国产精品| 超碰av人人做人人爽久久| 中文在线观看免费www的网站| 2021天堂中文幕一二区在线观| 岛国在线免费视频观看| 亚洲国产成人一精品久久久| 91久久精品电影网| 少妇的逼水好多| 国产亚洲精品av在线| 天天躁夜夜躁狠狠久久av| a级毛片免费高清观看在线播放| 精品一区二区三区视频在线| 亚洲国产精品久久男人天堂| 黄色一级大片看看| 国产精品,欧美在线| 熟女人妻精品中文字幕| 久久亚洲国产成人精品v| 午夜日本视频在线| 国产私拍福利视频在线观看| 美女cb高潮喷水在线观看| 国产精品美女特级片免费视频播放器| 啦啦啦观看免费观看视频高清| 欧美最新免费一区二区三区| 日本色播在线视频| 偷拍熟女少妇极品色| 国产综合懂色| 美女大奶头视频| 1024手机看黄色片| 欧美区成人在线视频| 成人漫画全彩无遮挡| 国产久久久一区二区三区| 日韩av不卡免费在线播放| 毛片一级片免费看久久久久| 看免费成人av毛片| 亚洲精品国产成人久久av| 国产综合懂色| 久久精品夜夜夜夜夜久久蜜豆| 欧美xxxx性猛交bbbb| av天堂中文字幕网| 精品免费久久久久久久清纯| 国产伦理片在线播放av一区| 九色成人免费人妻av| av黄色大香蕉| 亚洲成人中文字幕在线播放| 国产中年淑女户外野战色| 精品不卡国产一区二区三区| 免费大片18禁| 人人妻人人看人人澡| 日日撸夜夜添| 亚洲久久久久久中文字幕| 日韩av不卡免费在线播放| 一本久久精品| 国产高清不卡午夜福利| 天天躁日日操中文字幕| 两性午夜刺激爽爽歪歪视频在线观看| 欧美一区二区亚洲| 大又大粗又爽又黄少妇毛片口| 26uuu在线亚洲综合色| 人妻少妇偷人精品九色| 全区人妻精品视频| 一区二区三区免费毛片| 亚洲,欧美,日韩| 变态另类丝袜制服| 男女那种视频在线观看| 97超视频在线观看视频| 99久久人妻综合| 亚洲乱码一区二区免费版| 2021少妇久久久久久久久久久| 色综合色国产| 午夜亚洲福利在线播放| 亚洲在线观看片| 久久久久久久久大av| 免费看日本二区| 亚洲精品日韩在线中文字幕| 精品欧美国产一区二区三| 久久久久久久久久久免费av| 久久精品久久精品一区二区三区| 18+在线观看网站| 少妇的逼好多水| 高清毛片免费看| av播播在线观看一区| av女优亚洲男人天堂| 欧美最新免费一区二区三区| 国产精品久久久久久av不卡| 搞女人的毛片| 一级二级三级毛片免费看| 日韩视频在线欧美| 97超视频在线观看视频| 国产精品麻豆人妻色哟哟久久 | 亚洲熟妇中文字幕五十中出| 精品久久久噜噜| 久久99热这里只有精品18| 国产精品伦人一区二区| 欧美日本亚洲视频在线播放| 麻豆精品久久久久久蜜桃| 中文字幕av成人在线电影| 日本欧美国产在线视频| 九九在线视频观看精品| 18禁在线无遮挡免费观看视频| 日本熟妇午夜| 一区二区三区四区激情视频| 精品少妇黑人巨大在线播放 | 三级男女做爰猛烈吃奶摸视频| 日本爱情动作片www.在线观看| 国产午夜精品一二区理论片| 国产视频首页在线观看| 亚洲国产日韩欧美精品在线观看| 欧美性感艳星| 伦精品一区二区三区| 男女国产视频网站| 黄色欧美视频在线观看| 天堂√8在线中文| 久久久精品大字幕| 国产成年人精品一区二区| 中文字幕人妻熟人妻熟丝袜美| 热99在线观看视频| 精品99又大又爽又粗少妇毛片| 亚洲欧美精品专区久久| 成年免费大片在线观看| 欧美又色又爽又黄视频| 精品无人区乱码1区二区| 国产高潮美女av| 精品欧美国产一区二区三| 人体艺术视频欧美日本| 一区二区三区乱码不卡18| 在线a可以看的网站| 99国产精品一区二区蜜桃av| 日本熟妇午夜| 色播亚洲综合网| 色综合色国产| 欧美激情在线99| 国产成人freesex在线| 免费看美女性在线毛片视频| 欧美+日韩+精品| 嫩草影院精品99| 国产成人午夜福利电影在线观看| 亚洲精品国产成人久久av| 看片在线看免费视频| 久久人人爽人人爽人人片va| 成人欧美大片| 淫秽高清视频在线观看| 色5月婷婷丁香| 久久久精品欧美日韩精品| 人人妻人人看人人澡| 久久久精品94久久精品| 可以在线观看毛片的网站| 国产成人精品婷婷| 久久精品熟女亚洲av麻豆精品 | 欧美性猛交黑人性爽| 一卡2卡三卡四卡精品乱码亚洲| 国产亚洲av片在线观看秒播厂 | 深爱激情五月婷婷| 色视频www国产| 久久精品国产亚洲av涩爱| 22中文网久久字幕| 舔av片在线| 又爽又黄a免费视频| 日本黄大片高清| 国产黄片视频在线免费观看| 国产 一区 欧美 日韩| 亚洲av中文av极速乱| 综合色丁香网| 国产av一区在线观看免费| 亚洲精品乱码久久久久久按摩| av视频在线观看入口| 黄色欧美视频在线观看| 日本免费一区二区三区高清不卡| 亚洲国产精品成人综合色| 天堂av国产一区二区熟女人妻| 免费播放大片免费观看视频在线观看 | 亚洲国产精品成人久久小说| 久久久国产成人免费| 最近中文字幕2019免费版| 女人十人毛片免费观看3o分钟| 秋霞伦理黄片| 久久精品久久久久久噜噜老黄 | 人体艺术视频欧美日本| 国产成人a∨麻豆精品| 国产探花极品一区二区| 精品久久国产蜜桃| 国产91av在线免费观看| 国产免费一级a男人的天堂| 嫩草影院入口| 日日摸夜夜添夜夜添av毛片| 亚洲国产精品sss在线观看| 丝袜美腿在线中文| 亚洲丝袜综合中文字幕| 亚洲不卡免费看| 午夜福利高清视频| 噜噜噜噜噜久久久久久91| 午夜久久久久精精品| 精品人妻视频免费看| 美女cb高潮喷水在线观看| 亚洲欧美成人精品一区二区| 国产精品久久久久久精品电影| 国产乱来视频区| 国语自产精品视频在线第100页| 狂野欧美激情性xxxx在线观看| 午夜免费男女啪啪视频观看| 精品免费久久久久久久清纯| av福利片在线观看| 国产淫片久久久久久久久| 精品人妻偷拍中文字幕| 国产精品人妻久久久影院| 亚洲av成人精品一区久久| 少妇的逼好多水| 国产午夜精品一二区理论片| 日日干狠狠操夜夜爽| 免费无遮挡裸体视频| 欧美成人一区二区免费高清观看| av在线蜜桃| 非洲黑人性xxxx精品又粗又长| 国产av在哪里看| 欧美成人精品欧美一级黄| 九九久久精品国产亚洲av麻豆| 国产亚洲av片在线观看秒播厂 | 中国国产av一级| 欧美日韩国产亚洲二区| 一二三四中文在线观看免费高清| 国产精品久久久久久久电影| 成人性生交大片免费视频hd| 婷婷色麻豆天堂久久 | 高清日韩中文字幕在线| 日韩视频在线欧美| 偷拍熟女少妇极品色| 免费无遮挡裸体视频| 免费观看的影片在线观看| 久久热精品热| 欧美成人免费av一区二区三区| 免费av毛片视频| 欧美潮喷喷水| 搡老妇女老女人老熟妇| 一级黄色大片毛片| 精品一区二区三区视频在线| 欧美日韩在线观看h| 亚洲国产精品国产精品| 国产av不卡久久| 久久久久免费精品人妻一区二区| 亚洲高清免费不卡视频| 国产精品综合久久久久久久免费| 女人被狂操c到高潮| 热99在线观看视频| 啦啦啦韩国在线观看视频| 日韩强制内射视频| 少妇高潮的动态图| 久久久色成人| 91aial.com中文字幕在线观看| 高清毛片免费看| 中文精品一卡2卡3卡4更新| 91精品一卡2卡3卡4卡| 日日摸夜夜添夜夜爱| 亚洲中文字幕日韩| 日韩一区二区三区影片| 99久久精品热视频| 久久这里有精品视频免费| 午夜亚洲福利在线播放| 欧美日韩一区二区视频在线观看视频在线 | 亚洲精品自拍成人| 精品久久久久久成人av| 国产成人精品婷婷| 久久精品国产亚洲av涩爱| 18禁裸乳无遮挡免费网站照片| 精品人妻偷拍中文字幕| 欧美另类亚洲清纯唯美| 变态另类丝袜制服| 久久久久国产网址| 毛片女人毛片| 日产精品乱码卡一卡2卡三| 午夜精品国产一区二区电影 | 欧美一区二区精品小视频在线| 亚洲一区高清亚洲精品| 国产成人91sexporn| 麻豆一二三区av精品| 岛国在线免费视频观看| av在线蜜桃| 国产精品伦人一区二区| 国产精品人妻久久久影院| 亚洲无线观看免费| 国模一区二区三区四区视频| 亚洲精品亚洲一区二区| av线在线观看网站| 变态另类丝袜制服| 大香蕉久久网| 国产人妻一区二区三区在| 婷婷色综合大香蕉| 夜夜爽夜夜爽视频| 91在线精品国自产拍蜜月| 老司机影院毛片| 最近2019中文字幕mv第一页| 最新中文字幕久久久久| 精品国内亚洲2022精品成人| 久久6这里有精品| 午夜福利网站1000一区二区三区| 如何舔出高潮| 亚洲精品色激情综合| 国产午夜精品论理片| 国模一区二区三区四区视频| 汤姆久久久久久久影院中文字幕 | 自拍偷自拍亚洲精品老妇| 欧美成人午夜免费资源| 午夜日本视频在线| 女的被弄到高潮叫床怎么办| 中文精品一卡2卡3卡4更新| 伊人久久精品亚洲午夜| 久久久久网色| 91午夜精品亚洲一区二区三区| 午夜福利成人在线免费观看| 国产黄a三级三级三级人| 亚洲国产日韩欧美精品在线观看| 久久久久九九精品影院| 三级国产精品欧美在线观看| 女人被狂操c到高潮| 国产一区亚洲一区在线观看| 欧美变态另类bdsm刘玥| 丝袜美腿在线中文| 午夜福利高清视频| 色视频www国产| 亚洲av免费在线观看| 精品国内亚洲2022精品成人| 久久6这里有精品| ponron亚洲| 成人鲁丝片一二三区免费| 最近最新中文字幕免费大全7| 午夜视频国产福利| 美女内射精品一级片tv| 日日啪夜夜撸| 高清毛片免费看| or卡值多少钱| 日本黄色视频三级网站网址| 狂野欧美白嫩少妇大欣赏| 国产极品天堂在线| 一区二区三区免费毛片| 大又大粗又爽又黄少妇毛片口| 精品国产露脸久久av麻豆 | 不卡视频在线观看欧美| 校园人妻丝袜中文字幕| 国产亚洲精品av在线| 中文字幕熟女人妻在线| 高清在线视频一区二区三区 | 亚洲自偷自拍三级| 亚洲在线自拍视频| 五月玫瑰六月丁香| av在线老鸭窝| 精品久久久久久久久av| 99热这里只有是精品在线观看| 啦啦啦啦在线视频资源| 成年女人看的毛片在线观看| 久久婷婷人人爽人人干人人爱| 国产精品国产三级国产专区5o | 亚洲精品成人久久久久久| 亚洲精品一区蜜桃| 久久久成人免费电影| 精品久久久久久久久久久久久| 天堂av国产一区二区熟女人妻| 美女国产视频在线观看| 亚洲国产精品久久男人天堂| 精品久久久久久久久亚洲| a级毛色黄片| 小说图片视频综合网站| 天天躁日日操中文字幕| 国产黄片美女视频| 六月丁香七月| 亚洲在久久综合| 视频中文字幕在线观看| 国产黄片美女视频| 男人和女人高潮做爰伦理| 三级国产精品欧美在线观看| 精品无人区乱码1区二区| 午夜日本视频在线| 蜜桃亚洲精品一区二区三区| av免费观看日本| 尤物成人国产欧美一区二区三区| 久久婷婷人人爽人人干人人爱| 国产免费又黄又爽又色| 中文字幕熟女人妻在线| 中文字幕av在线有码专区| 在线观看美女被高潮喷水网站| 欧美日韩国产亚洲二区| 嫩草影院入口| 亚洲精品aⅴ在线观看| 精品人妻熟女av久视频| 亚洲乱码一区二区免费版| 久久精品夜色国产| 国产私拍福利视频在线观看| 国内精品一区二区在线观看| 中文乱码字字幕精品一区二区三区 | 欧美+日韩+精品| 国产美女午夜福利|