Performance of Four Optical Methods in Estimating Leaf Area Index at Elementary Sampling Unit of Larix principis-rupprechtii Forests

Optical methods are frequently used as a routine method to obtain the elementary sampling unit (ESU) leaf area index (LAI) of forests. However, few studies have attempted to evaluate whether the ESU LAI obtained from optical methods matches the accuracy required by the LAI map product validation community. In this study, four commonly used optical methods, including digital hemispherical photography (DHP), digital cover photography (DCP), tracing radiation of canopy and architecture (TRAC) and multispectral canopy imager (MCI), were adopted to estimate the ESU (25 m × 25 m) LAI of five Larix principis-rupprechtii forests with contrasting structural characteristics. The impacts of three factors, namely, inversion model, canopy element or woody components clumping index ( Ω e or Ω w ) algorithm, and the woody components correction method, on the ESU LAI estimation of the four optical methods were analyzed. Then, the LAI derived from the four optical methods was evaluated using the LAI obtained from litter collection measurements. Results show that the performance of the four optical methods in estimating the ESU LAI of the five forests was largely affected by the three factors. The accuracy of the LAI obtained from the DHP and MCI strongly relied on the inversion model, the Ω e or Ω w algorithm, and the woody components correction method adopted in the estimation. Then the best Ω e or Ω w algorithm, inversion model and woody components correction method to be used to obtain the ESU LAI of L. principis-rupprechtii forests with the smallest root mean square error (RMSE) and mean absolute error (MAE) were identified. Amongst the three typical woody components correction methods evaluated in this study, the woody-to-total area ratio obtained from the destructive measurements is the most effective method for DHP to derive the ESU LAI with the smallest RMSE and MAE. In contrast, using the woody area index obtained from the leaf-off DHP or DCP images as the woody components correction method would result in a large LAI underestimation. TRAC and MCI outperformed DHP and DCP in the ESU LAI estimation of the five forests, with the smallest RMSE and MAE. All the optical methods, except DCP, are qualified to obtain the ESU LAI of L. principis-rupprechtii forests with an MAE of <20% that is required by the global climate observation system. None of the optical methods, except TRAC, show the potential to obtain the ESU LAI of L. principis-rupprechtii forests with an MAE of <5%.

. [23]. Ω e and corrected needle-to-shoot area ratio (γ c ) estimation formulae used in this study. γ is the needle-to-shoot area ratio. F m (0, θ) is the measured total canopy element gap fraction at θ, F mr (0, θ) is the total canopy element gap fraction after removing the large gaps resulting from the non-random distribution of the canopy element [39,40], p e (θ) is the canopy element gap fraction at θ, p e (θ) is the mean canopy element gap fraction of all segments at θ, ln[p e (θ)] is the mean logarithmic canopy element gap fraction for all segments at θ [37], n is the segment number, p e_k (θ) is the canopy element gap fraction of segment k and Ω e_CC_k (θ) is the Ω e of segment k [32]. ∅ and f f are the crown porosity and foliage cover, respectively [41]. A n is half the total needle area in a shoot. A p (0 • , 0 • ), A p (45 • , 0 • ), and A p (90 • , 0 • ) are the shoot projection areas at azimuth angle 0 • and zenith angles 0 • , 45 • , and 90 • , respectively [42].

.2. PAI, WAI, and LAI
The inversion model is a key factor that affects the PAI or LAI estimation from optical methods due to the large differences observed between the PAIs or LAIs obtained using different inversion models [19,23,24,35]. Two fundamental inversion models were proposed previously and are widely used in the ESU LAI estimation of forests, i.e., Beer's Law (Beer) (Equations (8) and (12)) [44] and Miller theorem (Miller) (Equations (9) and (13)) [45]. However, these two fundamental inversion models are unsuitable for all optical methods because the zenith angle ranges covered by these methods differ from those of the Miller theorem (0 • -90 • ), such as LAI-2000/LAI-2200 (0 • -74 • ) and DCP (0 • -15 • ) methods. Therefore, several inversion models were proposed with the development of these optical methods [14,19,41,46]. For example, Equation (10) or (14) is the modified Miller integration adopted by the LAI-2200 method (LAI-2200) in estimating the PAI or LAI of forest canopies [23,46]. Similarly, Equation (11) is also the modified Miller integration adopted by the MCI method (MCI_0-85) for estimating the PAI of forest canopies [14,23]. Table 2 lists the inversion models used in the present study. Table 2. The inversion model formulae used for the plant area index (PAI), wood area index (WAI), and leaf area index (LAI) estimations in this study. γ c is the corrected needle-to-shoot area ratio. G e (θ) is the canopy element projection coefficient (G e ) at θ. G e_i is the G e of the i th annulus. Ω e (θ) is the Ω e at θ derived from the gap size distribution (CC), logarithmic averaging (LX), and combination of gap size and logarithmic averaging (CLX) using Equations (1)-(3), respectively. p e (57) and Ω e (57) are the canopy element gap fraction and Ω e measurements obtained at the zenith angle range of 52 • -62 • , respectively. θ i , p e_i (θ i ), Ω e_i , and W i are the center zenith angle, canopy element gap fraction, Ω e and weight factor of the i th annulus of digital hemispherical photography (DHP) or multispectral canopy imager (MCI), respectively. For the LAI-2200 method, the zenith angle ranges of the five annuli are 0 • -13 • , 16 • -28 • , 32 • -43 • , 47 • -58 • , and 61 • -74 • , respectively [19,46]. The W i of the five annuli of LAI-2200 are 0.041, 0.131, 0.201, 0.29, and 0.337, respectively [19,46]. For MCI, the zenith angle ranges of the nine annuli are 0  (11)) [23]. The W i of the nine annuli of MCI_0-85 are 0.0038, 0.0303, 0.0596, 0.0872, 0.1120, 0.1335, 0.1510, 0.1640, and 0.2580, respectively [23]. f c is the crown cover. α is the woody-to-total area ratio. α i is the α of the i th annulus of MCI. α DCP is obtained by averaging the α of the first and second annuli of MCI derived using CC.

PAI =
−ln(p e (θ))cos(θ)γ c G e (θ)Ω e (θ) (7) [44] PAI Beer = −2ln(p e (57))cos(57)γ c /Ω e (57) (8) [19,44] PAI Miller = −2 π 2 0 ln[p e (θ)]γ c Ω e (θ) cos(θ)sin(θ)dθ (9) [45] [19,46] ln[P e_i (θ i )]γ c cos(θ i )W i G e_i Ω e_i (11) [23] LAI Beer = −2ln(p e (57)) cos(57)γ c (1 − α 7 )/Ω e (57) (12) [19,44] [41] Zou et al. [19] reported that the impact of the assumption of spherical angle distribution of canopy element or woody components (canopy element or woody components projection coefficient is equal to be 0.5) on the ESU PAI and WAI estimations of the leaf-on and leaf-off forests from DHP can be reduced to a low level (<4%) if the three inversion models of Miller, LAI-2200 and Beer are adopted in the estimation. Therefore, that assumption was made in this study. The equations of the three inversion models (e.g., Beer, Miller, and LAI-2200) to estimate the WAI of the leaf-off forests from DHP are the same as those used for estimating the PAI of the leaf-on forests, respectively (Equations (8)- (10)). Similarly, the equations of the two inversion models (e.g., Beer and MCI_0-85) to estimate the WAI of the leaf-off forests from MCI are the same as those used for estimating the PAI of the leaf-on forests, respectively (Equations 8 and 11). For DCP, the equations used to estimate the PAI and WAI of the leaf-on and leaf-off plots are the same to derive LAI (Equation (15)) by assuming that the woody-to-total area ratio value is equal to zero. If the PAI and WAI measurements are simultaneously available, then the LAI of leaf-on forests also can be calculated as follows: Forests 2020, 11, 30 6 of 28

