Jie Zou ,Wei Hou,Ling Chen,Qianfeng Wang,Peihong Zhong,Yong Zuo,Shezhou Luo and Peng Leng

Abstract

Keywords:Sampling scheme,Elementary sampling unit,Clumping index,Leaf area index,Digital hemispherical

Introduction

Leaf area index (LAI) is an important forest structure parameter that is defined as one half of the total canopy element green area per unit ground area of forest plots(Watson 1947; Chen and Black 1992). LAI is required in most forest biophysical and physiological process models(Running and Hunt 1993; Zou et al. 2018b). It can be measured at the elementary sampling unit (ESU) scale via direct or indirect methods, and can be retrieved from remote sensing images at the regional, national or even global scale. ESU-scale LAI estimates are typically used to evaluate satellite-derived LAI products due to their relatively high accuracies and reliabilities.

LAI is one of the 16 essential terrestrial climate variables required by the global climate observing system(GCOS). GCOS has specified that the estimation error of LAI products must be limited to 20% and should be improved to 5% for future applications (Fernandes et al.2014; Woodgate 2015). However, large differences have been reported amongst commonly used LAI products(Garrigues et al. 2008; Fang et al. 2013; Chen 2014).Therefore, accurate ESU-scale LAI estimates are required to consistently evaluate LAI products in order to pick out qualified LAI products and improve those unqualified ones.

Optical methods are typical indirect methods to obtain the ESU LAI of forests with high efficiency and low cost(Jonckheere et al. 2004; Zou et al. 2009), including LAI-2000/LAI-2200(LI-COR,Lincoln,NE,USA),digital hemispherical photography (DHP) (Jonckheere et al. 2004;Weiss et al. 2004), multiband canopy imaging methods such as multiband vegetation imager(Zou et al.2009)and multispectral canopy imager(Kucharik et al.1997),as well as the tracing radiation of canopy and architecture method (3rd Wave Engineering, Winnipeg, Manitoba,Canada) (Leblanc and Chen 2002). Previous studies reported that the LAI estimation error of optical methods in the ESU LAI estimation of forests are mainly caused by inversion models (Liu et al. 2015; Zou et al. 2018b), clumping effects (Leblanc and Chen 2002; Leblanc et al. 2005;Chen et al. 2006; Zou et al. 2018b), nonphotosynthetic components (Weiss et al. 2004; Zou et al. 2009; Ercanl? et al.2018;Zou et al.2018a),terrain slopes(Gonsamo and Pellikka 2008; Cao et al. 2015), sampling schemes(Nackaerts et al. 2000; Weiss et al. 2004; Majasalmi et al. 2012; Pfeifer et al. 2018) and observation conditions (Leblanc and Chen 2001; Jonckheere et al. 2004;Zhang et al. 2005). In recent years, great efforts have been made to reduce the ESU LAI estimation errors of optical methods, and numerous advancements have been achieved through the involvement of appropriate inversion models (Leblanc and Fournier 2014;Woodgate 2015; Zou et al. 2018b), terrain slope correction methods (Cao et al. 2015) and woody components (nonphotosynthetic components) correction methods (Zou et al. 2018a; Zou et al. 2019). However,to the best of our knowledge, no attempt has been made to evaluate the impact of sampling schemes on the ESU LAI estimation of forests via optical methods. Sampling schemes with various numbers and spatial arrangements of sample points have been adopted by optical methods such as DHP and LAI-2000/LAI-2200, to estimate the ESU effective plant area index (PAIe), effective woody area index and LAI of forest plots. For example, sample points were spatially arranged in certain patterns such as circle(Soto-Berelov et al. 2015; Woodgate 2015), square(Neumann et al. 1989; Baret et al. 2005; Macfarlane et al. 2007; Ryu et al. 2010a; Pisek et al. 2011;Woodgate et al. 2012; Liu et al. 2015), transect(Cutini et al. 1998; van Gardingen et al. 1999; Hyer and Goetz 2004; Abuelgasim et al. 2006; Ryu et al.2010a; Majasalmi et al. 2012), cross (Leblanc 2008;Woodgate et al. 2012; Leblanc and Fournier 2014;Zou et al. 2018a) and dispersed (Abuelgasim et al.2006; Majasalmi et al. 2012). Moreover, different sample sizes ranging from 1 to 176 have been applied at the ESU scale (Cutini et al. 1998; van Gardingen et al. 1999; Nackaerts et al. 2000; Macfarlane et al.2007; Gonsamo et al. 2010; Ryu et al. 2010a; Pisek et al. 2011; Woodgate 2015; Calders et al. 2018; Zou et al. 2018a; Zou et al. 2018b).

Previous studies have reported obvious differences in the PAIeestimates from different sampling schemes(Nackaerts et al. 2000; Majasalmi et al. 2012; Calders et al. 2018). Similar to the PAIeobtained indirectly from gap fraction measurements, canopy element clumping index (Ωe), which can be estimated by current available algorithms such as gap-size analysis (CC) (Chen and Cihlar 1995a), finite-length averaging (LX) (Lang and Yueqin 1986) and combination of gap-size and finitelength averaging (CLX) (Leblanc et al. 2005), also depends on gap fraction and gap size measurements from the optical method.Therefore, Ωeestimation should also be affected by sampling schemes.

