999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Improving precision of field inventory estimation of aboveground biomass through an alternative view on plot biomass

2021-01-11 07:06:50ChristophKleinnSteenMagnussenNilslkePaulMagdonJuanGabriellvarezGonzlezLutzFehrmannandsarrezCruzado
Forest Ecosystems 2020年4期

Christoph Kleinn ,Steen Magnussen,Nils N?lke,Paul Magdon,Juan Gabriel álvarez-González,Lutz Fehrmann and César Pérez-Cruzado,

Abstract We contrast a new continuous approach (CA) for estimating plot-level above-ground biomass (AGB) in forest inventories with the current approach of estimating AGB exclusively from the tree-level AGB predicted for each tree in a plot, henceforth called DA (discrete approach). With the CA, the AGB in a forest is modelled as a continuous surface and the AGB estimate for a fixed-area plot is computed as the integral of the AGB surface taken over the plot area. Hence with the CA,the portion of the biomass of in-plot trees that extends across the plot perimeter is ignored while the biomass from trees outside of the plot reaching inside the plot is added. We use a sampling simulation with data from a fully mapped two hectare area to illustrate that important differences in plot-level AGB estimates can emerge. Ideally CA-based estimates of mean AGB should be less variable than those derived from the DA. If realized, this difference translates to a higher precision from field sampling, or a lower required sample size. In our case study with a target precision of 5%(i.e. relative standard error of the estimated mean AGB), the CA required a 27.1% lower sample size for small plots of 100 m2 and a 10.4%lower sample size for larger plots of 1700 m2.We examined sampling induced errors only and did not yet consider model errors. We discuss practical issues in implementing the CA in field inventories and the potential in applications that model biomass with remote sensing data. The CA is a variation on a plot design for above-ground forest biomass; as such it can be applied in combination with any forest inventory sampling design.

Keywords: Branch biomass, Foliage biomass, Stem biomass, Biomass surface plots, Sampling surfaces, Standard error of estimation

Background

Fixed-area field plots have a long history in forest inventory and they are commonly used in national forest inventories. In the early nineteenth century or even earlier,fixed-area plots were used to expand values recorded per-plot to entire stands (e.g. K?nig 1835; Heyer 1861).That was long before statistical sampling had been established, called at its time the “representative method” (Kiaer 1895-1896).

When sampling for tree attributes, there are at least two different views on fixed-area plots. Each carries a different analytical approach: (1) The plot cuts out a sample area from a population (defined as the horizontally projected surface of a defined area of forest land) with recordings of all variables exactly above that area;recordings may include individual objects but also continuously distributed variables. Henceforth, we call this approach the “continuous approach to plot-level biomass”, CA. (2) The plot serves as a selection tool and defines the discrete set of sample trees to be included into this particular plot from the population of all trees in a defined area of forest land.We will call this approach the“discrete approach to plot-level biomass”,DA.In the CA, the per-plot value is the sum of the target variable exactly above the horizontal projection of the plot area(for circular plots,the portion of the target variable within a cylinder on the plot perimeter); in the DA the per-plot value comes from the sum of the target variable over all included sample objects(i.e.trees).

In remote sensing applications, when it is possible to match remotely sensed data used as predictors of biomass to the specific field plot location, the CA has intuitive appeal because the plot area is matched as accurately as possible with the image data (i.e. pixels).N?sset (2002) named this technique the “area-based approach” (ABA), which, in remote sensing is contrasted to the “individual tree detection” approach. In field sampling, the DA is commonly applied, in which the tree variables for all sample trees within the plot are recorded, regardless of whether a sample tree has a portion of its stem and crown outside of the plot area.

Both approaches, the CA and the DA, and their corresponding analytical approaches have been applied in plot design comparisons when sampling for dead wood: Gove and Van Deusen (2011) introduced the instructive terms“chainsaw method”, when only those parts of the deadwood that are exactly within the plot perimeter are recorded (as if cut out by a chainsaw; this corresponds to our CA), and “stand-up method”, when all those deadwood pieces with their thicker end within the plot area are included and recorded (as if these dead wood pieces would stand up within the plot; this corresponds to our DA). In forest inventory field sampling, both approaches may be applied. However, they differ in the required efforts and the precision of estimation; design-unbiased estimators are available for both approaches.