Woody-to-Total Area Ratio (α)
The α of leaf-on forests can be derived based on the WAI and PAI as follows [14]: If the WAI and PAI in Equation (17) are derived using Equation (7) based on the canopy element gap fraction (p e_i (θ i )), woody components gap fraction (p w_i (θ i )), Ω e_i , and Ω w_i of the i th annulus of the MCI, respectively, then the estimate derived from Equation (17) is α i .

Plot Description
Five forest plots with varying LAIs, stand densities, mean tree heights and tree ages (Table 3) were selected as the long-term observation ESUs with a size of 25 m × 25 m for validating LAI map products. The selected plots are located in the Saihanba National Forest Park in Hebei Province, China [23]. The dominant tree species in the park is L. principis-rupprechtii, which is a tree species that is spread widely in the northern area of China [47]. The five plots are single-species plots with the tree species of L. principis-rupprechtii, and the tree ages covered by the five even-aged plots are the typical tree ages of the park. The plots are at least 120-300 m away from the forest border with a flat terrain [47]. The canopy structure characteristics (e.g., tree height, diameter at breast height (DBH), tree age, stand density, and LAI) of the canopy around the plots are similar to those of the canopies of the plots [47]. Forest inventory was collected during the field campaign in 2017. Table 3 presents the site description of the five plots. The TRAC measurements were collected in the five plots using TRAC-II (Nanjing Huiming Instrument Inc., Nanjing, China). The leaf-on TRAC measurements were collected between 11 August and 12 September 2017 at the maximum PAI (maximum LAI) of the plots [48][49][50]. The leaf-off TRAC measurements were collected between 30 October and 15 November 2017 at the minimum PAI (LAI = 0) of the plots. For the leaf-on plots, approximately six TRAC measurements were collected with the zenith angles within the maximum range of local solar elevation angle (0 • -70 • ) with an interval of 10 • . For the leaf-off plots, approximately four TRAC measurements were collected with the zenith angles within the maximum range of local solar elevation angle (0 • -40 • ) with an interval of 10 • . The TRAC instrument was operated at a height of approximately 1.2 m. The transect length of each TRAC measurement ranged from 100 m to 170 m and it was divided into several 10-m sub-transects using the flags. All TRAC measurements were collected under clear days. The TRAC measurements were processed using TRACWin software (version 5.5.4). Some TRAC measurements were not correctly stored by the instrument because of the Bluetooth module failure. Thus, these measurements were not used in this study. The measurements with the photosynthetic photon flux density equal to zero were not used in this study because no gaps were detected by TRAC. Only the TRAC measurements with zenith angles near 57 • were adopted to derive the ESU LAI of the five forests.

DCP
The DCP images were collected at a resolution of 5472 pixels × 3648 pixels using a Canon EOS 6D camera equipped with a Canon 24-70 mm lens. The focal length of the camera was fixed at 24 mm. The camera height was approximately 1.2 m. Manual mode was used to avoid overexposure and underexposure of the image collection. All images were collected under an overcast sky, before sunrise or after sunset. A sampling scheme with 16 sampling points that were evenly distributed within the plot with a 5-m distance was adopted for DCP image collection. The leaf-on DCP images were collected between 11 August and 3 September 2017 at the maximum PAI of the forests [48][49][50]. The leaf-off DCP measurements were collected between 30 October and 9 November 2017 at the minimum PAI of the forests. A total of 160 DCP images were obtained from the five forests. The center areas of the original DCP images were clipped into a new image with a resolution of 2227 pixels × 2061 pixels (a field of view of approximately 15 • × 15 • ) before further processing. Then, the clipped images were processed using Adobe PhotoShop 7.0 (Adobe Inc., USA) following the procedures described by Macfarlane et al. [41]. Subsequently, the Ω e and Ω w of the five leaf-on and leaf-off forests for DCP were calculated using Equation (4) following the same procedures described by Macfarlane et al. [41], respectively. Two LAI estimates were obtained for each plot. One LAI estimate was derived using Equation (15) from the leaf-on DCP images. Another LAI estimate was obtained based on the derived PAI and WAI measurements using Equation (16). The PAI and WAI of leaf-on and leaf-off plots were obtained using Equation (15) by assuming that the woody-to-total area ratio value is equal to zero based on the leaf-on and leaf-off DCP images, respectively.

DHP
The sampling scheme of DHP was the same as that of DCP. DHP images were collected using a Canon 6D camera equipped with a Sigma 8 mm fisheye lens [23]. The DHP image resolution was 5472 pixels × 3678 pixels. The camera height was approximately 1.2 m. Manual mode was used to avoid overexposure and underexposure of the DHP image collection [23]. The DHP images were collected before sunrise, after sunset, or under overcast sky conditions. The leaf-on DHP images were collected between 11 August and 2 September 2017 at the maximum PAI of the five forests [48][49][50], whilst the leaf-off DHP images were collected between 30 October and 15 November 2017 at the minimum PAI of the five forests. A total of 160 DHP images were obtained from the five plots.
The procedures described in Gonsamo and Pellikka [51] were adopted to process the DHP images, including selecting the blue channel, applying a correction to the gamma and thresholding to classify the canopy element and sky. The classified images were then inputted into MTVSP (version 2018) [23,52] to calculate p e (θ), p w (θ), Ω e (θ), Ω w (θ), WAI, PAI, and LAI. The two segment sizes of 5 • and 15 • were selected for LX and CLX in the Ω e or Ω w estimations, respectively. The detailed explanation of the choice of the two segment sizes can be found in the works of Gonsamo and Pellikka and Zou et al. [47].
Three Ω e (θ) or Ω w (θ) estimates were derived by using CC, LX, and CLX at each zenith angle in the zenith angle range of 10 • to 90 • with an interval of 1 • in each forest (one estimate for each Ω e or Ω w algorithm). To obtain 91 Ω e (θ) and Ω w (θ) estimates that match the p e (θ) and p w (θ) measurements at the same zenith angle range of 0 • to 90 • with an interval of 1 • in each forest, the Ω e (θ) or Ω w (θ) at the zenith angles ranging from 0 • to 9 • were assumed to be equal to Ω e (10) or Ω w (10) [19]. Two LAI groups (9 estimates of each group for the nine combinations of the three Ω e algorithms and inversion models) were calculated for each forest using Equations (12)-(14) based on p e (θ), p e_i (θ i ), Ω e (θ), Ω e_i , corrected needle-to-shoot area ratio, and destructive or MCI woody-to-total Forests 2020, 11, 30 8 of 28 area ratio measurements. In addition to the LAI estimates derived using Equations (12)- (14) based on the leaf-on DHP images, another alternative to derive the ESU LAI of the five forests was based on the PAI and WAI estimates of the leaf-on and leaf-off forests, respectively. One PAI group (9 estimates for the nine combinations of the three Ω e algorithms and inversion models) was calculated for each leaf-on forest using Equations (8)-(10) based on the p e (θ), p e_i (θ i ), Ω e (θ), Ω e_i , and corrected needle-to-shoot area ratio measurements. One group of WAI (9 estimates for the nine combinations of the three Ω w algorithms and inversion models) was calculated for each leaf-off forest using Equations (8)-(10) based on the p w (θ), p w_i (θ i ), Ω w (θ), and Ω w_i measurements (corrected needle-to-shoot area ratio was assumed to be equal to 1). Then, another LAI group (9 estimates) could be derived using Equation (16) based on the two groups of PAI and WAI. Thus, 27 ESU LAI estimates were obtained for each forest.