Different zenith angles (θ) have been used by different inversion models. For instance, the effective zenith angle range for the Nilson model (Nilson 1999) and the Miller theorem (Miller) (Miller 1967) is from 0° to 90°, while the zenith angle range for the LAI calculation method of the LAI-2200 instrument is from 0° to 74° (LAI-2200)(LI-COR 2009), and zenith angles near 57.3° are used for Beer’s law (Beer) (Ross 1981; Leblanc and Fournier 2014;Zou et al. 2018b). Variations in gap fraction or Ωemeasurements with zenith angles ranging from 0° to 90° are usually large for forest plots (Leblanc et al. 2005; Liu et al.2015;Woodgate 2015;Zou et al.2018b).Therefore,LAI estimation is expected to be influenced by inversion models which depend on gap fraction and Ωemeasurements. For example, the plant area index (PAI) estimates of forest plots were reported to be largely affected by inversion models (Zou et al. 2018b). Therefore, accurate ESU-scale LAI should be obtained by the joint consideration of inversion models, Ωealgorithms and sampling schemes.

The following three aspects of sampling schemes need further investigation to improve LAI estimates from DHP. First, the PAIeestimated from DHP is 35%-75%of the PAI of forest canopies if Ωeis ignored in the estimation (Leblanc and Chen 2002; Zou et al. 2018b).Therefore, the impact of sampling schemes on PAI estimates should be started by investigating the impact of sampling schemes on Ωeestimation. However, no attention has been paid to such impact yet. Second, large differences were observed in the PAI or LAI of the same forest plot obtained from optical methods with different inversion models (Ryu et al. 2010a; Zou et al. 2018a; Zou et al. 2018b). However, only one inversion model (i.e.Beer or LAI-2200) was used to evaluate the impact of sampling schemes on the PAIeestimation of forest plots(Nackaerts et al. 2000; Majasalmi et al. 2012; Calders et al. 2018). Thus, whether the impact of sampling schemes on LAI estimation varies with different inversion models should be investigated. Third, recommendations should be given to select optimal sampling schemes for accurate ESU LAI estimation of forests.

In this study, DHP images were collected at 89 unique sample points of 13 sampling schemes in five typical Larix principis-rupprechtii forest plots. Three commonly used inversion models (i.e. Miller, LAI-2200 and Beer)and Ωeestimation algorithms (i.e. CC, LX and CLX)were adopted to estimate the ESU LAI of forests. The main purposes of this study are to: (1) analyse the impact of sampling schemes on gap fraction, Ωe, PAIeand LAI estimations; (2) evaluate the performance of sampling schemes on the ESU LAI estimation of forests using litter collection LAI measurements; and (3) identify the best sampling scheme for the ESU LAI estimation of forests.

Theory

The LAI of forest canopies can be obtained on the basis of the PAI from DHP and woody-to-total area ratio (α)measurements:

In the current study, three models (i.e. Miller, LAI-2200 and Beer) were used because they were reported to outperform other commonly used inversion models in the PAI estimation of forest canopies (Zou et al. 2018b).Several algorithms have been proposed to estimate the Ωeof forest canopies using gap fraction or gap size measurements collected from DHP. They include CC, LX,CLX, modified gap-size analysis (Pisek et al. 2011),modified finite-length averaging (Gonsamo and Pellikka 2009) and Pielou’s coefficient of spatial segregation(Walter et al. 2003). Pielou’s coefficient of spatial segregation algorithm was claimed to produce inaccurate Ωeestimates (Walter et al. 2003; Pisek et al. 2011; Zou et al.2018b). High similarities were found between CC and modified gap-size analysis as well as between LX and modified finite-length averaging algorithms (Zou et al.2018b). Moreover, amongst CC, LX and CLX, no single algorithm is better than the other two in the Ωeestimation of all forest plots (Zou et al. 2018b). Therefore, all three algorithms (i.e. CC, LX and CLX) were adopted in the PAI estimation in this study. The details of the formula to obtain the Ωeand PAI of forest plots can be found in Additional file 1.

Materials and methods

Plot description

The study area is located in the Saihanba National Forest Park in Hebei Province, China. Five even-aged plots of 25 m×25 m (ESU) that covered typical ages of L.principis-rupprechtii in the study area were used. L.principis-rupprechtii is one of the dominant and widespread tree species in northern China. The five plots are at least about 120 m away from the forest border. The structure characteristics of the canopy around the plots are similar to those of the plots. The majority of branches below the live canopy in plot 1 and 2 were harvested by local forest managers before our field experiment. Dead branches at heights below 2.5 m in plots 3 and 4 were harvested in the low and medium levels before the field campaign (Zou et al. 2018a). The main characteristics of the five plots are listed in Table 1.

Data acquisition and processing DHP