To the best of our knowledge, the DA and the CA have not been compared for their influence on AGB estimation, which has become a core variable in forest monitoring. The role of forests in global carbon cycling accentuates the importance of quantifying both the status and trends of forest biomass (IPCC 2006) and points to the demand for “credible” estimates.

Above-ground biomass is defined as “all living plant biomass above the soil … including stem, stump, branches,bark,seeds and foliage”(IPCC 2006,Annex 4A.1,4.72).To bypass issues of destructive biomass sampling,tree biomass is predicted with models linking biomass to easy-tomeasure tree variables,namely the diameter at breast height(dbh) given its correlation with stem basal area, and tree height (Brown et al. 1989; Fehrmann et al. 2008; Picard et al. 2012; Chave et al. 2015; Magnussen and Carillo Negrete 2015).

In this study, we compare the precision of estimation of AGB from fixed-area field inventory plots under the two approaches CA and DA. The DA is the approach detailed in most forest inventory textbooks in which plot-level AGB is predicted as the sum over the sample trees in a plot of the biomass predicted for each tree using an appropriate model (Kershaw et al. 2016). A DA estimate of plot-level biomass may include crown parts of sample trees that extend outside of the plot perimeter and biomass in oblique stems leaning outside of the plot area.

We consider the CA with tree biomass viewed as a variable that is continuously distributed horizontally. Its application requires, for each inventory plot, the prediction of the plot-level biomass to include all parts with a horizontal projection within the area defined by the plot.A continuous view is also adopted in forest monitoring when defining an infinite population of sample elements(i.e. plots) in an area sampling frame (Mandallaz 1991),but applying a continuous view to the inventory variable biomass is new and modelling the continuous horizontal biomass distribution (HBD) to produce plot-level predictions of biomass has not been reported before. Aside from forest field inventories, a continuous view is not uncommon, for example, for topographic variables such as height above sea level, and remotely-sensed canopy height models.

Figure 1 illustrates the key difference between the DA and CA for the estimation of tree variables.

While the vertical distribution of forest biomass has been researched intensively (e.g. Tahvanainen and Forss 2008; Ruiz-González and álvarez-González 2011; Nemec et al. 2012; Jiménez et al. 2013), studies that address the horizontal forest or tree biomass distribution are sparse.Kershaw and Maguire (1996) and Xu and Harrington(1998) modelled the horizontal distribution of leaf biomass - but not AGB. Affleck et al. (2012) noted an absence of HBDs for the modelling of fuel load where HBDs would be required to fully utilize next-generation fire behaviour simulators, and to improve fuel management strategies to reduce fire hazard. Mascaro et al.(2011) considered a uniform HBD across the crown projection area in a study on uncertainty in carbon estimation. They may have been the first to integrate the DA/CA issue in the context of airborne laser scanning (ALS)for forest carbon. They named the approaches “stem-localized” (corresponding to our DA and Gove and Van Deusen’s (2011) “stand-up” approach) and “crown-distributed” (corresponding to our CA and Gove and Van Deusen’s (2011) “chainsaw” method). Mascaro et al.(2011) investigated the issue in a remote-sensing based modelling context. Their main finding was that the root mean square errors (RMSE) were lower for the “crowndistributed” approach when modelling biomass from ALS data. As expected, the results varied with field plot size. Packalen et al. (2015) suggested a different approach to the edge tree issue in the ABA in ALS-based modelling; they modified the plot area so that the crowns of in-trees were fully contained in an extended plot while overlapping crown parts of out-trees were discounted. This technique allowed them to compute the predictor variables tree-wise, while area-based predictions were still made for the original plot.

This study is a first step towards integrating the horizontal distribution of forest biomass into AGB estimation from forest inventory field sampling. To that end,we provide a case study where DA estimates of AGB are compared to CA estimates predicted by considering AGB as a continuous variable across an area of interest.Our case study considers the sampling induced standard error exclusively; we are not yet dealing here with a fully realized error propagation including measurement and model errors. The findings of our case allow us to draw first conclusions about the CA’s statistical performance,and identify pathways for further developments.

Materials and methods

Study area

Tree heights and maximum crown diameters(CDmax= diameter of crown projection area) were predicted from dbh via local models (Guericke 2001). A circular tree crown projection was assumed.

Construction of a tree-level HBD

The cornerstone in the CA is tree-specific HBDs for all trees in a plot, and for trees outside of the plot whose branches and foliage extend inside the plot perimeter.The HBD allocates a predicted tree biomass over a tree’s crown projection area.