MCI
The detailed procedures of collecting and processing MCI images can be found in the work of Zou et al. [23]. We provide a brief description here. Both leaf-on and leaf-off MCI images were collected from the same two representative sampling points in each forest. The MCI image pairs of each sampling point were collected in zenith directions ranging from 0 • to 80 • with an interval of 10 • [23]. The MCI instrument height was approximately 1.2 m. Six azimuthal measurements were collected at each zenithal direction with an interval of 60 • [14,23]. The MCI image resolution was 3488 pixels × 2616 pixels. The leaf-on and leaf-off MCI images of the five forests were collected at the same or similar period of the leaf-on and leaf-off DHP or DCP image collection, respectively. A total of 1080 MCI image pairs were collected from the five forests. The MCI image pairs were processed using ENVI 4.7 (Harris Geospatial Solutions Inc., USA) with the "IsoData" classification method [14]. The MCI image pairs were eventually classified into three classes of sky, woody components, and shoots, respectively [23]. The full size of each pixel of the MCI images was calculated based on the field of view of MCI, the image resolution, mean tree height of the plot and zenith angles where the images were collected [14,23]. Then, the pixel number in each segment which is required as the key parameter for LX and CLX to derive Ω e or Ω w can be obtained based on the derived values of the full size of each pixel. The classified MCI images were inputted into MTVSP (version 2018) [23,52] to calculate p e (θ), P e_i (θ i ), p w (θ), P w_i (θ i ), Ω e (θ), Ω e_i , Ω w (θ), Ω w_i , WAI, PAI, and LAI of each forest.
For MCI, both MCI_0-85 and Beer were used for the LAI estimation. The seventh annulus (55 • -65 • ) of MCI was adopted by the Beer inversion model in the LAI estimation. One group of six PAI estimates (one estimate for each combination of the Ω e algorithm and the inversion model) was calculated using Equations (8) and (11) based on the p e_i (θ i ), p e (57), Ω e_i , and Ω e (57), which were calculated from the classified leaf-on MCI images, and corrected needle-to-shoot area ratio measurements. Two groups of twelve WAI estimates (one estimate for each combination of the Ω w algorithm and inversion model) were derived using Equations (8) and (11) based on the p w_i (θ i ), p w (57), Ω w_i , and Ω w (57) measurements that were calculated from the classified leaf-on and leaf-off MCI images, respectively. Then, 12 LAI estimates were derived from the MCI in each forest using Equation (16) based on the one PAI group and two WAI group measurements.
For MCI, only the MCI_0-85 inversion model was used for the woody-to-total area ratio determination of the five forests. The woody-to-total area ratio values of the five forests were calculated using Equation (17) based on the obtained PAI and WAI measurements, which were derived using Equation (11). Detailed procedures for obtaining the woody-to-total area ratio of the five forests from MCI can be found in Zou et al. [23].

Mean Element Width, Litter Collection LAI, Corrected Needle-to-Shoot Area Ratio and Woody-to-Total Area Ratio
The mean element width, litter collection LAI, corrected needle-to-shoot area ratio and destructive woody-to-total area ratio of the five plots were obtained during the field campaign. The detailed procedure of obtaining the mean element width, corrected needle-to-shoot area ratio, and destructive woody-to-total area ratio measurements of the five plots can be found in the work of Zou et al. [23]. Here, we provide a brief description. The mean element width of each plot was obtained using typical shoot samples clipped from the three height classes of forest canopies (i.e., top, middle, and bottom) (two to four typical shoot samples for each height class). Then, the method described by Leblanc et al. [53] was used to calculate the mean element width of each typical shoot. The mean element width of each plot is the average of the mean element width values of all typical shoots of each plot [23]. Another group of 2-4 typical shoots was clipped from each height class of the canopy (i.e., top, middle, and bottom) for the needle-to-shoot area ratio calculation in each plot. The projection areas of each typical shoot were obtained using a Canon 6D camera equipped with a Canon 24-70 mm lens by rotating the shoot's main axis at azimuth angle 0 • and zenith angles 0 • , 45 • , and 90 • , respectively [23]. The half the total needle area of each typical shoot was obtained using the volume displacement method described by Chen et al. [36]. Then, the needle-to-shoot area ratio of each typical shoot was calculated using Equation (5) based on the three projection areas and half the total needle area measurement of the shoot. The needle-to-shoot area ratio of each plot was obtained by averaging the needle-to-shoot area ratio of all typical shoots of each plot [23]. Thereafter, the corrected needle-to-shoot area ratio of each plot was calculated using Equation (6) based on the needle-to-shoot area ratio and woody-to-total area ratio measurements.
The destructive method was used to measure the destructive woody-to-total area ratio values of five plots. Two or three representative trees were harvested in each plot. Then, the woody components (i.e., stem, branch, and fruit) area of each harvested tree was measured using a tape measure or digital caliper by assuming that the stem or branch sections were circular truncated cones and the fruits were spheroids [23]. The leaf area of each harvested tree was calculated using the dry mass of all needles and the specific leaf area of the harvested tree. The specific leaf area of the harvested tree was determined based on the dry mass and the hemisurface area of a number of typical needles (approximately 300-350) was estimated using the volume displacement method [36]. Afterward, the woody-to-total area ratio value of each harvested tree was calculated based on the woody components and needle areas of all needles of each tree. Finally, the destructive woody-to-total area ratio value of each plot was obtained based on the woody-to-total area ratio values of the harvested trees and the DBH measurements of all trees of the plot [23].
The litter collection LAI of the five plots was obtained using the litter collection measurements. Nine litter traps with a size of 50 cm × 50 cm were evenly placed in each plot with a distance of 6.25 m. The trap height was approximately 0.5 m. Six litter collection measurements were collected from 1 September 2017 to 28 October 2017 with a time interval of approximately one to two weeks [47]. The specific leaf area of the collected typical needles was derived using the same method applied in the woody-to-total area ratio estimation of the representative trees. Then, the litter collection LAI of each plot was obtained based on the specific leaf area and dry mass of all collected needles. Details of the determination of litter collection LAI for each plot can be found in the work of Zou at al. [47].

Gap Fraction
Large differences were observed between the p e (θ) or p w (θ) of the four optical methods in the five plots (Figures 1 and 2). For example, the difference between the p w (0) of DHP and MCI was as large as 0.45 in leaf-off plot 5 ( Figure 2). However, no large differences were observed between the p e (θ) or p w (θ) of TRAC and MCI even though the estimates of TRAC tended to be slightly smaller than those of MCI in most cases (Figures 1 and 2). We also found that the p e (θ) or p w (θ) of MCI was smaller than those of DHP and DCP in the five leaf-on or leaf-off plots (Figures 1 and 2). Similarly, the p e (0) or p w (0) of DHP tended to be larger than those of DCP in most cases (Figures 1 and 2). The differences between the p e (θ) or p w (θ) of the four optical methods tended to decrease with the increasing zenith angles due to the p e (θ) or p w (θ) estimates of the four methods were decreased with zenith angles obviously.
( ) or ( ) of TRAC and MCI even though the estimates of TRAC tended to be slightly smaller than those of MCI in most cases (Figures 1 and 2). We also found that the ( ) or ( ) of MCI was smaller than those of DHP and DCP in the five leaf-on or leaf-off plots (Figures 1 and 2). Similarly, the (0) or (0) of DHP tended to be larger than those of DCP in most cases ( Figures  1 and 2). The differences between the ( ) or ( ) of the four optical methods tended to decrease with the increasing zenith angles due to the ( ) or ( ) estimates of the four methods were decreased with zenith angles obviously.