Thirteen sampling schemes (SS1-SS13) (predefined sampling schemes) comprising 89 unique sample points were used in this study (Fig. 1). Figures 1a-d show four square sampling schemes, amongst which 9, 16 and 25 sample points were evenly distributed in the plots for SS2, SS3 and SS4,respectively.Figures 1e-g show three cross sampling schemes (i.e. SS5, SS6 and SS7), amongst which the sample points were located at or close to the two perpendicular centre lines of each plot. Compared with SS6, one additional sampling point was added in each quadrant 2.5 m away from the centre lines for SS7 to increase the canopy sampling of the centre areas of the plots (Leblanc 2008). Two dispersed sampling schemes are shown in Fig.1h and i.For SS8,three sample points were randomly distributed within the plot along with one central point and eight others were distributed evenly on the boundary lines. For SS9, four sample points were distributed evenly along one of the diagonal lines and 12 sample points were distributed evenly on the four boundary lines. Figure 1j-l illustrate three transect sampling schemes, in which sample points were evenly distributed along the diagonal lines.A total of 13 sample points in a circle sampling scheme(i.e. SS13) were evenly distributed with a distance interval of 5 m along three intersecting lines with the same intersection angle of 60° (Woodgate 2015). The number of sample points for these 13 sampling schemes ranged from 3 (i.e. SS10) to 25 (i.e. SS4). To increase the number of sample points, a sampling scheme (SS14)(simulated sampling scheme) was generated using all the sample points from these 13 sampling schemes to represent the sample scheme with 89 sample points(Fig.1n).For SS14,the distances between the neighbour points for all sampling points ranged from 0.7 to 5.9 m (Fig. 1n). The position of each sampling point of all sampling schemes in the field was established with tape and marked with wooden stakes with flags to facilitate the DHP image collection (Zou et al. 2018a).Parts of the 89 sampling points were slightly moved from its designed positions to new adjacent positions in the field in order to make the positions of all sampling points were at least about 0.6 m away from its nearest large tree stems. This movement is necessary in order to avoid the greatly increased visible proportion of stems in DHP images.

DHP images were collected using a Canon 6D camera with a Sigma 8 mm fisheye lens before sunrise,after sunset or under overcast conditions (Zou et al.2018a). The DHP image resolution was 5472×3648 pixels. Manual exposure mode was used in the DHP image collection. The exposure was determined using the method described in Leblanc et al. (2005) and Woodgate et al. (2015) to achieve good contrast between canopy element and sky. The camera was set approximately 1.2 m above the ground. The DHP images were acquired one by one at all 89 unique sample points of the 13 sampling schemes in each plot. A total of 445 DHP images were obtained from the five plots. The DHP images were taken in 2017 from August 11 to September 2, which was the peak season with maximum LAI of the year for the five plots(Zou et al. 2018a; Zou et al. 2019).

The DHP images were manually processed using Adobe Photoshop 7.0 (Adobe Systems, San Jose, CA,USA). The image preprocessing procedures included image cropping, blue channel selection, gamma correction (the gamma of the blue channel was corrected to about 1.0) and image classification (Gonsamo andPellikka 2008; Zou et al. 2018a). The threshold values for the DHP images to be classified into canopy element and sky were determined by referring to the mean digital numbers of the canopy element and sky peaks in the histogram of the blue channel images.

Table 1 Characteristics of the Larix principis-rupprechtii plots (Zou et al. 2018a)

In this study, a narrow zenith angle range of 52°-62°was used by the Beer inversion model. Segment size is a key parameter required by LX and CLX in Ωeestimation. This parameter was used to divide each annulus of DHP into several small segments and make the canopy element at the segment scale approach random spatial distribution which was assumed by LX. The 5° segment size recommended by Gonsamo et al. (2010) was adopted for LX (hereinafter, LX_5, in which ‘5’ represents the segment size of 5°); this size would result in 72 segments for each annulus of DHP. As for CLX, no large differences were found between the root mean square error and mean absolute error of the PAI of forest plots estimated using the same inversion model and CLX but with three segment sizes of 15°, 30° and 45° (Zou et al.2018b). Therefore, for simplicity, only the 15° segment size was chosen for CLX (CLX_15) to calculate Ωein this study. For CC, two transect processing schemes (i.e.CCW and CCS) were used in the Ωeestimation. CCW merges the transects of all DHP images into a whole transect to calculate the Ωefor each sampling scheme.By contrast, CCS treats the transect of each DHP image individually and the final Ωeis obtained by averaging the Ωeof all DHP images for each sampling scheme. The same estimation procedure to obtain Ωeat zenith angles ranging from 0° to 90° with a 1° interval, as suggested by our previous study (Zou et al. 2018b), was also used in the current study. MTVSP software (version 2018) (Zou et al. 2015; Zou et al. 2018a) was operated to estimate the canopy element gap fraction at θ (pe(θ)), Ωe, PAIe,PAI and LAI of the five plots from the preprocessed DHP images. PAIeand PAI of each plot were calculated using the mean pe(θ) of all processed DHP images of each sampling scheme, as suggested by Ryu et al.(2010b). For each sampling scheme of each plot, 3 PAI estimates were calculated using Equations A6-A8 on the basis of the pe(θ), canopy element gap fraction of the ith annulus of LAI-2200(pe_i), canopy element gap fraction at 57° (pe(57)) and Ωeestimate, which was obtained from one Ωealgorithm amongst the four Ωealgorithms(i.e., CCS, CCW, LX_5 and CLX_15). Then 12 PAI estimates were obtained for each sampling scheme of each plot. Next, 12 LAI estimates were calculated for each sampling scheme of each plot using Equation (1) based on the calculated PAI and α measurements, which were obtained with destructive measurements (detailedly described in section below).

The pe(θ) measurements of annuli with zenith angles close to 90° were tended to approach or equal to zero.In this study, PAI or LAI under the condition of zero gap fraction was defined as 10 (Zou et al. 2018b). The upper limit of the estimated LAI for very small pe(θ)measurements was set to 10 to avoid obtaining unrealistically high LAI estimates (Leblanc et al. 2005; Pisek et al. 2011; Yan et al. 2019).