Mascaro et al. (2011) proposed a uniform distribution by assigning (simply but quite unrealistically) the same biomass density to the inner and outer parts of the tree crown. To more closely reflect reality, we developed a per-tree crown horizontal biomass distribution by combining a simple stem biomass model with both branchwood models and foliage biomass models, derived from allometry and horizontal foliage distribution.

As an approximation, we applied one single general HBD model for all trees.

In the first step, we predicted the AGB for each tree separately for the stem, the branches and foliage using models provided in Bartelink (1997). In the second step,we distributed this AGB horizontally along our model of HBD over the crown projection area.

The compartment-wise development of our HBD is given in Additional file 1 (Max and Burkhart 1976).Figure 2 illustrates this HBD by giving a 2D profile of this three-dimensional distribution.

The 3D HBD for each biomass component (stem,branches and foliage) was generated by rotating the corresponding functions around the stem axis. We standardized the X-axis (distance from the tree axis in Fig. 2)of the distribution to the CDmax/2 of each tree, and distributed the volume under the rotation solid over thehorizontal plane. Hence, we assumed for the time being that the HBD is isotropic in the horizontal plane,and that all trees have the same standardized distribution. Results of an as of yet unpublished empirical study on beech trees confirm the general shape of this distribution.

Table 1 Characteristics of the 2 ha study area.N=stand density,G=basal area,and dg=quadratic mean diameter

Having built the model for the individual tree HBD,the basic difference between the two approaches becomes obvious: in the DA, the total tree predicted AGB is located at the stem position whereas in the CA the same total tree AGB is distributed over the assumed circular crown projection area along the HBD model. That means that the considered and evaluated AGB per tree remains the same. The differences between the DA and the CA are (1) how this AGB is horizontally distributed and (2) which portion of this AGB is subsequently recorded for a defined plot area.

Construction of an area-level HBD

To illustrate the difference between AGB estimates from the DA and the CA for an area of interest, we needed an HBD model for the entire area (here: our 2-ha study area). In principle this is done by predicting the HBD for each tree and summing the predicted biomass for each location in the study area. For practicality, we chose to rasterize the study area into cells of 2 cm×2 cm, resulting in a total of 50×106cells. This fine spatial resolution of the horizontal distribution of biomass over the study area provided an approximation to the continuous biomass distribution.

To facilitate a comparison between DA and CA estimates of AGB, and to better illustrate their differences,we also produced an HBD model for the study area as it would emerge with the DA. In contrast to the CA, the DA specific HBD model assigns the predicted AGB of a tree to a single 2 cm×2 cm cell containing the stem position; resulting in an HBD map where the individual tree biomass is concentrated at the tree positions and displays as peaks whose height depends on the respective biomass.

The DA and CA area-level HBD models were then used in a sample simulation as a basis for comparison.

Sample simulation

To capture effects of plot-size and shape, we carried out simulations with two common plot shapes (circles as preferred in forest inventory and squares as preferred in ecological surveys), and 17 plot sizes (100, 200, …, 1700 m2), hereby referred to as plot designs. To ensure that our DA-CA comparison would not be influenced by Monte-Carlo errors (Rao 1973; Koehler et al. 2009) we opted for simulating all possible samples through the construction of a sampling surface (Van Deusen et al.1986; Roesch et al. 1993; Hradetzky 1995). A sampling surface is specific to a plot design, and provides - for each point within a defined sampling frame (area) - the predicted AGB value that would be realized (from either a DA or CA HBD area level model) had the point been selected for sampling. Our sampling frame was the grid with the total of 50×106square units of 4 cm2. Upon completion we have, for each plot design, two sampling surfaces of AGB, one for the DA and one for the CA.

In our simulations we eliminated boundary effects by excluding sample locations with a distance less than or equal to 23.26 m from the study area border. The width of this sample exclusion zone corresponds to the radius of the largest circular plot. In actual applications,boundary effects will be different for the DA and the CA, but boundary effects were beyond the scope of our study.

The standard deviation of all AGB values across all 4 cm2units in a sampling surface for a specific plot design and approach (DA or CA), corresponds to the standard error under equal probability sampling and a sample size of n=1. Standard errors for larger sample sizes are not needed, as they are covered by the central limit theorem(Johnson 2004) and may be obtained by first principles.