Canopy Element and Woody Components Clumping Indices
Large differences were observed between the ( ) or ( ) obtained from DHP and MCI using the same or algorithm as well as from DHP or MCI using different or algorithms in the five leaf-on or -off plots (Figures 3 and 4). For example, the differences between the (60) of DHP and MCI in leaf-on plot 2 are 0.25, 0.16, and 0.05 for the three algorithms of CC, LX, and CLX, respectively ( Figure 3); the differences between the (60) of CC and LX as well as CC and CLX in leaf-on plot 2 are 0.36 and 0.47 for DHP, respectively ( Figure 3). However, the trend of the ( ) or ( ) obtained from TRAC was nearer to that of MCI derived using CC with a similar zenith angle in the five leaf-on or leaf-off plots, except leaf-on plot 5 in most circumstances (Figures 3 and 4). Similarly, no large differences were found between the (0) or (0) measurements obtained from DCP and MCI using CC in the five leaf-on plots, except plots 1, 2, and 5, or five leaf-off plots, except plots 1 and 2 (Figures 3 and 4). Interestingly, no large differences (<0.07)

Canopy Element and Woody Components Clumping Indices
Large differences were observed between the Ω e (θ) or Ω w (θ) obtained from DHP and MCI using the same Ω e or Ω w algorithm as well as from DHP or MCI using different Ω e or Ω w algorithms in the five leaf-on or -off plots (Figures 3 and 4). For example, the differences between the Ω e (60) of DHP and MCI in leaf-on plot 2 are 0.25, 0.16, and 0.05 for the three Ω e algorithms of CC, LX, and CLX, respectively ( Figure 3); the differences between the Ω e (60) of CC and LX as well as CC and CLX in leaf-on plot 2 are 0.36 and 0.47 for DHP, respectively ( Figure 3). However, the trend of the Ω e (θ) or Ω w (θ) obtained from TRAC was nearer to that of MCI derived using CC with a similar zenith angle in the five leaf-on or leaf-off plots, except leaf-on plot 5 in most circumstances (Figures 3 and 4). Similarly, no large differences were found between the Ω e (0) or Ω w (0) measurements obtained from DCP and MCI using CC in the five leaf-on plots, except plots 1, 2, and 5, or five leaf-off plots, except plots 1 and 2 ( Figures 3 and 4). Interestingly, no large differences (<0.07) were observed between the Ω e (θ) or Ω w (θ) obtained from MCI using LX and CLX in the five leaf-on or leaf-off plots, except leaf-off plot 5 (Figures 3  and 4). For DHP, CC tended to produce Ω e (θ) or Ω w (θ) estimates equal to 1 at the zenith angle range of approximately 40 • -90 • in the five leaf-on or leaf-off plots (Figures 3 and 4). That indicates that CC was not performed well in estimating the Ω e or Ω w of the five leaf-on or leaf-off plots from DHP at the zenith angle range of 40 • -90 • .

Canopy Element and Woody Components Clumping Index Algorithms
Tables 4-6 show that the LAI estimations of DHP and MCI were greatly affected by the or algorithm adopted in the estimation. Large differences were observed between the LAIs obtained from DHP or MCI using the same inversion model and woody components correction method but with different or algorithms in the five plots (Tables 4 and 6). For example, the differences in the LAIs obtained from DHP using the same inversion model and the woody components correction method but with CC and LX or CC and CLX ranged from 0.73 to 1.54 (25%-47%) when destructive