α, γ and LAI

The α for each plot was estimated with destructive measurements of two or three representative trees that were selected and harvested. The stems and branches of the harvested trees were divided into height classes with a height interval of 1 m and a starting height of 1.2 m. The stem,branch and fruit areas of each height class for each harvested tree were measured with a tape or a digital calliper, with the assumption that the stem or branch sections were circular truncated cones and the fruits were spheroids (Zou et al. 2018a). Approximately 300-350 typical needles were randomly selected from the branches of each height class or several height classes with similar needle characteristics. The needle cross section was assumed to be rectangular in shape. The hemisurface area of typical needles was measured with the volume displacement method described by Chen et al.(1997). The dry mass of typical needles and all needles of each height class of each tree were obtained by drying the needles at 65°C for 48 h until the weights of the needles remained almost unchanged. For these typical needles, the specific leaf area was obtained by dividing the hemisurface area by the dry mass. The leaf area of a height class was calculated by multiplying the specific leaf area with the dry mass of all needles in that height class. The leaf area of each harvested tree was the sum of the leaf area of every height class for that tree. The α of each harvested tree was then calculated by dividing the area of the woody components (stem, branches and fruits) by the sum of the woody components and leaf areas of that tree. The α of each plot was calculated using the α of the representative trees and the DBH measurements of each plot. Details of the α estimation of the five plots are available in our previous work (Zou et al. 2018a).

For each plot, two to four typical shoots were randomly selected from each height class of forest canopy(i.e. bottom, middle and top) for needle-to-shoot area ratio (γ) determination. Incomplete and abnormal (untypical) shoots were discarded before the estimation of γ. For each typical shoot, the projection areas were recorded using a Canon 6D camera equipped with a Canon 24-70 mm lens and a flat, levelled white panel with two rulers laid on its top surface (Zou et al. 2018a).Three shoot projection images were taken right above each shoot sample at a distance of approximately 0.6 m by rotating the shoot main axis to three zenith angles(i.e. 0°, 45° and 90°) and one azimuth angle (i.e. 0°).Afterwards, the three projection areas of the typical shoot (Ap(0°,0°), Ap(45°,0°) and Ap(90°,0°)) were calculated according to the three shoot projection images(Zou et al. 2018a). The half of the total needle area in a typical shoot (An) was measured by the volume displacement method described by Chen et al. (1997). Then, the effective needle-to-shoot area ratio (γe) of each shoot could be obtained with Equation (A1). The γeof each plot was obtained by averaging the γeof all typical shoots. Then the γ of each plot was obtained with Equation (A2) on the basis of γeand α.

In each plot, nine litter traps of 50 cm×50 cm in size were placed 50 cm aboveground at each sample point of SS2. The litter collection was started before defoliation and ended after it. The measurements were performed six times during the entire period, on September 1, 18 and 29 and October 13, 18 and 24 in each plot. During each measurement time, approximately 300-350 typical needles were randomly selected from the litter collections of each plot to measure the specific leaf area. The method for the calculating the specific leaf area of the typical needles was the same as those for calculating α. The dry mass of all needles of each measurement time for each plot was obtained by drying the needles at 65°C for 48 h until the weights of the needles remained almost unchanged. Totally six specific leaf area and dry mass measurements corresponding to the six measurement times, respectively, were obtained for each plot. In each plot, the one-half of the area of fallen needles for each measurement time was obtained by multiplying the specific leaf area by the dry mass of all fallen needles of that measurement time. The LAI of each plot was then obtained by dividing one-half of the total area of the fallen needles collected during the six measurement times by the area of the litter traps.

Statistical analysis

Mean relative error (MRE), which is similar to the mean relative deviation from Calders et al. (2018), was used in this study because it corresponds to the accuracy evaluation required by GCOS. The MRE of the LAI for each combination of inversion model, Ωealgorithm and sampling scheme was calculated using the litter collection LAI as the reference.

Results

Effect of sampling schemes on gap fraction (pe(θ))estimation

Obvious differences were observed amongst the pe(θ)from different predefined sampling schemes at the four typical zenith angles, except at 90°, and the differences seem to decrease with the zenith angle (Fig. 2). For example, the maximum pe(θ) differences amongst the 13 predefined sampling schemes in the five plots ranged from 0.29 to 0.44 for zenith angle 0°, 0.07-0.16 for 30°,0.03-0.09 for 57° and 0- 6.81×10?4for 90°. Note that the pe(90) of all predefined sampling schemes in the five plots, except for plots 4 and 5, were equal to 0 (Fig. 2);pe(90) of all predefined sampling schemes, except for SS8-SS9, were equal to 0 in plot 4 (Fig. 2d); and the pe(90) of all predefined sampling schemes ranged from 3.07×10?5to 7.12×10?4in plot 5 (Fig. 2e).

Effect of sampling schemes on canopy element clumping index (Ωe) estimation