To illustrate our results also in more practical terms,we defined a target relative standard error of 5%, and compared the expected sample sizes required to reach this target with the DA and the CA under the various plot designs. For the evaluation of sampling surfaces and determining the sample size needed to reach the precision target, we have assumed equal probability sampling.All statistical analyses and computations were made in R(R Core Team 2018).

Results

We first built the HBDs for the study area for each approach, the CA and the DA. The enlarged 20 m×20 m sections in Fig. 3 highlight the difference between the two, where the biomass distribution for the DA is characterized by high peaks at the tree positions while under the CA the tree biomass is distributed isotropically over each circular crown area as per our example in Fig. 2.

From the predicted HBD of the study stand we calculated the AGB sampling surfaces both for the DA and the CA (Fig. 4, Mg·ha?1). A sampling surface is specific for a defined plot design. It gives per point (i.e. cell of 2 cm×2 cm)the per-plot biomass prediction that results for a sample plot centred at this particular point. While the topography of a sampling surface map varies by the two approaches and by plot design, the overall (parametric)mean of AGB computed for the sampling frame of 4 cm2units remains stable for all sampling surfaces at about 279 Mg·ha?1; minor deviations from this value are caused by rasterizing instead of using a truly continuous surface,and by sampling from slightly different sampling frames due to different buffers to avoid boundary overlaps. We acknowledge that the population established by the sampling surface is built from predicted biomass values and not from the true but unknown biomass. As expected, for a given plot size, the sampling surfaces (Fig. 4) are smoother for the HBD of the CA than for the HBD of the DA. Equally,these differences diminish with an increase in plot size due to a lower within-surface variance, and due to a smaller proportion of boundary trees. With plot sizes above 900 m2the difference in the smoothness of the DA and the CA sampling surface were minor and without practical importance(not shown).

A smoother (less variable) sampling surface translates to a lower standard error in a mean estimated from a probability sample, at least when we set aside issues linked to the uncertainty in predictions of biomass and look only at the sampling induced uncertainty (as if plotlevel AGB had been observed without measurement nor model errors). Figure 5 depicts the standard deviations of the DA and the CA HBD sampling surfaces by plot size.The CA sampling surface consistently exhibits a smaller standard deviation for all plot sizes and the two plot shapes (circular, square). The relative reduction in variability is more pronounced for smaller plots and decreases as plot size increased. This observation can be explained by a higher perimeter-to-area ratio of smaller plots, so that, per-unit area, the frequency of crowns intersecting a smaller plot is higher than for larger plots. Figure 5 also indicates a small but consistently superior performance of the circular plot shape; this feature is in line with the minimum perimeter-to-area ratio of a circle.

We also used the DA and CA sampling surfaces of plot-level AGB estimates to gauge the required sample size for a target relative standard error of 5% in an estimate of the mean AGB in the studied sample frame. The results are shown in Fig. 6. As alluded to, the CA can achieve the target precision with a smaller sample size.In our study area, a reduction of 27.1% was achieved with our smallest plots of 100 m2. For the largest plots considered (1700 m2), the reduction was 10.4%. We noted a slightly but consistently greater (1.4%) reduction with circular plots than with square plots, in agreement with their more favourable perimeter-to-area ratio.

Discussion

In the context of field sampling, our case study indicates that, when assuming equal model errors, the CA approach to the estimation of AGB tends to be more efficient in terms of precision of estimation than the traditional DA. This advantage of the CA over the DA will depend, however, in a rather complex way, on tree species, tree and crown sizes, stand structures, stem densities, and the interactions of these factors with plot size and shape (in particular, perimeter length). In our case study, for a defined target precision of estimation, we observed a reduction in the required sample size of 10.4% for larger plots (1700 m2) and 27.1% for small plots (100 m2). It will be instructive to investigate all factors influencing precision and producing an error budget. Part of that research would need to be the explicit consideration of model errors when comparing the CA to the DA. In the DA, the only model error comes from the AGB model, while in the CA, there are more models in use: (1) predicting AGB from an allometric model (like with the DA), (2) predicting the circular crown projection area as a function of dbh and (3)predicting the HBD to estimate which share of the AGB of a sample tree is assigned to the particular plot.