Canopy Element and Woody Components Clumping Index Algorithms
Tables 4-6 show that the LAI estimations of DHP and MCI were greatly affected by the Ω e or Ω w algorithm adopted in the estimation. Large differences were observed between the LAIs obtained from DHP or MCI using the same inversion model and woody components correction method but with different Ω e or Ω w algorithms in the five plots (Tables 4 and 6). For example, the differences in the LAIs obtained from DHP using the same inversion model and the woody components correction method but with CC and LX or CC and CLX ranged from 0.73 to 1.54 (25%-47%) when destructive woody-to-total area ratio measurements were adopted in the estimation ( Table 4). The differences between the RMSE and MAE of the LAI obtained from DHP using the same inversion model but with CC and LX or CC and CLX ranged from 0.94-1.17 (21%-25%) and 1-1.22 (23%-28%) when destructive woody-to-total area ratio measurements were used in the estimation (Table 5). By contrast, small differences were observed between the RMSE (0.01-0.24, 0%-5%) and MAE (0.03-0.23, 1%-5%) of the LAI obtained from DHP using the same inversion model and woody components correction method but with LX and CLX (Table 5). A similar trend of differences between the LAI obtained using the same inversion model and woody components correction method but with different Ω e or Ω w algorithms for DHP was also observed for MCI (Tables 6 and 7).
For DHP, CC tended to underestimate the ESU LAI of the five forests obviously compared with LX and CLX (Tables 4 and 5). LX and CLX outperformed CC in ESU LAI estimation with smaller RMSE and MAE for DHP regardless of the inversion model and woody components correction methods used in the estimation (Tables 4 and 5). Compared with DHP, CC outperformed LX and CLX in the ESU LAI estimation with the smallest MAE for MCI when the leaf-on MCI images were used in the estimation (Table 7). Table 4. LAI derived from the leaf-on DHP images using the three canopy element or woody components clumping index (Ω e and Ω w ) algorithms, three inversion models and destructive or MCI woody-to-total area ratio (α) or WAI obtained from the leaf-off DHP images in the five plots.  Table 5. Correlation statistics between the litter collection LAI and the LAI derived from the leaf-on DHP images using the three canopy element or woody components clumping index (Ω e or Ω w ) algorithms, three inversion models, and destructive or MCI woody-to-total area ratio (α) or WAI obtained from the leaf-off DHP images in the five plots (two-tailed Student's t-test with 95% confidence level). The root mean square error (RMSE) and mean absolute error (MAE) are expressed in LAI units (m 2 /m 2 ).  Table 6. LAI derived using the three canopy element or woody components clumping index (Ω e or Ω w ) algorithms and inversion models from the leaf-on MCI images or both the leaf-on and leaf-off MCI images in the five plots. When only the leaf-on MCI images were used in the LAI estimation, the WAIs of the five plots were obtained using Equations (8) and (11) based on the leaf-on MCI images. When both the leaf-on and leaf-off MCI images were used in the LAI estimation, the WAIs of the five plots were obtained using Equations (8) and (11) based on the leaf-off MCI images. Then, the LAI of the five plots was derived using Equation (16)

Inversion Model
Large differences were observed between the LAIs obtained from DHP using the same Ω e or Ω w algorithm and woody components correction method but with different inversion models in the five plots (Table 4). For example, the difference between the LAI of Miller and LAI-2200 for DHP was as large as 0.81 (27%) in plot 3 when CC and MCI woody-to-total area ratio were adopted in the estimation ( Table 4). The differences between the RMSE and MAE of the LAI obtained from DHP using the same Ω e algorithm and woody components correction method but with different inversion models were in the range of 0.0-0.45 (0%-10%) and 0.01-0.46 (0%-10%), respectively ( Table 5). The LAI of Miller tended to be larger than LAI-2200 or Beer for DHP in the five plots in most cases (Table 4). Miller is amongst the three inversion models to derive the LAI with the smallest MAE and RMSE for DHP (Table 5).
Similar to DHP, a large difference (1.61, 46%) was also observed for MCI between the LAI of Beer and MCI_0-85 in plot 4 when CC and the WAI obtained from the leaf-off MCI images were adopted in the estimation ( Table 6). The differences between the RMSE and MAE of the LAI obtained from MCI using the same Ω e or Ω w algorithm and woody components correction method but with different inversion models were in the range of 0.22-0.90 (5%-20%) and 0.19-0.89 (3%-20%) ( Table 7). For MCI, the best inversion model for estimating the ESU LAI of the five forests changed with the Ω e or Ω w algorithm and woody components correction method adopted in the estimation (Table 7).

Woody Component Correction Method
The LAIs obtained from DCP, DHP, and MCI in the five plots were remarkably affected by the woody components correction methods used in the estimation (Tables 4-9). For example, the difference was as large as 1.87 (49%) in plot 1 between the LAIs obtained from the leaf-on DHP images using CLX and Miller but with two woody components, correction methods of MCI woody-to-total area ratio and the WAI obtained from the leaf-off DHP images, respectively (Table 4). For DHP, large differences were observed between the RMSE or MAE of the LAIs obtained using the same Ω e or Ω w algorithm and inversion model but with two woody components correction methods of MCI woody-to-total area ratio and the WAI derived from the leaf-off DHP images in the range of 0.72-1.04 (16%-23%) or 0.74-1.14 (17%-27%) ( Table 5). However, small differences were found for DHP between the RMSE or MAE of the LAIs obtained using the same Ω e or Ω w algorithm and inversion model but with MCI and destructive woody-to-total area ratio values (Table 5). Like DHP, large differences were also observed for MCI between the RMSE or MAE of the LAIs obtained using the same Ω e or Ω w algorithm and inversion model but with the two woody components correction methods, ranging from 0.04-1.13 (1%-26%) ( Table 7). The LAI obtained from DHP or MCI or DCP using the woody components correction method of WAI derived from the leaf-off DHP or MCI or DCP images was obviously smaller than that derived using the same inversion model and Ω e or Ω w algorithm but with other woody components correction methods (Tables 4, 6, and 8). Table 8. LAI derived from DCP using two calculation schemes. One scheme calculates the LAI from the leaf-on DCP images (Equation (15)), and the other scheme calculates the LAI based on the PAI and WAI derived from the leaf-on and leaf-off DCP images, respectively (Equation (16)).

LAI Estimation Methods
Tables 4, 6, 8, and 10 indicate large differences between the LAI obtained using the same or similar inversion models, Ω e or Ω w algorithm and woody components correction method but with different optical methods in the five plots. For example, the difference between the LAI of DCP derived from the leaf-on DCP images and that of MCI obtained from the leaf-on MCI images using CC and Beer is 3.48 (71%) in plot 1 (Tables 6 and 8). Amongst the four optical methods, TRAC exhibited the best performance in terms of obtaining the ESU LAI of the five forests with the smallest RMSE (0.31, 6%) and MAE (0.24, 4%) followed by MCI (Tables 5, 7, 9 and Figure 5). For MCI, the ESU LAIs with the smallest RMSE (0.64, 14%) and MAE (0.43, 12%) were obtained from leaf-on MCI images using CC and Beer (Table 7). DHP tended to underestimate the ESU LAI of the five forests regardless of the inversion model, Ω e or Ω w algorithm and woody components correction method used in the estimation (Tables 4 and 5). For DHP, the ESU LAIs with the smallest RMSE (1.28; 28%) and MAE (1.06; 21%) were obtained using Miller, CLX and destructive woody-to-total area ratio ( Table 5). The trend of LAI underestimation for DHP in the five plots was also observed for DCP (Tables 8 and 9). MCI and TRAC did not systematically underestimate the LAI in the five plots (Tables 6 and 10). For plot 5 with a high litter collection LAI of 6.69, DHP, DCP, and MCI tended to underestimate the LAI (Tables 4, 6, and 8).

Which Canopy Element or Woody Component Clumping Index Algorithm(s) or Inversion Model(s) is (are) More Reliable to Be Adopted in the ESU LAI Estimation of L. principis-rupprechtii Forests from DHP and MCI?
Based on Tables 4-7, the ESU LAI estimation of DHP and MCI was largely affected by the canopy element or woody components clumping index ( or ) algorithm and inversion model adopted in the estimation. This finding is consistent with the conclusions drawn in previous studies [19,24,32,41]. As expected, large differences were found between the LAIs obtained from DHP or MCI using the same inversion model and woody components correction method but with different or algorithms or the same or algorithm and woody components correction method but with different inversion models (Tables 4 and 6). Three factors contributed to the large LAI differences: the obvious variation of ( ) or ( ) and (θ) or (θ) at the zenith angle range of 0°-90° (Figures 1-4); the large differences between the (θ) or (θ) derived using different algorithms (Figures 3 and 4) and the different zenith angle ranges covered by the four inversion models ( Table  2).
The best or algorithm for estimating the ESU LAI of the five forests changed with the optical methods. For DHP, CLX, and LX outperformed CC in obtaining LAIs with the smallest RMSE and MAE (Table 5). This finding is consistent with the finding of previous studies which reported that CLX and LX outperformed CC for DHP in estimating the or of the leaf-on or leaf-off forest plots [19,33,41]. For MCI, CC outperformed LX and CLX in LAI estimation with the smallest RMSE and MAE ( Table 7). The different conclusions on the best or algorithm for DHP and MCI can be attributed to the following reasons: (1) DHP has the disadvantage of insufficient sampling at the zenith and a coarse resolution of DHP images at large zenith angles [54]. By contrast, MCI has the advantage of sampling the canopy with high resolution at all zenith angles. (2). The "P" method described by Chen et al. [40], which was adopted by DHP to derive the mean element width values, tended to overestimate the mean element width values and result in or overestimation for CC. The good performance of CC in the ESU LAI estimation for MCI was further confirmed by the

Which Canopy Element or Woody Component Clumping Index Algorithm(s) or Inversion Model(s) Is (Are)
More Reliable to Be Adopted in the ESU LAI Estimation of L. principis-rupprechtii Forests from DHP and MCI?
Based on Tables 4-7, the ESU LAI estimation of DHP and MCI was largely affected by the canopy element or woody components clumping index (Ω e or Ω w ) algorithm and inversion model adopted in the estimation. This finding is consistent with the conclusions drawn in previous studies [19,24,32,41]. As expected, large differences were found between the LAIs obtained from DHP or MCI using the same inversion model and woody components correction method but with different Ω e or Ω w algorithms or the same Ω e or Ω w algorithm and woody components correction method but with different inversion models (Tables 4 and 6). Three factors contributed to the large LAI differences: the obvious variation of p e (θ) or p w (θ) and Ω e (θ) or Ω w (θ) at the zenith angle range of 0 • -90 • (Figures 1-4); the large differences between the Ω e (θ) or Ω w (θ) derived using different algorithms (Figures 3 and 4) and the different zenith angle ranges covered by the four inversion models ( Table 2).
The best Ω e or Ω w algorithm for estimating the ESU LAI of the five forests changed with the optical methods. For DHP, CLX, and LX outperformed CC in obtaining LAIs with the smallest RMSE and MAE (Table 5). This finding is consistent with the finding of previous studies which reported that CLX and LX outperformed CC for DHP in estimating the Ω e or Ω w of the leaf-on or leaf-off forest plots [19,33,41]. For MCI, CC outperformed LX and CLX in LAI estimation with the smallest RMSE and MAE ( Table 7). The different conclusions on the best Ω e or Ω w algorithm for DHP and MCI can be attributed to the following reasons: (1) DHP has the disadvantage of insufficient sampling at the zenith and a coarse resolution of DHP images at large zenith angles [54]. By contrast, MCI has the advantage of sampling the canopy with high resolution at all zenith angles. (2). The "P" method described by Chen et al. [40], which was adopted by DHP to derive the mean element width values, tended to overestimate the mean element width values and result in Ω e or Ω w overestimation for CC. The good performance of CC in the ESU LAI estimation for MCI was further confirmed by the small RMSE and MAE of the LAI obtained from TRAC ( Figure 5). This point illustrates that when the canopy was sufficiently sampled and an accurate mean element width value was obtained, CC could provide relatively more accurate Ω e or Ω w estimates than those of LX and CLX. CC has two advantages in deriving Ω e or Ω w over LX and CLX. Firstly, CC requires no spatial distribution assumption of canopy element. By contrast, an assumption of the randomly distribution of canopy element at the segment scale with a length that is 10 times mean element width value was made for LX or CLX [37]. This assumption seems to almost hold true for MCI based on the close agreement between the Ω e or Ω w obtained from LX and CLX in most cases (Figures 3 and 4). The agreement illustrates that the canopy element at the segment scale approached random distribution. However, a reasonable guideline for defining the segment lengths of LX and CLX for DHP was unavailable. Obvious differences were found between the Ω e or Ω w of LX or CLX obtained with varying segment sizes [33,54]. Secondly, the large between-crown gaps, which dominated the clumping effects of forest canopies [39,40], can be effectively captured by CC, whereas LX or CLX was not qualified to detect and utilize such kinds of large gaps to derive Ω e or Ω w due to the small segment size.
The best inversion model for estimating the ESU LAI of the five forests was different between DHP and MCI (Tables 4-7). For DHP, Miller outperformed Beer and LAI-2200 in the ESU LAI estimation with the smallest RMSE and MAE (Table 5). For MCI, Beer outperformed MCI_0-85 in the LAI estimation with the smallest RMSE and MAE when CC was used in the estimation; by contrast, MCI_0-85 outperformed Beer in the LAI estimation when LX or CLX was adopted in the estimation (Table 7). This finding is consistent with the finding of Zou et al. [19] who reported that the best inversion model for estimating the PAI and WAI of leaf-on and leaf-off forest plots from optical methods is the function of the Ω e or Ω w estimation algorithm, PAI, LAI and WAI and the plant functional types of forests. When CC was adopted in the ESU LAI estimation, the best inversion model was expected with differences between DHP and MCI as the contrasting performance of CC in estimating the Ω e or Ω w of the five forests for DHP and MCI (Figures 3 and 4).

Which Woody Components Correction Method(s) Is (Are) Better in ESU LAI Estimation for the Four Optical Methods?
Destructive woody-to-total area ratio was the best woody components correction method amongst the three woody components correction methods in obtaining the LAI with the smallest RMSE and MAE, followed by MCI woody-to-total area ratio (Tables 5, 7, and 9). By contrast, the woody components correction method of obtaining WAI from leaf-off DHP, MCI and DCP images was the worst method for deriving an LAI with the largest RMSE and MAE in most cases (Tables 5, 7, and 9). As expected, destructive woody-to-total area ratio was the best woody components correction method as they were obtained based on the destructive woody-to-total area ratio measurements of representative trees of each plot, which were usually regarded as the highly accurate woody-to-total area ratio estimates compared with those obtained from indirect methods [14,23,55]. The small differences between the RMSE or MAE of the LAIs obtained from DHP using the same Ω e or Ω w algorithm and inversion model but with MCI and destructive woody-to-total area ratio measurements (Table 5) illustrate that MCI is an effective and relatively accurate solution to correct the overestimation of the woody components on the LAI estimation from DHP. Compared with the destructive method, MCI has the advantage of low cost, high efficiency, and non-destructiveness to canopies. Thus, MCI can be adopted as the routine option for deriving the woody-to-total area ratio values of L. principis-rupprechtii forests.
The PAI obtained from the leaf-on DHP images only using Miller and the three Ω e algorithms was systematically smaller than those obtained by summing the LAIs derived from the leaf-on DHP images using the same inversion model and Ω e algorithm with MCI woody-to-total area ratio and the WAI obtained from the leaf-off DHP images using the same inversion model and Ω w algorithm ( Figure 6). This finding is consistent with the finding of Calders [56], although different tree species were covered in the two studies. The differences between the derived PAI of the two calculation methods were mainly attributed to the preferential shading of woody components by the shoots in the canopies. The preferential shading of woody components by shoots would increase the p w (θ) and reduce the WAI and PAI obtained from the leaf-on DHP images as they should be [56]. Zou et al. [23] reported that the proportions of woody components shaded by shoots were large at the zenith angles ranging from 0 • to 80 • with a 10 • interval (7%-73%) in the five leaf-on plots of this study. Therefore, if the WAI derived from the leaf-off DHP images was used as the woody components correction method, then the contribution of woody components to the LAI estimation would be overcorrected as the PAI, which was obtained from the leaf-on DHP images containing a certain degree of WAI underestimation due to the preferential shading of woody components by the shoots in the canopies. This condition also explains why the LAIs obtained from leaf-on and leaf-off measurements were obviously smaller than those derived from the leaf-on measurements regardless of the optical methods, inversion models, and Ω e or Ω w algorithm used in the estimation (Tables 4, 6, and 8). Therefore, the WAI obtained from the leaf-off measurements is not recommended as the woody components correction method in the LAI estimation. Caution is required when adopting the WAI obtained from the leaf-off measurements as the woody components correction method in the LAI estimation even this method was frequently adopted in previous studies, such as Cutini et al. [57], Breda [17], Ryu et al. [35], and Toda et al. [58].
Forests. 2019, 10, x FOR PEER REVIEW 19 of 28 the canopies. The preferential shading of woody components by shoots would increase the ( ) and reduce the WAI and PAI obtained from the leaf-on DHP images as they should be [56]. Zou et al. [23] reported that the proportions of woody components shaded by shoots were large at the zenith angles ranging from 0° to 80° with a 10° interval (7%-73%) in the five leaf-on plots of this study. Therefore, if the WAI derived from the leaf-off DHP images was used as the woody components correction method, then the contribution of woody components to the LAI estimation would be overcorrected as the PAI, which was obtained from the leaf-on DHP images containing a certain degree of WAI underestimation due to the preferential shading of woody components by the shoots in the canopies. This condition also explains why the LAIs obtained from leaf-on and leaf-off measurements were obviously smaller than those derived from the leaf-on measurements regardless of the optical methods, inversion models, and or algorithm used in the estimation (Tables 4, 6, and 8). Therefore, the WAI obtained from the leaf-off measurements is not recommended as the woody components correction method in the LAI estimation. Caution is required when adopting the WAI obtained from the leaf-off measurements as the woody components correction method in the LAI estimation even this method was frequently adopted in previous studies, such as Cutini et al. [57], Breda [17], Ryu et al. [35], and Toda et al. [58].

Which Optical Method(s) Is (Are) More Reliable to Obtain the LAI of L. principis-rupprechtii Forest Plots?
On the basis of TAbles 5,7,9 and Figure 5, TRAC and MCI are recommended for use in obtaining the ESU LAI of L. principis-rupprechtii forests followed by DHP. DCP is not recommended because it results in obvious LAI underestimation. Given that DCP can sample the canopies sufficiently due to its high image resolution and the large sample size of the sampling scheme (16), accurate p e (0) or p w (0) and gap measurements can be obtained from DCP. Therefore, the assumption of the spherical angle distribution of the canopy element made in this study would be the key reason for the LAI underestimation of DCP. On the other hand, the LAI underestimation for DCP indicates that the canopy element and woody components projection coefficients of L. principis-rupprechtii forests at 0 • obviously deviated from 0.5, which was also observed for forests with other tree species [59,60]. Macfarlane [41] and Chianucci et al. [34] reported that accurate ESU LAI estimates can be obtained from DCP in broad-leaved forests because of two reasons. Firstly, the leaf angle distribution of the forests in the studies of Macfarlane [41] and Chianucci et al. [34] approached spherical distribution. Secondly, the estimates derived from DCP in the two studies were PAIs and not LAIs. Good agreement was found between the PAI and LAI obtained through the litter collection method or allometric equation in the two studies [34,41]. Therefore, if woody components were considered in the LAI estimation, then DCP would underestimate the ESU LAI of forests in the two studies [34,41]. Thus, when the field collected leaf or shoot angle distributions are not available, caution is needed if DCP is adopted in estimating the ESU LAI of forests.
The p e (θ) or p w (θ) of MCI tends to be smaller than those of DHP and DCP in the five plots (Figures 1 and 2) because of two reasons. Firstly, the isolated shoots and branches in the visible band of the MCI or DHP images tended to be easily overexposed due to the high contrast between the canopy element and the sky (Figure 7a). However, the contrast between the canopy element and the sky or cloud in the near-infrared band of MCI is relatively small (Figure 7b) [61]. Secondly, the shoots have a low reflectance in the visible band due to the strong absorption by chlorophyll, a relatively high reflectance in the near-infrared band due to the internal leaf scattering, and no absorption [62]. Consequently, the shoot sizes in the near-infrared band were larger than those in the visible band for the same shoot (Figure 7a,b). This issue was not carefully considered by the "IsoData" classification algorithm adopted in the MCI image classification.
MCI was originally designed to obtain the woody-to-total area ratio values of forests [14,23]. Compared with DHP, MCI has two advantages in the ESU LAI estimation of forests. Firstly, the ESU LAI of forests can be derived from this single method without the combination of any other optical method due to its ability to discriminate between leaves or shoots, woody components, and sky [14]. Secondly, the forest canopies can be sufficiently sampled by MCI at all zenith directions of the hemisphere due to its high image resolution and hemispherical sampling scheme [14]. MCI also has two disadvantages in ESU LAI estimation. Firstly, MCI is a proprietary device and is unavailable as a commercial instrument for scientists compared with LAI-2200 and TRAC. Secondly, great effort and a large amount of resources are needed to collect and process the MCI images. For example, 54 or 60 measurements (6 measurements for each zenithal direction) and 1 measurement must be collected at each sampling point for MCI [14,23] and DHP, respectively. It is encouraging that the two disadvantages of MCI can be overcome as follows: (1) The charge couple devices of consumer-level cameras have an imaging capability in the near-infrared band after removing the infrared blocking filter before the charge couple devices [63]. Then, consumer-level cameras can be modified following the procedures described in Chapman [63] to obtain the visible and near-infrared band images of forests similar to MCI with the combination of filters. (2) Given that relatively accurate woody-to-total area ratio measurements [23] and LAI can be derived from MCI using Beer (Table 7), the Beer inversion model can be used to derive the ESU LAI of forests, and then only those MCI images with a zenith angle near 57 • should be collected to reduce the effort and resources needed for MCI.  Figure 7c shows the corresponding classified MCI image of the MCI image pair.
The ( ) or ( ) of TRAC tended to be smaller than those of MCI, DHP, and DCP ( Figures  1 and 2). This finding is consistent with the finding of Raabe et al. [64], Leblanc et al. [32], Pisek et al. [33], and Ryu et al. [35]. This can be attributed to the limitation of TRAC sensor resolution as the small within canopy gaps were underestimated by TRAC [64]. For TRAC, a problem occurred when the LAI-2200 or Miller inversion model was adopted by TRAC in the ESU LAI estimation because only one ( ) or ( ) and ( ) or ( ) estimate can be obtained from each TRAC measurement. In practice, the TRAC measurement datasets of each plot with the same zenith angle range covered by LAI-2200 or Miller are difficult to gather, especially for plots in high latitude areas where the solar zenith angles of daytime are usually small. An alternative option for TRAC to be adopted as the routine method in the LAI estimation is to derive the LAI with the combination of TRAC and Beer. Similar estimates were obtained from MCI and TRAC (0.90 and 0.94) at a zenith angle of 64° in the leaf-off plot 3, indicating that TRAC remains effective in deriving at this zenith angle (Figure The p e (θ) or p w (θ) of TRAC tended to be smaller than those of MCI, DHP, and DCP (Figures 1  and 2). This finding is consistent with the finding of Raabe et al. [64], Leblanc et al. [32], Pisek et al. [33], and Ryu et al. [35]. This can be attributed to the limitation of TRAC sensor resolution as the small within canopy gaps were underestimated by TRAC [64]. For TRAC, a problem occurred when the LAI-2200 or Miller inversion model was adopted by TRAC in the ESU LAI estimation because only one p e (θ) or p w (θ) and Ω e (θ) or Ω w (θ) estimate can be obtained from each TRAC measurement. In practice, the TRAC measurement datasets of each plot with the same zenith angle range covered by LAI-2200 or Miller are difficult to gather, especially for plots in high latitude areas where the solar zenith angles of daytime are usually small. An alternative option for TRAC to be adopted as the routine method in the LAI estimation is to derive the LAI with the combination of TRAC and Beer. Similar Ω w estimates were obtained from MCI and TRAC (0.90 and 0.94) at a zenith angle of 64 • in the leaf-off plot 3, indicating that TRAC remains effective in deriving Ω w at this zenith angle ( Figure 8). However, Law et al. [65] suggested avoiding taking TRAC measurements at zenith angles larger than 60 • due to the difficulty of TRAC to distinguish small gaps at these zenith angles. An error Ω e estimate of 0.75 was derived at zenith angle of 75 • for TRAC in the leaf-on plot 5 as all the photosynthetic photon flux density values of this TRAC measurement was equal to zero µmol m −2 s −1 (not shown in Figure 4). Therefore, caution is needed for TRAC to derive the Ω e or Ω w and LAI of L. principis-rupprechtii forests at zenith angles larger than 70 • .
Forests. 2019, 10, x FOR PEER REVIEW 22 of 28 8). However, Law et al. [65] suggested avoiding taking TRAC measurements at zenith angles larger than 60° due to the difficulty of TRAC to distinguish small gaps at these zenith angles. An error estimate of 0.75 was derived at zenith angle of 75° for TRAC in the leaf-on plot 5 as all the photosynthetic photon flux density values of this TRAC measurement was equal to zero µ mol m -2 s -1 (not shown in Figure 4). Therefore, caution is needed for TRAC to derive the or and LAI of L. principis-rupprechtii forests at zenith angles larger than 70°. Two reasons contributed to the ESU LAI underestimation for DHP in the five forests [47]. Firstly, previous studies reported that the gap fraction measurements of optical methods tended to be saturated at forests with large LAIs due to canopy closure, especially for forests with LAIs ~ > 5.0-6.0 [16,47,55]. Secondly, the ( ) or ( ) of DHP was larger than that of MCI in the five leaf-on or leaf-off plots in most cases (Figures 3 and 4), indicating that DHP overestimated the or estimates in the five leaf-on or leaf-off plots because the or estimation of DHP suffered from the two error sources of canopy sampling and mean element width estimation.
Compared with TRAC, MCI, and DCP, DHP has the advantage of obtaining the hemispherical direction measurements of ( ) or ( ) and ( ) or ( ) simultaneously from single or several DHP images. Therefore, less effort and resources are required for DHP than the other three optical methods in the ESU LAI estimation of forests. DHP is frequently used as the routine method for obtaining the ESU LAI of forests. However, the performance of DHP in deriving the ESU LAI of the five forests was strongly affected by the inversion model, or algorithm and woody components correction method adopted in the estimation (Tables 4 and 5). This finding is consistent with the finding of Liu et al. [24], Zou et al. [19], and van Gardingen et al. [66] who reported large differences between the LAIs derived from DHP using various inversion models, or algorithms and woody components correction methods. Therefore, the best combination of inversion model, or algorithm and woody components correction method should be identified before DHP is selected as the routine method to estimate the ESU LAI of forests. Table 5 indicates that the LAI of DHP with the smallest MAE (21%) was obtained using Miller, CLX and destructive woody-to-total area ratio. This result indicates that DHP is not qualified to derive the ESU LAI of L. principis-rupprechtii forests with estimation errors of <5% or 20% required by GCOS. However, the LAI with the MAEs <20% can be obtained from DHP in the five L. principisrupprechtii forests using Miller, CLX and destructive woody-to-total area ratio by optimizing the sampling schemes reported by Zou et al. [47]. An LAI MAE range of 11%-16.4% was obtained from DHP when Miller, CLX, and destructive woody-to-total area ratio in addition to two sampling schemes with the sample sizes in the range of 3-9 were adopted in the ESU LAI estimation of the five Two reasons contributed to the ESU LAI underestimation for DHP in the five forests [47]. Firstly, previous studies reported that the gap fraction measurements of optical methods tended to be saturated at forests with large LAIs due to canopy closure, especially for forests with LAIs~> 5.0-6.0 [16,47,55]. Secondly, the Ω e (θ) or Ω w (θ) of DHP was larger than that of MCI in the five leaf-on or leaf-off plots in most cases (Figures 3 and 4), indicating that DHP overestimated the Ω e or Ω w estimates in the five leaf-on or leaf-off plots because the Ω e or Ω w estimation of DHP suffered from the two error sources of canopy sampling and mean element width estimation.

Do the Four Optical Methods Qualify to Obtain the ESU LAI of Forests with the Accuracy Match the Requirements of GCOS?
Compared with TRAC, MCI, and DCP, DHP has the advantage of obtaining the hemispherical direction measurements of p e (θ) or p w (θ) and Ω e (θ) or Ω w (θ) simultaneously from single or several DHP images. Therefore, less effort and resources are required for DHP than the other three optical methods in the ESU LAI estimation of forests. DHP is frequently used as the routine method for obtaining the ESU LAI of forests. However, the performance of DHP in deriving the ESU LAI of the five forests was strongly affected by the inversion model, Ω e or Ω w algorithm and woody components correction method adopted in the estimation (Tables 4 and 5). This finding is consistent with the finding of Liu et al. [24], Zou et al. [19], and van Gardingen et al. [66] who reported large differences between the LAIs derived from DHP using various inversion models, Ω e or Ω w algorithms and woody components correction methods. Therefore, the best combination of inversion model, Ω e or Ω w algorithm and woody components correction method should be identified before DHP is selected as the routine method to estimate the ESU LAI of forests.

Do the Four Optical Methods
Qualify to Obtain the ESU LAI of Forests with the Accuracy Match the Requirements of GCOS? Table 5 indicates that the LAI of DHP with the smallest MAE (21%) was obtained using Miller, CLX and destructive woody-to-total area ratio. This result indicates that DHP is not qualified to derive the ESU LAI of L. principis-rupprechtii forests with estimation errors of <5% or 20% required by GCOS. However, the LAI with the MAEs <20% can be obtained from DHP in the five L. principis-rupprechtii forests using Miller, CLX and destructive woody-to-total area ratio by optimizing the sampling schemes reported by Zou et al. [47]. An LAI MAE range of 11%-16.4% was obtained from DHP when Miller, CLX, and destructive woody-to-total area ratio in addition to two sampling schemes with the sample sizes in the range of 3-9 were adopted in the ESU LAI estimation of the five plots except plot 5 [47]. This point illustrates that DHP is qualified to derive the ESU LAI of L. principis-rupprechtii forests with estimation errors of <20% if the inversion model, Ω e or Ω w algorithm and woody components correction method in addition to sampling scheme are considered in the estimation. The LAI MAE range of 11%-16.4% is near that in the study by Leblanc and Fournier [20] who reported a minimum MAE of 11% for the PAI derived from DHP when an appropriate inversion model and Ω e algorithm are adopted in the PAI estimation of the simulated forest scenes. The LAI MAE range of 11%-16.4% is larger than the minimum PAI MAE in the study of Leblanc and Fournier [20] because the LAI estimation from the field-collected DHP images suffered from three additional estimation errors (i.e., observation conditions, DHP image classification, and woody components correction method) compared with those in the simulation study [47]. Table 5 indicates that DHP is not qualified to derive the ESU LAI of L. principis-rupprechtii forests with estimation errors of <5%. Table 9 shows that DCP is not qualified to derive the ESU LAI of L. principis-rupprechtii forests with estimation errors of <5% or 20%. Unlike DCP and DHP, MCI and TRAC are qualified to derive the ESU LAI of L. principis-rupprechtii forests with estimation errors of <20% (Table 7 and Figure 5). Moreover, TRAC shows the potential in obtaining the ESU LAI of L. principis-rupprechtii forests with estimation errors of <5% ( Figure 5).

Conclusions
In this study, the performance of four optical methods (i.e., DHP, DCP, TRAC, and MCI) in estimating the ESU LAI of L. principis-rupprechtii forests was evaluated using the LAIs obtained from litter collection measurements. The impact of three factors (i.e., inversion model, Ω e or Ω w algorithm, and woody components correction method) on the LAI estimation was analyzed. Results show that the ESU LAI estimation was largely affected by the three factors. To obtain accurate ESU LAIs of L. principis-rupprechtii forests, we suggest the following procedures: (1) using MCI and TRAC as the ESU LAI estimation methods; (2) using CC to derive the Ω e or Ω w ; (3) using the destructive or MCI woody-to-total area ratio as the woody components correction method; (4) using the Beer inversion model to derive the ESU LAI.
The accuracies of ESU LAIs obtained by the four optical methods were evaluated in terms of whether they matched the LAI accuracy target required by GCOS. Results show that the four optical methods, except for DCP, could obtain the ESU LAI of L. principis-rupprechtii forests with an MAE of <20% required by GCOS. Only TRAC shows potential in obtaining the ESU LAI of L. principis-rupprechtii forests with an MAE of <5%.
Given the limited efforts and resources available, only five plots were covered in this study. More plots should be covered in the future to improve the conclusions drawn in this study and evaluate whether management activities (e.g., branches harvesting and thinning) and the LAI range of the plots are factors that affect the performance of the four optical methods for deriving the ESU LAI of L. principis-rupprechtii forest. Future work could include efforts to evaluate the performance of a commonly used method of LAI-2200 in the ESU LAI estimation of L. principis-rupprechtii forests.