Similarly to the gap fraction, obvious variations were also found amongst the Ωefrom the 13 predefined sampling schemes, which were obtained using the same Ωealgorithm at four typical zenith angles (i.e. 10°, 30°, 57°and 90°) in plot 5. The variations also seem to decrease with the zenith angle, except at 90° (Fig. 3). For example,the ranges of maximum variations amongst the Ωefrom the 13 predefined sampling schemes for the four Ωealgorithms were 17%-28% (0.07-0.16) for the zenith angle of 0°, 9%-22% (0.07-0.13) for 30° and 5%-10% (0.05-0.08) for 57° (Fig. 3). Variations were also observed at 90°, especially for algorithms LX_5 and CLX_15. The ranges of variations amongst the Ωefrom the 13 predefined sampling schemes at 90° were 0%-0% (0-0), 0%-1% (0-0.01), 0%-37% (0-0.27) and 0%-35% (0-0.26) for CCW, CCS, LX_5 and CLX_15,respectively (Fig. 3d).

Figure 3 also shows that the Ωeestimated with CCW was greater than that by CCS, which was obtained with the same scheme and zenith angle, for all predefined sampling schemes at the four zenith angles, except at 90°. In addition, the differences between the Ωefrom CCW and CCS under the same sampling scheme tended to increase with the increasing sample size at almost all zenith angles, except at 90°. For example, the differences between the Ωefrom CCW and CCS at 30° increased from 0.02 (3%) for SS2 to 0.12 (15%) for SS3 and 0.15(18%) for SS4. An example of the gap removal procedure of the Ωeestimation for CCW and CCS was presented(Fig. 4) to investigate the cause of the difference between Ωefrom CCW and CCS. For CCW, the gap fraction of the largest gap size was small (8.62×10?4for the gap size of 132 in pixels in the Fig. 4a example) because the length of the whole transect was large (153, 225 pixels)(Fig. 4a).

Effect of sampling schemes on PAIe and LAI estimation

Obvious differences were observed amongst the PAIeof the 13 predefined sampling schemes, which were estimated using the same inversion model in the five plots.The strength of the sampling scheme’s impact on the PAIeestimation varied with the inversion model (Fig. 5).The widest range of differences (0%-43%) amongst the PAIeof the 13 predefined sampling schemes was found if Miller was used in the PAIeestimation. By contrast,the range of differences amongst the PAIeof the 13 predefined sampling schemes was narrowed to 0%-32% for Beer and 0%-30% for LAI-2200.

Fig. 6 shows that the LAI estimation was also largely affected by the sampling schemes. However, it is different from the PAIeestimation since no large difference of the impact on LAI estimation was found amongst the three inversion models when the Ωealgorithm of CCW was not used in the estimation. Specifically, the differences between the LAI obtained using the same inversion model and Ωealgorithm amongst the four Ωealgorithms, except CCW, under different predefined sampling schemes ranged from 0% to 36.2% (0-0.88) for Miller, 0%-35.1% (0-0.92) for LAI-2200 and 0%-37.6%(0-0.9) for Beer. When the Ωealgorithm was CCW, the differences between the LAI obtained using the same inversion model under different predefined sampling schemes became larger, which ranged from 0% to 46.8%(0-0.79) for Miller, 0%-37.9% (0-0.69) for LAI-2200 and 0%-64.5% (0-0.75) for Beer.

Fig.6 also shows that all the LAI estimates obtained using the three inversion models and four Ωealgorithms with 14 sampling schemes in the five plots were systematically smaller than the litter collection LAI. The underestimation was especially obvious if CCW was adopted in the estimation, and it was reduced if CCW was replaced by CCS in the estimation (Fig. 6 and Table 2). For example, as shown in Table 2, the MREs of LAI obtained using Beer with 13 predefined sampling schemes ranged from 47.2% to 59.7%for CCW and reduced to 44.2%-53.0%for CCS.The smallest MRE was 15.2%,which was obtained using the combination of Miller,CLX_15 and SS10.

Table 2 Mean relative errors (MREs)(%) of the LAI obtained with all the combinations of inversion model,canopy element clumping index(Ωe) algorithm and sampling scheme. For each combination of inversion model and Ωe algorithm, the four sampling schemes with the smallest MREs shown in red and the largest shown in green

As shown in Table 2, the MREs of the LAI obtained using Miller and CLX_15 with the scheme with the largest sample size (SS14) was larger than those obtained using the same inversion model and Ωealgorithm but with all predefined sampling schemes with smaller sample sizes (SS1-SS13), except for SS6 and SS7. The same tendency can be found for the cases of two sampling scheme groups with the same sampling scheme type but with different sampling sizes, such as SS2-SS4 and SS10-SS11. For instance, the MREs of the LAI obtained using Miller and CLX_15 were 19.9% for SS2 with a sample size of 9, 21.0% for SS3 with a sample size of 16 and 24.3% for SS4 with a sample size of 25 (Table 2).Similarly, the MREs were 15.2% for SS10 with a sample size of 3 and increased to 20.0% for SS11 with a sample sizes of 5 if Miller and CLX_15 were adopted in the LAI estimation (Table 2).

Table 2 also shows that the performance of sampling schemes in the LAI estimation changes with the inversion models. For example, SS2 and SS10 outperformed the other 11 predefined sampling schemes in the LAI estimation if the inversion model was LAI-2200 or Beer due to the MREs of SS2 and SS10 were amongst the four smallest MREs of the LAI obtained using the same inversion model and Ωealgorithm but with the 13 predefined sampling schemes (Table 2). Similarly, the MREs of SS10 and SS11 were amongst the four smallest MREs of the LAI obtained using the same inversion model and Ωealgorithm with the 13 predefined sampling schemes if Miller was used in the estimation (Table 2). Note that most of the MREs of the LAI obtained with SS5, SS6,SS7 and SS9 were greater than those obtained using the same inversion model and Ωealgorithm but with the other 9 predefined sampling schemes (Table 2).