We compared our findings regarding the precision of estimation with those of a dead wood sampling study by Gove and Van Deusen (2011), who compared the performance of their “stand-up method” (corresponding to our DA) to their “chainsaw method” (corresponding to our CA) for estimating the volume of downed coarse woody debris. They found the CA to be superior. In their study, the respective standard deviations for the sample surfaces were 135.48 m3and 120.18 m3, meaning that their (continuous) chainsaw approach was approximately 12% more precise than the (discrete) stand-up approach. Mascaro et al. (2011) also observed superiority from their version of the CA when they derived a remote sensing-based model for forest biomass prediction,even under the assumption of a very basic model of a uniform distribution for the individual tree HBD.

For this first-of-its-kind case study on the comparison of the CA and the DA for forest biomass, we worked with various simplifying assumptions regarding regular and symmetric (rotational) stem and crown shapes,branch allometry, profile distributions of foliage, isotropic distribution of biomass within the volumetric confinements of a biomass component, and perfectly upright stem axes. Combined, these assumptions may have influenced our plot-level uncertainty in biomass predictions; this will be subject of further research.

Refinements regarding our assumptions when building the tree-level HBD model are needed, and will be pursued. In particular, the use of more realistic models for each biomass component will be addressed to investigate their effect on the performance differential between the CA and the DA.

From a practical perspective, the DA for estimation of AGB has well-established advantages in terms of expediency because measurements and model-based predictions are straightforward and apply only to the sample trees within a plot. With the DA, there is neither a need for a map of tree locations in the plot, nor measurements (dbh and location) of trees in the surrounding vicinity of a plot as would be required by the CA. Additional measurements come at an often detrimental cost to existing budgets for already costly field work. However, one may realistically expect the future to bring forest mensuration devices that allow one to automatically and simultaneously measure the distance to sample trees, the height of 1.3 m, and the diameter at this reference height. This advancement would accelerate field work on mapped forest inventory plots and their surroundings and improve the appeal of the CA approach.Also, laser scanning may possibly allow us to directly model the horizontal volume distribution over a defined plot area, which could then serve as a proxy for HBD.

We compared the DA and the CA from the perspective of field plot sampling. The theoretical consistency of the CA in terms of estimating AGB over a specifically defined piece of forest (i.e. plot area) is, however, of particular relevance in model-assisted and model-based estimation problems in which remotely-sensed auxiliary variables are linked to per-plot predictions of AGB(Nelson et al. 2000; Mascaro et al. 2011; Magdon et al.2018). Conceptually, the remotely-sensed auxiliary variable only registers the biomass proxies within the confinements of a field plot while the field data from the DA - that are used to build the remote sensing based biomass model - include biomass of plot trees extending beyond the plot perimeter, and ignore the biomass from outside trees overhanging a plot. The extent of the misalignment depends (to a large degree) on the spatial resolution or density of the remotely-sensed data, plot size and in particular plot perimeter. This observation has also been described, for example, by Packalen et al.(2015) with respect to the ABA in lidar remote sensing.They proposed to extend the plot area to avoid truncations of crowns associated with plot trees.

Only from remote sensing data with a spatial resolution much finer than the size of most crowns, should we, a priori, expect the CA to generate a stronger correlation between the auxiliaries in an assisting model and AGB than the DA. This is because the misalignment of biomass, vis-à-vis the confinements of a plot, will act as a measurement error - with an attenuation of correlations (Fuller 1987; Carroll et al. 2006). In fact, Mascaro et al. (2011) in the context of an ALS study, addressed three major sources of errors: (1) GPS positional error,(2) temporal differences between field and lidar data,and (3) a lack of consistency between lidar and field plot measurements rooted in the uncertainty of deciding on whether a tree or any part of it resides inside the confinements of a field plot. With respect to the third source of error, which relates to our study, they found that the RMSE of the model that links the lidar with the ground data could be improved when a CA replaced the traditional DA. This ‘third’ error has so far not been explicitly recognized in otherwise detailed error budgets for AGB estimates in forest monitoring (Chave et al.2004; St?hl et al. 2014; Molto et al. 2013; Magnussen et al. 2014; Chen et al. 2016).

Conclusions

