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

    Local and generalized height-diameter models with random parameters for mixed,uneven-aged forests in Northwestern Durango,Mexico

    2014-02-24 07:33:52SacramentoCorralRivasJuanGabriellvarezGonzlezFelipeCrecenteCampoandJosJavierCorralRivas
    Forest Ecosystems 2014年1期

    Sacramento Corral-Rivas,Juan Gabriel álvarez-González,Felipe Crecente-Campoand José Javier Corral-Rivas

    Local and generalized height-diameter models with random parameters for mixed,uneven-aged forests in Northwestern Durango,Mexico

    Sacramento Corral-Rivas1*,Juan Gabriel álvarez-González2,Felipe Crecente-Campo2and José Javier Corral-Rivas3

    Background:We used mixed models with random components to develop height-diameter(h-d)functions for mixed,uneven-aged stands in northwestern Durango(Mexico),considering the breast height diameter(d)and stand variables as predictors.

    Conifer and broadleaves forests;h-d relationship;Mixed models;Calibration

    Background

    Most forests in Durango State(Mexico)are comprised of a mixture of species of the genera Pinus and Quercus with an irregular distribution of trees of all size classes.However, species of the genera Arbutus and Juniperus are also found in most of these forests(Wehenkel et al.2011).These forests,which cover an area of 5.4 million ha,are considered as the primary forest reserve at a national level,and they provide almost a quarter of the national forest production in Mexico(SRNyMA 2006).The forests also play an important role in providing environmental services,such as protection against soil erosion,biodiversity conservation, carbon capture and protection of water reserves;they also provide recreational areas and represent an important source of income for their owners and local inhabitants.

    Forest management requires prediction tools that provide detailed information about the development of mixed, uneven-aged stands.Growth and production models are the most commonly used tools for this purpose.When the breast height diameter(d)and total height(h)are known,application of these models is relatively easy (Sharma and Parton 2007).Measuring diameter is simple,accurate and inexpensive,whereas measuring height is relatively more complex,time-consuming and expensive.Therefore,height-diameter(h-d)functions are often utilized,so that the height of an individual tree can be predicted only from the diameter.These relationships are also very useful for estimating individual volume,site index and for describing growth and production in forest stands over time when the height is not measured(Curtis 1967).

    Most h-d functions have been developed for forest plantations(e.g.Soares and Tomé 2002;López Sanchez et al.2003).However,the relationship between the diameter and height of a tree varies between stands(Calama andMontero 2004)because it depends on stand characteristics such as density and site index(Sharma and Zhang 2004). Moreover,the h-d relationship also varies over time within the same stand(Curtis 1967).Such considerations indicate that stand variables should be used to construct generalized functions that represent all possible conditions in forest stands(Temesgen and Gadow 2004).This is particularly important in mixed,uneven-aged stands in which different species,ages,structures and levels of competition coexist (Vargas-Larreta et al.2009).

    The hierarchical structure of h-d data(i.e.trees grouped in plots and plots grouped in stands)results in a lack of independence between measurements because the observations in each sampling unit will be correlated(Gregoire 1987).Mixed models have been successfully used to address this type of problem(e.g.Lappi 1997;Calama and Montero 2004;Castedo Dorado et al.2006).This approach simultaneously estimates fixed parameters(parameters that are common to the entire population)and random parameters(parameters that are specific to each plot)within the same model and enables the variability between plots of the same population to be modelled.

    The objectives of this study were as follows:i)to compare different local h-d equations for the mixed,uneven-aged forests in north-western Durango;ii)to develop new generalized h-d equations for different groups of species based on the best local model previously fitted iii)to use the local and the generalized equations to study the capacity of mixed models to explain the variability in the h-d relationship;and iv)to determine the most suitable size and type of sample for calibrating the functions fitted with mixed models.

    Methods

    Study area

    The study was carried out in the Ejido San Diego de Tezains,Municipality of Santiago Papasquiaro,Durango State,Mexico(between 105°53′36″and 106°12′40″W and 24°48′16″and 25°13′32″N).The predominant vegetation in the area is mixed,uneven-aged forests of Pinus and Quercus.The altitude above sea level of the study area varies between 1,400 and 3,000 m.The prevailing climate is temperate:the annual precipitation ranges between 800 and 1,100 mm and the mean annual temperature varies between 8°C in the highest elevations and 24°C in the lowest elevations(García 1981).

    Data

    The data were obtained from 44 permanent plots used to monitor the growth and production of the forests in the Ejido San Diego de Tezains.These plots,which were established in 2008,were selected with the aim of representing all types of vegetation,site qualities and diameter distributions in managed stands.The plots, of size 50x50 m,are distributed under a systematic grid sampling approach that varies between 3 and 5 km, and will be remeasured at 5 year intervals.We recorded the following main variables:number of trees,species code,breast height diameter at 1.3 m(d,cm),total tree height(h,m),azimuth(°)and radius(m)from the centre of the plot(point where the diagonals cross)towards all trees of breast height diameter≥5 cm.

    The database included 25 species,which were classified on the basis of their growth patterns into the following 13 groups for posterior analysis:1(Pinus arizonica),2 (P.ayacahuite),3(P.durangensis),4(P.herrerae),5 (P.lumholtzii),6(P.teocote),7(P.douglasina),8(Quercus sideroxyla),9(other species of Quercus:Q.arizonica,Q. mcvaughii,Q.durifolia,Q.crassifolia,Q.jonesii,Q.rugosa and Q.laeta),10Pinus species(all species of the genus Pinus[codes 1 to 7]),11Quercus species(all species of the genus Quercus[codes 8 and 9]),12other conifer species(Juniperus deppeanna,J.durangensis and Cupresus lusitanica)and 13other broadleaf species(Arbutus arizonica,A.bicolor,A.madrensis,A.tesselata,A.xalapensis and Alnus firmifolia).

    We examined the distribution of the pairs of h-d data for each species or group graphically to identify any possible anomalies.As extreme data points were observed,a systematic approach,similar to the one proposed by Bi (2000)for detecting abnormal data points,was applied to increase the efficiency of the process.A local quadratic equation with a smoothing parameter of 0.25(selected after iterative fitting and visual examination of the smoothed curves for different smoothing parameters overlaid on the data),was fitted for each of the species or group.In this approach,the number of extreme values accounted for about 1%for all species together,which were excluded from the database used for fitting the equations.The main descriptive statistics for the breast height diameter and the total height of the main groups that included the species under study are shown in Table 1.

    The following stand variables were calculated from the trees registered in each plot:number of trees per hectare (N,trees ha-1),stand basal area(G,m2ha-1),mean square diameter(dg,cm),dominant height(estimated as the mean height of the 100 largest diameter trees per hectare,independently of the species[H0,m]),dominant diameter(estimated as the mean diameter of the 100 largest diameter trees per hectare,independently of the species[D0,cm])and Hart’s index(%)estimated as follows:

    Comparison of equations

    We selected a total of 27 local equations(Huang et al. 2000)for data fitting.We also studied the relationship between the stand variables and the parameters of thelocal equations that best described the h-d relationship, with the aim of improving the accuracy of the equation and developing new generalized functions.

    Table 1 Summary statistics of the database used in fitting the h-d equations

    For preliminary selection,we used ordinary non-linear least squares(ONLS)to fit each of the local equations to the data from the 13 established groups,with the MODEL procedure in SAS/ETS?statistical software package(SAS Institute Inc 2008).We evaluated the goodness of fit of the models by graphical analysis and by considering the following statistics,calculated from the residuals:root mean square error(RMSE),the coefficient of determination (R2),bias,and Bayesian information criterion(BIC;Schwarz 1978).We used the following formulae to calculate these statistics:

    As each local equation has different strengths and weaknesses,which may lead to different goodness-of-fit results for each group of species,we used a Qualification Index (QIt)to evaluate the goodness of fit by considering the values of R2(with high values representing good fits),Bias (with low absolute values representing good fits)and RMSE and BIC(with low values representing good fits). For this index,a value of 1 is assigned to the equation that was best for each group of species and a value of 0 to the others.The qualifications obtained for each equation and statistics were then summed as follows:QIijis the qualification for the j-th goodness of fit criterion in the i-th group of species. For the local equation for which the QItotalwas highest for the defined groups,we used graphical analysis and the CORR procedure in SAS(SAS Institute Inc 2008)to analyse the relation between each of the parameters and the main stand variables,with the aim of testing different forms of generalized equations.

    Effect of mixed models

    The h-d observations made in plots and stands may be highly correlated,thus violating the principle of independence of error terms(Calama and Montero 2004).One procedure used to deal with correlated observations is to fit mixed models,in which the variability between the sampling units can be explained by including random parameters,which are estimated at the same time as the fixed parameters(Lappi 1997;Calama and Montero 2004). Basically,the parameter vector of a non-linear mixed model can be defined as follows(Pinheiro and Bates 1998):

    where Φjis the parameter vector r×1(where r is the total number of parameters in the model)specified for the j-th plot,λ is the vector p×1 of the common fixed parameters for the whole population(p is the number of fixed parameters in the model),bjis the vector q×1 of the random parameters associated with the j-th plot(q is the number of random parameters in the model),Ajand Bjare matrices of size r×p and r×q for specific and random effects for the j-th plot,respectively. The basic theory of non-linear mixed models says that the residual vector(^eij)and the random effects vector (bj)are often assumed to be uncorrelated and normally distributed with mean zero and variance-covariance matrices Rjand D,respectively.The residual vector represents within subject(e.g.,plot)variability and the random effects vector represents between subject variability(Littell et al.1996).

    We constructed the non-linear mixed effects model by selecting the local and generalized equations that yieldedthe best fits for the species groups defined using the NLMIXED procedure in SAS/ETS(SAS Institute Inc 2008).We tested different combinations of fixed and random parameters and compared the fitting statistics (RMSE,R2,Bias and BIC),to determine which parameter (s)should be considered mixed.

    Calibration

    The inclusion of random parameters in h-d equations leads to two possible situations as regards prediction of the height of trees within a stand(Vonesh and Chinchilli 1997):i)a population mean response(PMR)when only diameters are measured(and the stand variables are included in the model in the case of generalized models) and the vector of random parameters is assumed to have an expected value of E(bj)=0;and ii)a calibrated response, when the height of a subsample of mjtrees is measured along with diameter measurement in each new plot j(and the stand variables in the case of generalized models)and is subsequently used to calculate the specific random parameters of the new sampling units(Calibrated Response; CR),i.e.vector bj,expressed as follows(Vonesh and Chinchilli 1997):

    Two sampling options were considered for selecting the subsample of trees to measure within each for calibration of the local and generalized equations:

    (i)CR1:Measuring the total height of between 1 and 5 randomly selected trees within each plot that are close(±10%)to the mean breast height diameter.

    (ii)CR2:Measuring the total height of the tree of mean breast height diameter,or measuring the height of two trees–the mean and minimum breast height diameters,or measuring the height of three trees–the mean,minimum and maximum breast height diameters within each plot.

    We evaluated these two alternatives in terms of the previously defined goodness-of-fit statistics(RMSE,R2and Bias),which we compared with the statistics obtained for the equations fitted by the ONLS and NLMIXED procedures.

    Results and discussion

    Local equations

    In the comparison of the goodness-of-fit statistics for the local h-d equations fitted to the data for the 13 previously defined groups of species,the Bertalanffy-Richards equation(Bertalanffy 1949;Richards 1959)consistently yielded the highest R2and lowest RMSE values;however,the equations proposed by Stage(1975)and Meyer(1940) yielded the lowest values of BIC,which gives preference to models containing few parameters over those containing several parameters(Table 2).

    The Bertalanffy-Richards equation yielded the highest R2values for 6 of the 13 species groupings and the lowest RMSE values for 4 of the groups.Finally,comparison of the BIC values indicated that this was the preferred equation only for the Quercus spp.grouping.

    In selecting the best local equation,we also examined graphs of the residuals,the significance of the parameters and the mean bias for each equation.In this respect,the Bertalanffy-Richards equation was the preferred model. The final structure of the local Bertalanffy-Richards equation used was as follows:

    where b0-b2are equation parameters and the rest of variables as defined in the data section.

    The value of the goodness-of-fit statistics for Eq.(7) and the species groups are shown in Figure 1.

    Considering that the fitting statistics for the different broad groupings(Pinus species,Quercus species, other conifer species and other broadleaf species)are similar to those obtained for each individual species and that some parameters were not significant in the individual fits for some of the species,we decided touse 4 different local equations,one for each broad grouping studied.

    Table 2 Qualification index for the 8 best local equations for the 13 groups of species

    Figure 1 Values of R2(solid line),RMSE(dashed line),and mean bias(bottom)obtained by fitting Eq.(7)to the h-d relationship for the 13 groups of species considered.

    The parameters in Eq.(7)and the goodness-of fit statistics,obtained by the ONLS method for the 4 broad groupings,are shown in Table 3.

    Table 3 Estimated parameters and fitting statistics obtained for the local model with and without mixed effects for the groups of species considered

    Generalized equations

    On relating the parameters in Eq.(7)to the stand variables, we found that in the 4 groups that included all species, parameter b0,representative of the asymptote of Eq.(7), was positively correlated(almost 52%)with H0and D0, whereas parameter b1,representative of“scale”was only positively correlated by more than 50%with Hart’s index (HI)in the other conifer species and other broadleaf species and,finally,parameter b2,representative of“shape”was also positively correlated(~51%)with dgand N for the four groups that include all species.

    To develop a generalized equation from Eq.(7),we tested various combinations of the stand variables to improve their efficacy in the fit,taking into account the previously mentioned correlations;this resulted in Eqs(8)and(9). Equation 8 yielded the best fitting statistics for the Pinus and Quercus groups,whereas Eq.(9)yielded the best fitting statistics for the species included in other conifer species and other broadleaf species.

    where b0-b2are equation parameters and the rest of variables as defined in the data section.

    On comparing the goodness-of-fit statistics for Eqs(8) and(9)with those obtained when fitting the 30 generalized equations used in other studies(e.g.,López Sanchez et al. 2003;Sharma and Zhang 2004;Sharma and Parton 2007), we found that the value of RMSE for Eq.(8)was slightly lower than those obtained by fitting the above-mentioned generalized equations,whereas the value of the RMSE for Eq.(9)was only lower than some of these.However,the latter equations included parameters that were not significant at the 0.05 level.Another advantage of Eq.(9)is that the value of the BIC was lower than that obtained for the 30 generalized equations compared.Analysis of the residuals also revealed that there were no anomalies associated with Eqs(8)and(9)that would indicate non-compliance with the underlying hypotheses of normality,homogeneity of variance or independence of errors.We thereforedecided to use Eqs(8)and(9)as generalized models for the four groups that included all species considered.The values of the parameters of Eqs(8)and(9)obtained for fitting each group of species are shown in Table 4.The signs and values of all parameters are consistent with their biological interpretation and visual examination of the graphs of the h-d relationship indicates that its performance was consistent with the theory of growth.

    Table 4 Estimated parameters and fitting statistics obtained for the generalized model with and without mixed effects for the groups of species considered

    The generalized h–d equations selected in this study included dominant stand height.This represents an advantage over equations that include the mean height because less effort is required in conventional inventories to estimate the dominant height than the mean height of the stand(López Sanchez et al.2003).These functions also include the density of the stand in terms of number of trees per unit of area and mean square diameter.Stand density is the most obvious factor affecting the h-d relationship in a stand(Zeide and VanderSchaaf 2002);in other words,trees of the same diameter are generally taller in denser stands.

    Various stand variables have been proposed as predictors of the h-d relationship:stand age(Curtis 1967;Soares and Tomé 2002;López Sanchez et al.2003);crown competition index(Temesgen et al.2007);geographic variables (Schmidt et al.2011);and wind speed(Meng et al.2008). Although the inclusion of other variables may improve the predictive capacity of the selected functions,this requires great sampling effort and limits the practical application of the functions and therefore we did not take such variables into account.

    Effect of mixed models

    Parameters b0,b1and b2in Eqs(7)and(9)determine the asymptote,the scale and the shape of the h-d curves,respectively,whereas in Eq.(8),the parameters b0and b1define the asymptote,b2is the scale parameter and b3and b4define the shape of the curve.In fitting Eq.(7)to the data from each plot,we found that the parameter that affected the asymptote was the most variable,followed by the scale parameter.Therefore,in a first step,we fitted Eqs(7),(8) and(9)to the h-d data by considering the parameters that define the asymptote and scale as mixed,in other words, with a random parameter added.Pinheiro and Bates(1998) obtained similar results and found that the best results were obtained when the asymptote and scale of Eq.(7)were considered as random parameters.In most cases,the mixed model did not converge,so that we tested the inclusion of only mixed parameters associated with those parameters that define the asymptote of the h-d curve until reaching convergence.Similar results have been reported by Sharma and Parton(2007)and Vargas-Larreta et al.(2009). Some of the parameters were scaled so that all were of the same order of magnitude and to prevent instability in the fitting function(Calama and Montero 2004).The expressions of the mixed models finally obtained are Eqs(10)to(15):

    where b0–b4are the fixed parameters of the model (common to all plots);(uj,vj)~N(0,τ)are the random parameters(specific to each plot);and^hijand eijare respectively the height and error estimated by the model for the i-th observation(tree)in the j-th plot.

    The values of the parameters and goodness-of-fit statistics for the local mixed models(Eqs 10 and 11)and for the generalized mixed models(Eqs 12 to 15)are shown in Tables 3 and 4,respectively.We compared the RMSE values obtained with the mixed effects equations with those obtained with fixed effects equations(fitted by ONLS);the values obtained with the local mixed model (Eq.10)and the generalized mixed model(Eq.12)for the Pinus grouping were 25.0%and 5.2%lower than those obtained with the local model(Eq.7)and the generalized model(Eq.8)without random parameters,respectively. For the group of Quercus species,the RMSE values obtained with the local mixed model(Eq.11)and the generalized mixed model(Eq.13)were 26.0%and 9.0% lower than those obtained with the local(Eq.7)and the generalized models without random parameters(Eq.8), respectively.For the other conifers,the RMSE values obtained with the local mixed model(Eq.11)and the generalized mixed model(Eq.14),were 9.8%and 14.3% lower than those obtained with Eqs(7)and(9),respectively. For the group comprising other broadleaf species,the RMSE values were 20.9%and 25.0%lower with the local mixed model(Eq.11)and the generalized mixed model (Eq.15)than with Eqs(7)and(9).The results obtained for BIC and R2were similar to those obtained for RMSE.

    On inspecting the graphs of the residuals for the heights estimated by the models for each species grouping,we did not find any anomalies that would suggest non compliance of underlying hypothesis of independence of errors or homogeneity of variance.The magnitude of the bias in the residual values estimated by the two fitting methods (ONLS and NLMIXED)was consistent for all ranges and classes of heights observed by the defined species groupings.

    Calibrated response

    Calibration option CR1 for Eq.(12)in Pinus species and Eq.(13)in Quercus species was the most accurate when the total height of a subsample of 3 trees close(±10%) to the mean breast height diameter for the plot was measured(Figure 2),as indicated by the slight decrease in the RMSE by 0.4%and 3.0%respectively for the two groups relative to the values estimated by the generalized model(Eq.8)fitted without random parameters.For the other broadleaf species,calibration option CR2 for Eq.(14)was the most accurate(in terms of RMSE)when a subsample of 3 trees of mean,minimum and maximum breast height diameter were measured in each plot,as the RMSE was 21%lower than that estimated with the generalized model fitted without random parameters(Eq.9). Finally,in the group of other conifer species,the RMSE value obtained with Eq.(15)and calibration option CR2 was 13%lower than that obtained with the generalized model fitted without random parameters(Eq.9).Vargas-Larreta et al.(2009)found that the decrease in RMSE varied between 3.7 and 13.3%on calibrating the model of Sharma and Parton(2007)with data from 1 tree selected at random in the plot.However,Calama and Montero(2004) observed that use of a subsample of the 5 trees with the largest diameters significantly decreased the RMSE value; Castedo Dorado et al.(2006)observed that the bias estimated with Schnute’s generalized equation was lower when a subsample of the 3 trees with the smallest diameters was used for calibration than when the estimate was made without random parameters.

    Figure 2 Value of RMSE for the calibrated equations(CR1 and CR2)fitted by NLMIXED and ONLS for the four groups of species shown.

    The variation in the value of RMSE with respect to the number of trees used with the two calibration options for the four main groups of species studied is shown in Figure 2.This statistic was also compared with those values obtained when fitting the equations by the NLMIXED (minimum value of RMSE reached only using all trees as a calibration subsample)and ONLS(maximum value of RMSE using only fixed parameters)methods.

    In the calibration process,the reduction in the RMSE value was particularly evident with the generalized mixed models for the other broadleaf species and other conifer species(21.0%and 13.0%respectively)compared with the generalized model fitted without random parameters; however,for the Pinus and Quercus groupings,the decrease in the value of this statistic was lower.Both calibration options resulted in an important reduction of RMSE for the local mixed model compared to the same model fitted without random parameters for all the groups analized.In accordance withTrincado et al.(2007),the use of a local mixed model in forest inventories with a subsample of trees to calibrate and then predict the total height of all trees not used in calibration allows retention of a simple model structure(i.e.without the need to include stand predictor variables)and may be an useful alternative to generalized mixed models when there is a lack of data to calculate stand variables.

    Conclusions

    Two generalized equations(Eqs 8 and 9)were derived from a local equation(Eq.7)and used to estimate total tree height from breast height diameter and stand variables for the 25 species identified in the sample by using mixed models.The variability between plots is explained in terms of the random effect of each plot and from the stand variables included in the generalized models.

    For species in the Pinus and Quercus groups,inclusion of the height measurements of 3 trees close(±10%)of the mean breast height diameter from each plot improved the predictive capacity of the calibrated model.For the species included in other broadleaf species and other conifer species,the predictive capacity of the model was improved by including the total height measured in a subsample of 3 trees of minimum,mean and maximum breast height diameter.The possibility of using complementary data from the stands to calibrate the mixed models provides a clear advantage over models developed by other procedures,which require large amounts of data or are less accurate.

    Competing interests

    The authors declare that they have no competing interests.

    Authors'contributions

    SCR provied the experimental data,wrote the article and analyzed the data with JGAG and FCC.JJCR supervised the work and coordinated the research project.All authors read and approved the final manuscript.

    Acknowledgements

    The present investigation was financially supported by the“Programa de Mejoramiento del Profesorado”(project:Seguimiento y Evaluación de Sitios Permanentes de Investigación Forestal y el Impacto Socioeconómico del Manejo Forestal en Norte de México).The study was conducted during the doctoral studies of the first author at the Universidad de Santiago de Compostela USC(Spain),supported by“Programa Banco Santander–USC”(becas para estancias predoctorales destinadas a docentes e investigadores de América Latina).

    Author details

    1Instituto Tecnológico de El Salto,Mesa del Tecnológico s/n.Apdo.Postal No.2,El Salto P.N.34950,Durango,México.2Departamento de Ingeniería Agroforestal,Universidad de Santiago de Compostela.Escuela Politécnica Superior,Campus Universitario,27002 Lugo,Espa?a.3Instituto de Silvicultura e Industria de la Madera,Universidad Juárez del Estado de Durango,Blvd.del Guadiana#501.,Cuidad Universitaria,34160 Durango,México.

    Received:15 May 2013 Accepted:19 September 2013

    Published:26 February 2014

    Bates DM,Watts DG(1980)Relative curvature measures of nonlinearity.J R Stat Soc 42:1–16

    Bertalanffy LV(1949)Problems of organic growth.Nature 163:156–158

    Bi H(2000)Trigonometric variable-form taper equations for Australian eucalyptus. For Sci 46:397–409

    Calama R,Montero G(2004)Interregional nonlinear height–diameter model with random coefficients for stone pine in Spain.Can J For Res 34:150–163

    Castedo Dorado F,Diéguez-Aranda U,Barrio Anta M,Sánchez Rodríguez M, Gadow KV(2006)A generalized height–diameter model including random components for radiata pine plantations in northwestern Spain.For Ecol Manage 229:202–213

    Curtis RO(1967)Height–diameter and height–diameter–age equations for second-growth Douglas-fir.For Sci 13:365–375

    García ME(1981)Modificaciones al Sistema de Clasificación Climática de K?ppen, 4ath edition.Instituto de Geografía,Universidad Nacional Autónoma de México,México D.F

    Gregoire TG(1987)Generalized error structure for forestry yield models.For Sci 33:423–444

    Hossfeld JW(1822)Mathematik für Forstm?nner,?konomen und Cameralisten, 4th edition Gotha,Hennings p 472.

    Huang S,Price D,Titus SJ(2000)Development of ecoregion-based height–diameter models for white spruce in boreal forests.For Ecol Manage 129:125–141

    Lappi J(1997)A longitudinal analysis of height/diameter curves.For Sci 43:555–570

    Littell RC,Milliken GA,Stroup WW,Wolfinger RD(1996)SAS System for Mixed Models.SAS Institute Inc.,Cary

    López Sanchez CA,Varela JG,Dorado FC,Alboreca AR,Soalleiro RR,Alvarez Gonzalez JG,Rodriguez FS(2003)A height–diameter model for Pinus radiata D.Don in Galicia(Northwest Spain).Ann For Sci 60:237–245

    Meng SX,Huang S,Lieffers VJ,Nunifu T,Yang Y(2008)Wind speed and crown class influence the height-diameter relationship of lodgepole pine:Nonlinear mixed effects modeling.For Ecol Manage 256:570–577

    Meyer HA(1940)A mathematical expression for height curves.J For 38:415–420

    Pinheiro JC,Bates DM(1998)Model Building for Nonlinear Mixed Effects Model. Department of Biostatistics,University of Wisconsin,Madison,Wis..11

    Richards FJ(1959)A flexible growth function for empirical use.J Exp Biol 10:290–300

    SAS Institute Inc(2008)SAS/ETS?9.2 User’s Guide.SAS Institute Inc,Cary,NC, p 2861

    Schmidt M,Kiviste A,Gadow KV(2011)A spatially explicit height-diameter model for Scots pine in Estonia.Eur J For Res 130:303–315

    Schwarz G(1978)Estimating the dimension of a model.Ann Stat 5(2):461–464

    Sharma M,Parton J(2007)Height–diameter equations for boreal tree species in Ontario using a mixed-effects modeling approach.For Ecol Manage 249:187–198

    Sharma M,Zhang SY(2004)Height–diameter models using stand characteristics for Pinus banksiana and Picea mariana.Scand J For Res 19:442–451

    Soares P,Tomé M(2002)Height–diameter equation for first rotation eucalypt plantations in Portugal.For Ecol Manage 166:99–109

    SRNyMA(2006)Programa Estratégico Forestal 2030.Secretaría de Recursos Naturales y Medio Ambiente del Estado de Durango,Victoria de Durango, Dgo,p 242

    Stage AR(1975)Prediction of height increment for models of forest growth. USDA For Serv Res Pap INT-164:20

    Temesgen H,Gadow KV(2004)Generalized height-diameter models–an application for major tree species in complex stands of interior British Columbia.Eur J For Res 123:45–51

    Temesgen H,Hann DW,Monleon VJ(2007)Regional height-diameter equations for major tree species of southwest Oregon.West J Appl For 22:213–219

    Trincado G,VanderSchaaf CL,Burkhart HE(2007)Regional mixed-effects height–diameter models for loblolly pine(Pinus taeda L.)plantations.Eur J For Res 126:253–262

    Vargas-Larreta B,Castedo-Dorado F,álvarez-González JG,Barrio-Anta M,Cruz-Cobos F(2009)A generalized height-diameter model with random coefficients for uneven-aged stands in El Salto,Durango(Mexico).Forestry 84(2):445–462

    Vonesh EF,Chinchilli VM(1997)Linear and Nonlinear Models for the Analysis of Repeated Measurements.Marcel Dekker Inc.,New York.560 p

    Wehenkel C,Corral-Rivas JJ,Hernández-Díaz JC,Gadow KV(2011)Estimating balanced structure areas in multi-species forests on the Sierra Madre Occidental,Mexico. Ann For Sci 68:385–394

    Weibull W(1951)A statistical distribution function of wide applicability.J Appl Mech 18:293–297

    Wykoff WR,Crookston NL,Stage AR(1982)User’s Guide to the stand prognosis model.USDA For Serv Gen Tech Rep INT-133:122

    Zeide B,Vanderschaaf C(2002)The Effect of Density on the Height-Diameter Relationship.In:Outcalt KW(ed)Proceedings of the 11th Biennial Southern Silvicultural Research Conference.2001 March 20–22.USDA Forest Service, Gen.Tech.Rep.SRS–48,Asheville,NC,Knoxville,TN,pp 463–466

    Cite this article as:Corral-Rivas et al.:Local and generalized height-diameter models with random parameters for mixed,uneven-aged forests in

    Northwestern Durango,Mexico.Forest Ecosystems 2014 1:6.

    10.1186/2197-5620-1-6

    *Correspondence:juangabriel.alvarez@usc.es

    1Instituto Tecnológico de El Salto,Mesa del Tecnológico s/n.Apdo.Postal No.2,El Salto P.N.34950,Durango,México

    Full list of author information is available at the end of the article

    ?2014 Corral-Rivas et al.;licensee Springer.This is an Open Access article distributed under the terms of the Creative

    Commons Attribution License(http://creativecommons.org/licenses/by/2.0),which permits unrestricted use,distribution,and reproduction in any medium,provided the original work is properly cited.

    Methods:The data were obtained from 44 permanent plots used to monitor stand growth under forest management in the study area.

    Results:The generalized Bertalanffy-Richards model performed better than the other generalized models in predicting the total height of the species under study.For the genera Pinus and Quercus,the models were successfully calibrated by measuring the height of a subsample of three randomly selected trees close to the mean d,whereas for species of the genera Cupressus,Arbutus and Alnus,three trees were also selected,but they are specifically the maximum, minimum and mean d trees.

    Conclusions:The presented equations represent a new tool for the evaluation and management of natural forest in the region.

    看十八女毛片水多多多| 99热网站在线观看| 国产99久久九九免费精品| 久久人人爽人人片av| 一个人免费看片子| 51午夜福利影视在线观看| 18禁黄网站禁片午夜丰满| 国产亚洲午夜精品一区二区久久| 午夜91福利影院| 午夜视频精品福利| 国产成人一区二区在线| 国产黄色免费在线视频| 国产精品 欧美亚洲| 欧美性长视频在线观看| 亚洲成人免费电影在线观看 | 啦啦啦在线免费观看视频4| 18禁裸乳无遮挡动漫免费视频| 国产激情久久老熟女| 国产1区2区3区精品| 美国免费a级毛片| 国产精品免费视频内射| 97在线人人人人妻| 国产91精品成人一区二区三区 | 国产免费一区二区三区四区乱码| 亚洲午夜精品一区,二区,三区| www.自偷自拍.com| 亚洲av美国av| 最近中文字幕2019免费版| 视频在线观看一区二区三区| 日韩大片免费观看网站| 国产精品一区二区在线不卡| 制服人妻中文乱码| 日日摸夜夜添夜夜爱| 亚洲国产精品成人久久小说| 婷婷成人精品国产| 一区二区日韩欧美中文字幕| 90打野战视频偷拍视频| 一级,二级,三级黄色视频| 操美女的视频在线观看| 久久青草综合色| 亚洲欧美一区二区三区国产| 亚洲精品一区蜜桃| 一区二区av电影网| 香蕉丝袜av| 国产精品一国产av| 国产极品粉嫩免费观看在线| 中文字幕最新亚洲高清| 亚洲欧美清纯卡通| 国产成人欧美| 久久久欧美国产精品| 肉色欧美久久久久久久蜜桃| 欧美人与性动交α欧美精品济南到| 狂野欧美激情性xxxx| 中文字幕色久视频| 日日爽夜夜爽网站| 又紧又爽又黄一区二区| www.自偷自拍.com| 亚洲久久久国产精品| 欧美xxⅹ黑人| 国产国语露脸激情在线看| 搡老乐熟女国产| 日日爽夜夜爽网站| 国产视频一区二区在线看| 伊人久久大香线蕉亚洲五| 巨乳人妻的诱惑在线观看| 亚洲第一av免费看| 在线天堂中文资源库| 国产av一区二区精品久久| 亚洲国产精品999| 一级黄色大片毛片| av线在线观看网站| 亚洲国产精品999| 人成视频在线观看免费观看| 精品福利永久在线观看| 夜夜骑夜夜射夜夜干| 2021少妇久久久久久久久久久| 中文乱码字字幕精品一区二区三区| 一边摸一边做爽爽视频免费| 久久精品国产a三级三级三级| 国产91精品成人一区二区三区 | 熟女少妇亚洲综合色aaa.| e午夜精品久久久久久久| 国产成人91sexporn| 男女边吃奶边做爰视频| xxxhd国产人妻xxx| 啦啦啦视频在线资源免费观看| 欧美性长视频在线观看| 色婷婷久久久亚洲欧美| 国产亚洲欧美精品永久| 免费久久久久久久精品成人欧美视频| svipshipincom国产片| 美女扒开内裤让男人捅视频| 成人影院久久| 国产免费视频播放在线视频| 久久国产亚洲av麻豆专区| 亚洲精品一卡2卡三卡4卡5卡 | 黄色a级毛片大全视频| 欧美乱码精品一区二区三区| 国产精品99久久99久久久不卡| 伊人久久大香线蕉亚洲五| 国产精品久久久av美女十八| 久久久国产精品麻豆| 两性夫妻黄色片| 看免费av毛片| 国产在线观看jvid| av网站免费在线观看视频| 国产高清videossex| 亚洲三区欧美一区| 成人亚洲精品一区在线观看| av片东京热男人的天堂| 夜夜骑夜夜射夜夜干| 精品视频人人做人人爽| 色精品久久人妻99蜜桃| 真人做人爱边吃奶动态| 欧美97在线视频| 老鸭窝网址在线观看| 又粗又硬又长又爽又黄的视频| 曰老女人黄片| 尾随美女入室| 七月丁香在线播放| 丝袜脚勾引网站| 亚洲国产毛片av蜜桃av| 免费少妇av软件| 久久久久久久国产电影| 国产亚洲av高清不卡| 91字幕亚洲| 日韩电影二区| 1024香蕉在线观看| 90打野战视频偷拍视频| 青草久久国产| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲情色 制服丝袜| 亚洲情色 制服丝袜| 国产成人a∨麻豆精品| 99国产精品一区二区三区| 少妇粗大呻吟视频| 国产亚洲av片在线观看秒播厂| 欧美变态另类bdsm刘玥| 国产伦人伦偷精品视频| 日本色播在线视频| 国产精品亚洲av一区麻豆| 在线观看免费高清a一片| 丰满饥渴人妻一区二区三| 日本91视频免费播放| 精品人妻在线不人妻| 亚洲,欧美,日韩| 欧美日韩亚洲综合一区二区三区_| 欧美精品人与动牲交sv欧美| 亚洲国产av新网站| 各种免费的搞黄视频| 男女边吃奶边做爰视频| 99久久综合免费| 在线观看www视频免费| 少妇精品久久久久久久| 黄色视频不卡| 欧美在线一区亚洲| 99热国产这里只有精品6| 熟女少妇亚洲综合色aaa.| 亚洲精品av麻豆狂野| 亚洲熟女毛片儿| 国产野战对白在线观看| 亚洲,欧美,日韩| 69精品国产乱码久久久| 精品一区二区三区av网在线观看 | 两个人看的免费小视频| 丁香六月天网| 18在线观看网站| 黄色毛片三级朝国网站| 少妇被粗大的猛进出69影院| 一级毛片 在线播放| 777久久人妻少妇嫩草av网站| 91精品三级在线观看| 亚洲欧美成人综合另类久久久| 老司机亚洲免费影院| 国产91精品成人一区二区三区 | 亚洲美女黄色视频免费看| 亚洲av片天天在线观看| 1024香蕉在线观看| 可以免费在线观看a视频的电影网站| 水蜜桃什么品种好| av有码第一页| 亚洲国产精品一区三区| 巨乳人妻的诱惑在线观看| 成人手机av| 一区二区三区乱码不卡18| 午夜福利一区二区在线看| 69精品国产乱码久久久| 97精品久久久久久久久久精品| 国产精品.久久久| 一级片'在线观看视频| 欧美日韩黄片免| 久久中文字幕一级| 亚洲 欧美一区二区三区| 波多野结衣av一区二区av| 中国美女看黄片| 久久 成人 亚洲| 精品国产一区二区三区四区第35| 久久精品久久久久久噜噜老黄| 老司机影院成人| 中文字幕人妻丝袜一区二区| 免费久久久久久久精品成人欧美视频| 亚洲专区中文字幕在线| 最新的欧美精品一区二区| 国产免费一区二区三区四区乱码| 亚洲成av片中文字幕在线观看| 最新在线观看一区二区三区 | 大话2 男鬼变身卡| 丝袜美足系列| 一本一本久久a久久精品综合妖精| 亚洲伊人久久精品综合| 亚洲熟女毛片儿| 夫妻性生交免费视频一级片| 色视频在线一区二区三区| 丁香六月天网| 国产精品成人在线| 国产精品久久久久久精品电影小说| 亚洲久久久国产精品| 国产精品人妻久久久影院| 欧美日韩黄片免| 欧美xxⅹ黑人| 久久精品国产综合久久久| 精品人妻在线不人妻| 一区二区av电影网| 美女视频免费永久观看网站| 亚洲精品日韩在线中文字幕| 日韩免费高清中文字幕av| 黄片小视频在线播放| 十分钟在线观看高清视频www| 亚洲欧美激情在线| 亚洲国产精品一区二区三区在线| 久久精品国产综合久久久| 黄片小视频在线播放| 日韩大片免费观看网站| 亚洲 国产 在线| 国产精品 国内视频| 精品少妇久久久久久888优播| 国产日韩欧美在线精品| 极品人妻少妇av视频| 国产日韩欧美视频二区| 国产精品av久久久久免费| 黄色a级毛片大全视频| 欧美97在线视频| 女警被强在线播放| 久久久精品区二区三区| 1024视频免费在线观看| 老司机亚洲免费影院| av有码第一页| 一级毛片女人18水好多 | 校园人妻丝袜中文字幕| 亚洲自偷自拍图片 自拍| 搡老乐熟女国产| 亚洲精品一区蜜桃| 麻豆国产av国片精品| 老司机在亚洲福利影院| 免费观看av网站的网址| 精品少妇内射三级| 99re6热这里在线精品视频| 大片电影免费在线观看免费| 久久久国产欧美日韩av| 国产成人av教育| 国产日韩欧美视频二区| 久久久亚洲精品成人影院| 午夜日韩欧美国产| 高清视频免费观看一区二区| av一本久久久久| av在线app专区| 午夜av观看不卡| a级毛片黄视频| 在线天堂中文资源库| www.熟女人妻精品国产| 精品视频人人做人人爽| xxx大片免费视频| 成年av动漫网址| 大香蕉久久成人网| 少妇 在线观看| 免费观看人在逋| 人人妻人人添人人爽欧美一区卜| 男女午夜视频在线观看| 久久久久久久久久久久大奶| 大码成人一级视频| 美女国产高潮福利片在线看| 成年av动漫网址| 99国产综合亚洲精品| av福利片在线| 丝袜美足系列| 男女无遮挡免费网站观看| 精品免费久久久久久久清纯 | 国产99久久九九免费精品| 国产精品免费视频内射| 国产精品麻豆人妻色哟哟久久| 欧美日韩综合久久久久久| 成人国产av品久久久| 亚洲成色77777| 一区二区av电影网| 久久精品国产综合久久久| 午夜精品国产一区二区电影| 日本91视频免费播放| 99国产综合亚洲精品| 亚洲色图 男人天堂 中文字幕| 久久综合国产亚洲精品| 久久精品aⅴ一区二区三区四区| 成人国语在线视频| 新久久久久国产一级毛片| 国产熟女午夜一区二区三区| 国产无遮挡羞羞视频在线观看| 日韩人妻精品一区2区三区| 欧美 亚洲 国产 日韩一| 女人久久www免费人成看片| 黑人欧美特级aaaaaa片| 人妻 亚洲 视频| 国产成人欧美在线观看 | 午夜久久久在线观看| 少妇精品久久久久久久| 狂野欧美激情性xxxx| 久久精品国产a三级三级三级| 美女主播在线视频| 午夜视频精品福利| 最黄视频免费看| 啦啦啦啦在线视频资源| 成年人黄色毛片网站| 亚洲欧美日韩另类电影网站| 日韩一本色道免费dvd| 女性被躁到高潮视频| 黑丝袜美女国产一区| 亚洲 欧美一区二区三区| 亚洲av美国av| 国产成人av激情在线播放| 热99久久久久精品小说推荐| 欧美精品av麻豆av| 国产欧美日韩一区二区三 | 亚洲精品久久午夜乱码| 十八禁网站网址无遮挡| 亚洲色图综合在线观看| 成年人免费黄色播放视频| 中国国产av一级| 80岁老熟妇乱子伦牲交| 国产片内射在线| 午夜免费观看性视频| 大片电影免费在线观看免费| 日韩精品免费视频一区二区三区| 高清av免费在线| 天天影视国产精品| 中文精品一卡2卡3卡4更新| 波多野结衣av一区二区av| 两人在一起打扑克的视频| 老熟女久久久| 精品一区二区三区四区五区乱码 | 一区二区三区精品91| 又紧又爽又黄一区二区| 欧美日韩成人在线一区二区| 国产亚洲一区二区精品| 成年人黄色毛片网站| 日韩中文字幕欧美一区二区 | 久久人人97超碰香蕉20202| 别揉我奶头~嗯~啊~动态视频 | 久久久久精品人妻al黑| 日日夜夜操网爽| 99国产精品免费福利视频| 两个人看的免费小视频| 精品人妻熟女毛片av久久网站| 中国国产av一级| 欧美老熟妇乱子伦牲交| 人人妻人人添人人爽欧美一区卜| 啦啦啦 在线观看视频| 我的亚洲天堂| 国产在线免费精品| 亚洲成人手机| 天天添夜夜摸| 99国产精品一区二区三区| 肉色欧美久久久久久久蜜桃| 国产片特级美女逼逼视频| 亚洲国产精品一区二区三区在线| 丝袜脚勾引网站| 国产精品久久久久久人妻精品电影 | 亚洲成av片中文字幕在线观看| 亚洲av电影在线观看一区二区三区| 宅男免费午夜| 天堂8中文在线网| 人妻一区二区av| 久久国产精品人妻蜜桃| 日本黄色日本黄色录像| 91精品伊人久久大香线蕉| 亚洲av日韩精品久久久久久密 | 国产成人一区二区在线| 亚洲精品久久午夜乱码| 色94色欧美一区二区| 91精品伊人久久大香线蕉| 国产精品一区二区在线观看99| 精品少妇内射三级| 中文字幕人妻丝袜制服| 夜夜骑夜夜射夜夜干| 色网站视频免费| 亚洲av片天天在线观看| 香蕉丝袜av| av国产久精品久网站免费入址| 大香蕉久久网| 真人做人爱边吃奶动态| 大陆偷拍与自拍| 99久久综合免费| 性高湖久久久久久久久免费观看| 99国产精品99久久久久| 五月开心婷婷网| 国产99久久九九免费精品| 99re6热这里在线精品视频| 久久热在线av| 一本一本久久a久久精品综合妖精| 久久ye,这里只有精品| 久久久精品免费免费高清| www.精华液| 精品视频人人做人人爽| 国产亚洲午夜精品一区二区久久| 赤兔流量卡办理| 黑人猛操日本美女一级片| 久久精品aⅴ一区二区三区四区| 午夜激情久久久久久久| 国产午夜精品一二区理论片| 久久精品久久久久久噜噜老黄| 999久久久国产精品视频| 亚洲国产av影院在线观看| www.自偷自拍.com| 欧美黄色淫秽网站| 在线观看一区二区三区激情| 777米奇影视久久| 国产在线观看jvid| 一级片'在线观看视频| 日韩 欧美 亚洲 中文字幕| 成人三级做爰电影| 精品人妻在线不人妻| 国产精品 欧美亚洲| 午夜福利一区二区在线看| 国产男人的电影天堂91| 午夜福利乱码中文字幕| 欧美另类一区| 黄片小视频在线播放| 亚洲情色 制服丝袜| 国产精品二区激情视频| 啦啦啦中文免费视频观看日本| 亚洲中文av在线| 国语对白做爰xxxⅹ性视频网站| 满18在线观看网站| 亚洲精品成人av观看孕妇| 如日韩欧美国产精品一区二区三区| 在线亚洲精品国产二区图片欧美| 在线av久久热| 欧美变态另类bdsm刘玥| 国产av国产精品国产| 精品久久久久久久毛片微露脸 | 亚洲精品在线美女| 久久这里只有精品19| 亚洲欧美色中文字幕在线| 日本猛色少妇xxxxx猛交久久| 久久亚洲国产成人精品v| 国精品久久久久久国模美| 黑丝袜美女国产一区| 美女福利国产在线| 在线精品无人区一区二区三| 国精品久久久久久国模美| 国产av精品麻豆| 黄色毛片三级朝国网站| 男的添女的下面高潮视频| 人人妻人人添人人爽欧美一区卜| 久久精品aⅴ一区二区三区四区| 少妇的丰满在线观看| 91字幕亚洲| 观看av在线不卡| 精品国产超薄肉色丝袜足j| 麻豆av在线久日| 9热在线视频观看99| 狠狠婷婷综合久久久久久88av| 王馨瑶露胸无遮挡在线观看| 亚洲精品美女久久av网站| 亚洲情色 制服丝袜| 久久久欧美国产精品| 国产亚洲午夜精品一区二区久久| 国产免费又黄又爽又色| 成人国语在线视频| 国产亚洲精品第一综合不卡| 亚洲国产av影院在线观看| 亚洲 国产 在线| 女人精品久久久久毛片| 热99久久久久精品小说推荐| 一区二区三区四区激情视频| 你懂的网址亚洲精品在线观看| 少妇 在线观看| 精品熟女少妇八av免费久了| 国产成人系列免费观看| 成人18禁高潮啪啪吃奶动态图| 日韩熟女老妇一区二区性免费视频| 老司机深夜福利视频在线观看 | 久久午夜综合久久蜜桃| 国产成人精品在线电影| 水蜜桃什么品种好| 亚洲一卡2卡3卡4卡5卡精品中文| 亚洲欧美色中文字幕在线| 国产精品一区二区在线不卡| 免费高清在线观看日韩| 久热爱精品视频在线9| 男女边吃奶边做爰视频| 大香蕉久久网| 亚洲一区中文字幕在线| 久9热在线精品视频| 午夜福利视频在线观看免费| 久久ye,这里只有精品| 久久久国产精品麻豆| 一边摸一边抽搐一进一出视频| 成年人黄色毛片网站| 91成人精品电影| 亚洲欧洲精品一区二区精品久久久| 欧美激情极品国产一区二区三区| 在线 av 中文字幕| 国产熟女午夜一区二区三区| 亚洲第一av免费看| 我的亚洲天堂| 成年美女黄网站色视频大全免费| 啦啦啦啦在线视频资源| 一本色道久久久久久精品综合| 91麻豆精品激情在线观看国产 | 午夜免费成人在线视频| 亚洲精品国产色婷婷电影| 国产在线免费精品| 下体分泌物呈黄色| 99热国产这里只有精品6| 亚洲国产看品久久| 国产99久久九九免费精品| 亚洲欧美一区二区三区黑人| 超碰97精品在线观看| 欧美黄色淫秽网站| 亚洲av美国av| 婷婷色综合www| 嫁个100分男人电影在线观看 | 丰满人妻熟妇乱又伦精品不卡| 女人久久www免费人成看片| 9热在线视频观看99| 操出白浆在线播放| 精品一品国产午夜福利视频| 久久 成人 亚洲| 亚洲黑人精品在线| 亚洲国产成人一精品久久久| 国产激情久久老熟女| 国产精品.久久久| 水蜜桃什么品种好| 国产麻豆69| 久久精品aⅴ一区二区三区四区| 国产免费视频播放在线视频| 日韩人妻精品一区2区三区| 成人国产一区最新在线观看 | 女性生殖器流出的白浆| 亚洲欧美精品综合一区二区三区| 国产精品秋霞免费鲁丝片| 午夜免费观看性视频| 热re99久久国产66热| 久久精品成人免费网站| 亚洲伊人色综图| 欧美大码av| 欧美激情高清一区二区三区| 国产精品久久久久久精品古装| xxx大片免费视频| 欧美日本中文国产一区发布| 亚洲国产精品成人久久小说| 又大又黄又爽视频免费| 久热这里只有精品99| 在线观看免费高清a一片| 男女边摸边吃奶| 国产日韩欧美亚洲二区| 两个人看的免费小视频| 日本a在线网址| 久久久久久人人人人人| 黑人猛操日本美女一级片| 久久国产精品大桥未久av| 黄网站色视频无遮挡免费观看| 国产成人欧美在线观看 | 超碰成人久久| 最黄视频免费看| 久久精品aⅴ一区二区三区四区| 亚洲成人手机| 大香蕉久久成人网| 巨乳人妻的诱惑在线观看| 中文字幕av电影在线播放| 天天躁夜夜躁狠狠久久av| 黑人巨大精品欧美一区二区蜜桃| 中国国产av一级| 成年美女黄网站色视频大全免费| 国产高清视频在线播放一区 | 久久久久国产精品人妻一区二区| 日韩一卡2卡3卡4卡2021年| 久久av网站| 欧美国产精品va在线观看不卡| 精品国产乱码久久久久久小说| 两个人免费观看高清视频| 久久精品国产综合久久久| 久久九九热精品免费| 男人添女人高潮全过程视频| 午夜福利视频精品| 亚洲人成电影观看| 日本欧美视频一区| 精品高清国产在线一区| 久久精品成人免费网站| 欧美精品高潮呻吟av久久| av国产精品久久久久影院| 亚洲激情五月婷婷啪啪| 波多野结衣av一区二区av| 亚洲色图综合在线观看| 视频在线观看一区二区三区| 国产精品麻豆人妻色哟哟久久| 啦啦啦视频在线资源免费观看| 18在线观看网站| 国产免费又黄又爽又色| 色播在线永久视频| 国产一区二区三区av在线| 久久亚洲精品不卡| 一区二区三区乱码不卡18| 别揉我奶头~嗯~啊~动态视频 | 一区二区三区精品91| 男女无遮挡免费网站观看| 亚洲欧美一区二区三区久久|