Discussion

Are the Ωe and LAI estimates markedly affected by sampling schemes?

Large differences were found amongst the Ωe(0%-37%)obtained with the same Ωealgorithm under different predefined sampling schemes (Figure 3) and amongst the LAI (0%-64.5%) estimated with the same inversion model and Ωealgorithm under different sampling schemes in all the plots (Fig. 6). This finding indicates that both the Ωeand LAI estimation are largely affected by the sampling scheme.

The large differences for Ωeunder different sampling schemes could be explained by the obvious differences for pe(θ) under different sampling schemes (Fig. 2) given that all the Ωealgorithms directly or indirectly relied on pe(θ) to calculate Ωe. The fact that the difference between Ωeof CCW and CCS increases with the sample size for all the typical zenith angles, except at 90°, is primarily attributed to the incapability of CCW to effectively remove the large gaps (Fig. 4), which resulting from the nonrandom distribution of the canopy element in the gap removal procedure of the Ωeestimation (Chen and Cihlar 1995b). This tendency indicated that the transect processing scheme is an important issue for CC in the Ωeestimation. With the increase of transect length or sample size, the gap fraction of the same large gaps would be smaller. The variations between the LAIs obtained on the basis of the gap fraction estimates of the whole transect with or without the large gap sizes would tend to be smaller than the trigger condition of the gap removal procedure of CCW (<0.01) (Chen and Cihlar 1995b). Larger Ωeestimates would thus be obtained by CCW because large gaps were not effectively removed by the gap removal procedure in CCW (Fig. 4a). On the contrary, large gaps of the transects could be effectively distinguished and removed by CCS (Fig. 4b). Therefore,the differences between the Ωeof CCW and CCS increase with the sample size. Actually, many studies have also reported that CCW offers larger Ωeestimates than other Ωealgorithms, even though different PAI and plant function types of forest plots were covered in these studies (Pisek et al. 2011; Leblanc and Fournier 2014;Woodgate 2015; Zou et al. 2018a, 2018b). CCS could offer more accurate Ωethan CCW due to the MREs of the LAI obtained using CCS were smaller than that using CCW under the condition of the same inversion model and sampling scheme (Table 2).

The obvious differences amongst the pe(θ) or Ωefrom different sampling schemes would further result in large differences in PAIeor LAI (Figs. 5 and 6, respectively)because pe(θ) and Ωeare the key parameters for PAIeor LAI estimation (Equations A6-A8). The impact of the sampling schemes on PAIeestimation was greater when the inversion model was Miller rather than the other two (i.e. LAI-2200 and Beer). An indication was that the variations of PAIeobtained with different predefined sampling schemes ranged from 0% to 43% for Miller but only 0%-30%for LAI-2200 and 0%-32%for Beer(Fig.5).One reason was identified in relation to the large range of variations for Miller. The pe(90) of all predefined sampling schemes in the five plots, except in plots 4 and 5,were equal to 0 (Fig. 2); the pe(90) of all predefined sampling schemes, except for SS8 and SS9, were equal to 0 in plot 4 (Fig. 2); and the pe(90) of all predefined sampling schemes ranged from 3.07×10?5to 7.12×10?4in plot 5 (Fig. 2). For Miller, the small variations between pe(θ) measurements with zenith angles close to 90°, especially for those between zero and non-zero measurements, would result in obvious differences in PAIeestimates. This outcome is attributed to the variations in pe(θ) measurements being significantly enlarged by the logarithm of pe(θ) in the PAIeestimation (Equation A6)and the relatively large weight values (cos(θ)sin(θ)) at zenith angles close to 90° (not equal to 90°) compared to those at zenith angles close to zenith for Miller (Zou et al. 2018b). Compared with Miller, LAI-2200 or Beer can avoid the impact of zero or very small pe(θ) measurements on the PAIeestimation owing to the small zenith angle range covered by LAI-2200 (0°-74°) or Beer(52°-62°)and the large zenith angle width of the annulus of LAI-2200 (12°-14°) or Beer (10°) (Zou et al. 2018b).

Different from PAIe, the impact of sampling schemes on LAI estimation was not greater for Miller than the other two inversion models if CCW was not used in the estimation (Fig. 6 and Table 2). The reason could be attributed to the different zenith angle ranges covered by the three inversion models and the tendency of the variations amongst the pe(θ) or Ωeof different sampling schemes decreases with the zenith angle, except at 90°(Figs. 2 and 3).

Does increasing the sample size always improve LAI estimation accuracy?