Our overall goal was to propose a possible improvement of the precision of forest biomass estimation from field sampling by suggesting a novel continuous approach to predict plot-level AGB. A key characteristic of the CA is the use of a HBD model applied to individual trees; thus,the CA may be applied to any forest inventory sampling design. From our study we learned that the CA produced a smoother (i.e. less variable) sampling surface.We may therefore conclude that the CA has a potential to improve the precision of field-estimated AGB - at least when we assume (as done here) that the model errors in per-plot AGB prediction are the same for the CA and the DA. However, the CA approach will not be fully operational until issues related to field measurement efforts are resolved. Moreover, the integration of the CA in forest inventories will not only depend on the availability of measurement devices for rapid mapping of trees around a sample point, recordings of dbh and delineation of crown projection areas, but also on an expansion of our arsenal of models for crown size, foliage distribution, and branch architecture.

Supplementary information

Supplementary informationaccompanies this paper at https://doi.org/10.1186/s40663-020-00268-7.

Additional file 1.Technical steps of the construction of a tree-levelHBD.

Acknowledgements

The authors are grateful to the Guest Editor, the two anonymous reviewers and Dr. Ron McRoberts for their detailed and thoughtful reviews which helped us to considerably improve the manuscript.The final language review by Ms. Mary Mulligan is highly appreciated.

Authors’contributions

CK developed the idea, drafted the research design and wrote the first draft,CPC: developed the models and programed them. JGAG, PM: supported in programming and writing; SM,LF, NN: supported implementation of the study, and data. All authors critically participated in internal review rounds and read the final manuscript and approved it.

Funding

This study was entirely self-funded,no external funding.

Availability of data and materials

Data and material will be made available upon reasonable request.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

There are no competing interests.

Author details

1Forest Inventory and Remote Sensing, Georg-August-Universit?t G?ttingen,37077 G?ttingen, Germany.2Canadian Forest Service, Natural Resources Canada, Victoria, BC V8Z 1M5,Canada.3Department of Agroforestry Engineering, University of Santiago de Compostela, E-27002 Lugo, Spain.

Received: 15 May 2020 Accepted: 30 September 2020

主站蜘蛛池模板: 国产综合欧美| 国产麻豆精品久久一二三| 欧美一道本| AV天堂资源福利在线观看| 老司国产精品视频| 在线精品欧美日韩| 国产精品主播| 不卡网亚洲无码| 亚洲欧洲国产成人综合不卡| 91精品国产情侣高潮露脸| 九九视频免费在线观看| 日本免费一级视频| 91成人免费观看| 国产aⅴ无码专区亚洲av综合网| 少妇极品熟妇人妻专区视频| 国产喷水视频| 亚洲精品波多野结衣| 超清无码熟妇人妻AV在线绿巨人 | 国产喷水视频| 伊人色天堂| 美女无遮挡免费视频网站| 一级毛片a女人刺激视频免费| 99热这里只有成人精品国产| 亚洲色图欧美在线| 久久香蕉国产线看观看精品蕉| 亚洲第一视频区| 亚洲资源在线视频| 国产精品妖精视频| 亚洲免费黄色网| 欧美亚洲激情| 亚洲区视频在线观看| 九九视频免费看| 在线观看的黄网| 97在线免费视频| 在线综合亚洲欧美网站| 国产精品午夜福利麻豆| 亚洲av成人无码网站在线观看| 日本午夜网站| 国产在线观看一区二区三区| 久久精品丝袜高跟鞋| 免费亚洲成人| 999精品视频在线| 亚洲熟女偷拍| 欧美色香蕉| 日韩小视频网站hq| 色香蕉影院| 农村乱人伦一区二区| 最新国语自产精品视频在| 亚洲综合第一区| 国产精品观看视频免费完整版| 高清乱码精品福利在线视频| 激情综合网址| 国产99视频精品免费视频7| 亚洲精品午夜天堂网页| 婷婷午夜影院| 成人第一页| 91激情视频| 99在线小视频| 亚洲综合第一页| 亚洲国产一成久久精品国产成人综合| 成人免费网站在线观看| 亚洲精品中文字幕午夜| 一级片免费网站| 在线观看网站国产| 精品国产美女福到在线不卡f| 日韩 欧美 小说 综合网 另类| 就去吻亚洲精品国产欧美| 综合久久久久久久综合网| 亚洲欧洲免费视频| 久久精品视频亚洲| 亚洲三级网站| 在线观看91香蕉国产免费| 99久久精品无码专区免费| 国产麻豆福利av在线播放| 综合天天色| 91丝袜在线观看| 亚洲精品人成网线在线 | 99精品在线视频观看| 日韩成人免费网站| 伊人无码视屏| 国产91精选在线观看| 日韩欧美综合在线制服|