The MREs of the LAI obtained with SS14, which has the largest sample size (i.e. 89), were not always smaller than those of the LAI obtained using the same inversion model and Ωealgorithm but with all predefined sampling schemes with smaller sample sizes (i.e. from 3 to 25) (Table 2). The reason might be the overlap of canopy sampling of sample points in SS14 owing to the large sample size and relatively sparse spatial distribution of the sample points for SS14. The degree of overlap depends on the distance between the sample points and the canopy height, which is difficult to quantify because of the irregular spatial distribution of the sample points of SS14 in the heterogeneous forest plots. Furthermore,the MREs of the LAI obtained with SS4 were not the smallest amongst its sampling scheme group (i.e. SS2-SS4) to estimate LAI using the same inversion model and Ωealgorithm, even though it had the largest sample size in the group. The same phenomenon was observed for the group of SS10 and SS11. Therefore, the accuracy of the LAI estimates obtained from DHP cannot always be improved with increasing sample sizes. Caution is needed if the LAI estimates obtained with the sampling schemes having large sample sizes are treated as the estimates with high accuracy or even as reference estimates.Previous studies (e.g. Nackaerts et al. 2000; Majasalmi et al. 2012; Calders et al. 2018) treated the PAIeor effective woody area index estimates from sampling schemes with large sample sizes as highly accurate results or even as reference for validation which should provide more specific analyses.

Which sampling scheme(s) is (are) more reliable to estimate the LAI of forest plots?

The best sampling schemes for estimating the ESU L AI of L. principis-rupprechtii forests varied with the different inversion models used (Table 2). The two sampling scheme groups of SS2, SS10 and SS10, SS11 are recommended for use to estimate the ESU LAI of L. principisrupprechtii forests if LAI-2200 or Beer and Miller were adopted in the estimation. A common characteristic of SS2, SS10 and SS11 is their small sample sizes (i.e. 3-9).This means that an accurate estimation of the ESU LAI of L. principis-rupprechtii forest plots does not need sampling schemes with large sample sizes, which are very cost-effective for the long time series field measurements of forest ESU LAI to validate LAI products.

Amongst the 13 predefined sampling schemes, SS5-SS7 and SS9 are not recommended for use owing to their relatively large MREs compared with the MREs of the LAI obtained using the same inversion model and Ωealgorithm, except for CCW, but with the other nine predefined sampling schemes (Table 2).SS6 and SS7 are likewise not recommendable because they have the largest MREs in their own sampling scheme groups (with the same or similar sample sizes but with different sampling scheme types). Specifically, SS6 belongs to the group that made up SS2, SS6 and SS12 with a sample size of 9, whilst SS7 belongs to the group of SS7, SS8 and SS13 with the same or similar sample sizes ranging from 12 to 13. The large MREs for SS5-SS7 may be caused by an undersampling or oversampling of the four corners or centre areas of the plots because the sampling points of those sampling schemes are distributed along or near the two perpendicular centre lines of the plots. The poor performance of SS5-SS7 in this study is inconsistent with the conclusions of Majasalmi et al.(2012), who concluded that cross sampling schemes similar to SS5-SS7 in this study outperformed other sampling schemes in the PAIeestimation of forest plots. Three reasons can be responsible for the conflicting conclusions from the two studies. First, only the impact of sampling schemes on the PAIewas investigated in Majasalmi et al. (2012), which should be different from that on LAI given that Ωeis additionally required for LAI estimation. Furthermore, Ωewas largely affected by the sampling schemes (Fig. 3). Second, different references were used to evaluate the accuracy of retrievals from different sampling schemes in the two studies. The PAIeestimates of sampling scheme with the largest sample size were used as reference in Majasalmi et al. (2012), whereas litter collection LAI measurements were used in the current study. However, the results of this study showed that LAI obtained from the sampling scheme with the largest sample size (SS14) is not the one with the highest accuracy (Fig. 6 and Table 2). Third, the sampling schemes used in the two studies are not completely the same, especially for the square and dispersed sampling schemes.

Can accurate ESU LAI with estimation errors below 5%or 20% be obtained from DHP using the 13 predefined sampling schemes?

LAI underestimation was observed for all the combinations of inversion model, Ωealgorithm and sampling scheme in the five plots (Fig. 6). This finding can be attributed to two reasons. First, previous studies reported that optical methods tend to underestimate the LAI of forest plots when the LAI (not PAI) is large, especially for a closed canopy with an LAI larger than approximately 5 or 6, due to the gap fraction saturation issue for optical methods (Gower et al. 1999; Jonckheere et al.2004). Evident LAI underestimation of optical methods in coniferous forest plots with large LAI values was also reported by Chen et al. (1997). They presented a large difference of 45% between the LAI of optical methods and the allometric measurements in an old black spruce plot with an LAI of 6.3. As one of the optical methods,DHP underestimated LAI in plot 5 with an LAI of 6.69 obviously in this study (Fig. 6). Second, most shoots whose petioles measure about 1-4 mm are located directly on the surfaces of the branches and stems of L.principis-rupprechtii forest plots. Therefore, the overlap of shoots and woody components will make DHP unable to sample canopies sufficiently, leading to Ωeoverestimation and further LAI underestimation. Moreover, because branches in the bottom level of plots 1 and 2 were harvested by management activities, most shoots that were detected by DHP sensor were located in the middle and upper level of the canopies. Long-distance detection reduces the effectiveness of canopy sampling and therefore, results in LAI underestimation.

Amongst all the combinations of inversion model,Ωealgorithm and sampling scheme,the LAI with the two smallest MREs (i.e. 10.7% and 13.5%) were obtained using the combinations of Miller and SS10 with the two Ωealgorithms of CLX_15 and LX_5, respectively (Table 2). The two MREs of our study are close to those of Leblanc and Fournier (2014), who reported a minimum MRE of 11%for the PAI obtained from DHP if the appropriate inversion model and Ωealgorithm were adopted in the PAI estimation of the simulated forest scenes. It was expected that one of the two MREs of this study would be slightly larger than those of the simulation study. Two reasons were identified in relation to the relatively large MRE of this study. First, compared with the simulation study, the field LAI measurements of DHP usually include two more LAI estimation error sources (i.e. observation conditions and DHP image classification). Second, the estimates obtained from the simulation study are PAI, not the LAI of this study. The conversion from PAI into LAI would introduce additional estimation errors from the α measurements.

Table 2 shows that the ESU LAI estimates obtained by our study did not match the maximum LAI estimation error threshold of 5% set by GCOS. The reason for this is that the MREs of the LAI obtained using all the combinations of inversion model, Ωealgorithm and sampling scheme are larger than 5% (Table 2). The LAI differences between litter collection LAI measurements and those obtained by DHP with all the possible combinations of inversion models, Ωealgorithms and sampling schemes are also larger than 5% in the five L. principisrupprechtii plots, except in plot 4 (Fig. 6). However, ESU LAI estimates with the MREs below 20% could be obtained if appropriate combinations of inversion model,Ωealgorithm and sampling scheme were adopted in the LAI estimation of L. principis-rupprechtii forest plots(Table 2). Specifically, ESU LAI estimates with relatively small MREs (ranging from 10.7% to 20.0%) can be obtained for all plots, except plot 5, using Miller and CLX_15 with SS1-SS13, with the exception of the four sampling schemes of SS5-SS7 and SS9. The LAI differences between litter collection LAI measurements and those obtained from DHP using Miller and CLX_15 with SS10 or SS11 are lower than 20% in the five plots, except plot 5 (Fig. 6). This result indicates that ESU LAI estimates with the maximum LAI estimation errors below 20%,which is required by GCOS, can be obtained from DHP if the appropriate inversion model, Ωealgorithm and sampling scheme were adopted. This conclusion is at least effective for L. principis-rupprechtii forest plots, except for those with extremely large LAI values (e.g. 6.69 for plot 5), because wide canopy structure characteristics were covered by the five plots in this study (Table 1).

Conclusions

The impact of sampling schemes on the ESU LAI estimation from DHP in the five L. principis-rupprechtii forests was evaluated in this study. Results showed obvious differences amongst the pe(θ), Ωe, PAIeand LAI from different sampling schemes. The differences in impact of the sampling schemes on LAI estimation were not obvious amongst the three inversion models (i.e. Miller, LAI-2200 and Beer) if the Ωealgorithm of CCW was not adopted in the estimation. The accuracy of the ESU LAI estimates was not always improved by increasing the sample sizes. The two sampling scheme groups of SS2,SS10 and SS10, SS11 (with sample sizes ranging from 3 to 9) outperformed other predefined sampling schemes to obtain the ESU LAI of L. principis-rupprechtii forests if Beer or LAI-2200 and Miller were adopted in the estimation, respectively. SS5-SS7 and SS9 are not recommended for use owing to their relatively large MREs.Moreover, ESU LAI estimates with maximum LAI estimation errors below 20%, which is required by GCOS,could be achieved by DHP if the appropriate inversion model, Ωealgorithm and sampling scheme were adopted in the LAI estimation of L. principis-rupprechtii forests,except for those with extremely large LAI values. However, DHP is still not qualified for obtaining ESU LAI estimates of L. principis-rupprechtii forests with maximum LAI estimation errors below 5%, regardless of which combination of inversion model, Ωealgorithm and sampling scheme is adopted in the estimation.

Future work can include efforts to investigate the impact of sampling schemes on the ESU LAI estimation of forests from other indirect LAI measurement methods(e.g. terrestrial laser scanner and LAI-2200). Since the differences between the canopy structures of different plant functional types, caution is needed if the conclusion of this study is applied to forests with tree species other than the L. principis-rupprechtii covered in this study. Therefore, forest plots with different plant functional types can also be included to consolidate the conclusions of this study. Furthermore, the complicated issue of canopy sampling overlap between different sample points is worth analysing in the future.

Nomenclature

Supplementary information

Supplementary informationaccompanies this paper at https://doi.org/10.1186/s40663-020-00262-z.

Additional file 1. PAI estimation.

Acknowledgements

Fundingfor this work came from the National Science Foundation of China(Grant Nos.41871233, 41371330 and 41001203).

Authors’ contributions

Jie Zou participated in study design,field measurements,data processing,interpreting the results and writing the manuscript.Peng Leng took part in the field measurements.Wei Hou,Yong Zuo and Peihong Zhong took part in the procedure of data processing.Ling Chen,Qianfeng Wang and Shezhou Luo participated in the writing of the paper.The author(s)read and approved the final paper.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable requests.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

1The Academy of Digital China (Fujian), Fuzhou University, Fuzhou 350116,China.2Key Laboratory of Data Mining and Information Sharing, Ministry of Education, Fuzhou 350116, China.3School of Forestry, Beijing Forestry University, Beijing 100083, China.4College of Environment and Resources,Fuzhou University, Fuzhou 350116,China.5College of Resources and Environment,Fujian Agriculture and Forestry University, Fuzhou 350002,China.

Received: 9 April 2020 Accepted: 21 July 2020