Estimating the Horizontal and Vertical Distributions of Pigments in Canopies of Ginkgo Plantation Based on UAV-Borne LiDAR, Hyperspectral Data by Coupling PROSAIL Model

: Pigments are the biochemical material basis for energy and material exchange between vegetation and the management and precision silvicuture for plantations.


Introduction
As an essential part of the forest ecosystem, planted forests play a significant role in maintaining the ecological environment, conserving biodiversity, and providing a variety of biochemical products for the ecosystem [1]. According to the China Forest Resources Report (2009)(2010)(2011)(2012)(2013), the area of planted forests is approximately 7.95 × 10 7 km 2 , accounting for approximately 32.94% of the national forest area [2]. Ginkgo (Ginkgo biloba L.) is a unique multi-purpose tree species in China, with important economic and ecological value [3]. As the principal component of the forest's biochemical traits, total chlorophylls (Cab) and carotenoids (Car) can reflect the growth status of trees and environmental stress, and can data with the PROSAIL model for estimating pigments in individual Ginkgo trees with optimal LUTs, and (3) to characterize the horizontal and vertical distribution characteristics of pigments for Ginkgo plantations at different ages.

Study Area and Workflow
As a key coastal planted forest, Dongtai Forest (120 • 49 32.2 E, 32 • 52 20.6 N) is an essential multi-purpose forest plantation for wood production, marine disaster prevention, etc. The field covers approximately 2239 km 2 , and the forest coverage rate is around 85%. The vegetation types are mainly deciduous broad-leaved and needle-leaved forests, including mainly Ginkgo (Ginkgo biloba L.), Poplar Populus deltoides), and Dawn redwood (Metasequoia glyptostroboides). The Dongtai forest is located in a flat terrain, with an elevation range of between approximately 8-20 m. The annual average temperature and annual rainfall are approximately 14.6 • C and 1050 mm, respectively [45].
Within the forest stands of different ages (i.e., 13-year-old (young-aged) and 22-yearold (middle-aged)) and planting densities of Ginkgo trees, LiDAR and hyperspectral data were obtained separately including two UAV flight stripes (i.e., a LiDAR point cloud strip and a hyperspectral data strip) on 17-18 Aug. 2019. A total of seven sample trees were chosen for the study of biochemical traits' extraction: (1)   This study coupled UAV-borne LiDAR and hyperspectral imaging data with the PROSAIL model to map the horizontal and vertical distributions of pigments within the canopies of Ginkgo plantations. An overview of the workflow employed in this study is shown in Figure 2 and is mainly divided into three parts: (1) field data collection and measurement; (2) remote sensing data acquisition, processing and fusion; and (3) multi-LUTs generation from PROSAIL model and LUT-based inversion. An overview of the analysis workflow for assessing the content of ginkgo pigments. Note: "N, S" represent the north and south orientation; "upper, middle, lower" represent the upper, the middle, and the lower vertical layer; "2 age-groups" represent age 13-year-old and age 22-year-old; "5 entries" represent the number of LUT entries, varying from 1000 to 16,000 with the step of multiplying as 2 (i.e., 1000, 2000, 4000, 8000, 16,000); "5 random repeats" represent running random parameter sampling 5 times; "30 random repeats" represent running random parameter sampling 30 time; Cab = Chlorophyll; Car = Carotenoid; LAI = Leaf area index; ALA = Average leaf angle; N = Structural coefficient; Cw = Water content; Cm = Dry matter content.

Field Data Collection and Measurement
The leaf samples were collected according to different heights and orientations (See Figure 2). According to the crown height of each sampling tree, the crown was divided into three equal vertical layers, namely, the upper, middle and lower canopy, respectively. Furthermore, each vertical layer was further separated into the northern crown and southern crown, thus a total of six parts were separated within the crown for analysis (3 vertical layers × 2 orientations: upper north, upper south, middle north, middle south, lower north and lower south). Seven sampling trees were divided into 42-part samples (7 sampling trees × 6 parts). In each part, approximately 100 leaves of the outside crown were collected using branch shears. For each part, 100 leaves were divided into two piles randomly (ten leaves: the pigment content was measured; other leaves: the water and dry matter content was measured). These collected leaves of 42 parts were placed separately on the sealed bag with labels and put in a box with ice, then sent to the laboratory for the extraction of pigment content. Pigment content of each part was measured by the organic solvent extraction method in the laboratory with ten leaves. We punched with a disk of 1.48 cm 2 , cut these disks into pieces and mixed then weighted 0.1 g pieces. These pieces were dipped in 10 mL of 96% ethanol to extract pigments until they turned completely white. The Specord ® 200 plus UV/Vis spectrophotometer [46] was used to measure the absorption values of chlorophyll a (Chla), chlorophyll b (Chlb) and carotenoids (Car) at 470, 649, and 665 nm wavelengths, respectively. According to the area of the disk (1.48 cm 2 ), the specific leaf area (SLA) of each leaf sample was measured. Then, the mass-based Cab (mg·g −1 ) and Car (mg·g −1 ) were transformed to area-based Cab (µg·cm −2 ) and Car (µg·cm −2 ) content according to SLA.

Remote Sensing Data Acquisition and Processing
LiDAR data ( Figure 3a) and hyperspectral images ( Figure 3b) were acquired by the GV1300 multi-rotor UAS platform on 17 and 18 Aug. 2019, respectively. The flight parameters and UAV-LiDAR system properties are shown in Table 1.  The LiDAR data were captured using a Velodyne Puck VLP-16 scanner at an elevation of 80 m above ground level with a flight speed of 6 m·s −1 on 18 Aug. The obtained LiDAR original data were processed by point cloud denoising, filtering, classification and normalization. Combining the IMU data with GPS data recorded by the base station and UAV platform, point coordinates were calculated. The digital surface model (DSM) and digital elevation model (DEM) with a resolution of 0.18 m were generated by the irregular triangulation network (TIN) algorithm for the pre-processing and hyperspectral images co-registration. The processing LiDAR point cloud adopted the point cloud segmentation method to segment the seven sampling trees and obtain the corresponding six part LiDAR point clouds for each sampling tree. In addition, the boundaries corresponding to each part-point-clouds for seven sampling tree were determined by manual delineation.
Hyperspectral images with 0.18 m spatial resolution of the 13-year-old and the 22-yearold Ginkgo plantations were obtained at 10.00 a.m. and 16.00 p.m. on Aug. 17, respectively. The images (spectral range = 403-929 nm, spectral resolution = 2.3 nm) were obtained by the ZK-VNIR-FPG480 imager at 80 m altitude above ground level with a flight speed of 4.8 m·s −1 . A standard whiteboard was used to convert image brightness value (DN value) to spectral reflectance from the imagery. Then, the Savitzky-Golay filter (S-G filter), a weighted moving average filtering system, was used for reflectance image smoothing and denoising [47]. Specifically, S-G filtering was carried out in an IDL environment with a frame size of 5 data points (2nd degree polynomial) [48]. Due to the influence of random noise in the process of hyperspectral image collection, the 172 bands within the spectral range of 450 nm to 850 nm were retained for further analysis.
Before carrying out the remote sensing data fusion, the pre-processing hyperspectral images co-registration was carried out using the DSM created from the processing of the LiDAR point cloud. To make sure the fusion performance of the hyperspectral images and the processing LiDAR point cloud had the highest spatial consistency, in the co-registration two ground control points were randomly selected for each area of 5 × 5 m 2 from both the hyperspectral images and the DSM data. The accuracies of co-registration for hyperspectral images were both less than 9 cm (half a hyperspectral image pixel).

Remote Sensing Data Fusion
The DSM-based approach, proposed by Shen et al. (2020) [45], was adopted to fuse the processing LiDAR point cloud and co-registered hyperspectral images. This fusion approach can accommodate the lack of spectral information in LiDAR and the lack of canopy structure information in hyperspectral data. Firstly, the LiDAR point clouds were divided into columns (0.18 m × 0.18 m × 30 m) consistent with the pixel grid of the hyperspectral images in the horizontal orientation. Secondly, for each column, the highest LiDAR point cloud was retained, while excluding other point clouds. Finally, the hyperspectral reflectance value of the corresponding pixel grid was distributed to the highest LiDAR point cloud within each column. Thus, the fused hyperspectral point clouds were obtained, containing spectra and structural information. The reflectance of fused hyperspectral point clouds within six parts for each sampling tree was simplified as the point-based reflectance.

Extracting Profile Characteristic Variables
Referring to the previous studies [49,50], we calculated the Foliage Profile (FP) to describe the vertical distribution of branches and leaves. FP was extracted from the processing of the LiDAR point cloud [51,52]. FP is expressed as the sum of the branch and leaf area of the unit horizontal area within the unit volume at canopy height (i.e., 0.3 m) above ground level [50].

PROSAIL Model
In this study, the LUT-based-inversion method was adopted for the estimation of Cab and Car content. The PROSAIL model was used to generate multiple LUTs (multi-LUTs), which are query tables consisting of multiple model parameters' combinations and their corresponding simulated reflectance value. Then, the estimation of Cab and Car content was implemented with the above multi-LUTs [53].

Local Sensitivity Analysis
In this study, following the procedure in the previous study [54], local sensitivity analysis (LSA) was conducted to evaluate the relative importance of the model input parameters. These parameters were adjusted for value and combination to generate the corresponding modified reflectance curves, compared with the actual reflectance curves. According to the deviation of the comparison of curves, different parameter effects were determined.
The numerical range of model input parameters was mainly set in regard to field measured data, previous relevant studies and actual field flight situation [55,56]. There was a total of 15 PROSAIL model input parameters (e.g., structural coefficient (N), etc.). Given that eight parameters have significant effects on reflectance simulation within the visible-near infrared region [54,57], these eight were selected for sensitivity analysis within 450-850 nm. The numerical range of analysis parameters is shown in Table 2. First, eight parameter values were sampled seven times to gain the same step. Second, we evaluated the effort of model parameters on the modified reflectance within the same unit step.  Figure 4 shows the results of LSA. In the region of visible and near-infrared, LAI, ALA, Cab and Car cause a large variation in the spectral reflectance, while those of N, hspot, and psoil are negligible. In particular, skyl had negligible effect on the variation of the spectral reflectance. In particular, the sensitive bands of Cab and Car are concentrated in the range of 525-750 nm and 485-550 nm, respectively. In particular, LAI and ALA had a dominant influence on the whole of the spectral domains between 450-850 nm. However, the actual values of LAI and ALA were not measured in the field data collection. Given that nearly 45% of the Ginkgo branch angels were concentrated between 36-55 • [58], ALA was fixed at 45 • . Meanwhile, the range of LAI values was 0.3-3 m 2 ·m −2 for Ginkgo trees [59]. Considering that LAI have a dominating influence on the whole of the spectral domains, uniform distribution was adopted for the LAI sampling. Similarly, random distribution was adopted for the Cab sampling. Regarding the sampling of Car, we followed the procedure in Li et al. (2020) [60], adopting the relationship between Cab and Car to generate the value of Car. The values of other parameters were fixed according to the mean of the measured data (Table 3).

LUT-Based-Inversion and Pigment Estimations
In this study, the LUT-based-inversion method was adopted for the estimation of Cab and Car content. The PROSAIL model was used to generate multiple LUTs (multi-LUTs), which are query tables consisting of multiple model parameter combinations and their corresponding simulated reflectance value. Then, the estimation of Cab and Car content was implemented with the above multi-LUTs [63]. In particular, to maintain consistency between simulated reflectance spectra derived from the PROSAIL model (at 1 nm interval) and the hyperspectral imaging data, the spectral interval for simulated reflectance spectra in the LUTs were resampled to actual wavelengths of hyperspectral imaging data (around 2.3 nm intervals).
The LUTs used were divided into two age-groups, 13-year-old LUTs and 22-year-old LUTs. Multi-LUTs were constructed according to different calibrated parameters set at different ages. Following the procedure for establishing multi-LUTs by Li et al. (2019) [63], multi-LUTs were constructed from two aspects: (1) to determine the suitable LUT entries of pigments for two ages using LUT-database-1; and (2) to find the optimal LUT of pigments for two ages using LUT-database-2.
LUT-database-1 was composed of multi-LUTs with combinations of five entries (i.e., number of entries as 1000, 2000, 4000, 8000 and 16,000). In particular, repetition of independent random sampling for the LUT-database-1 was undertaken five times to avoid the influence of sampling effects (i.e., 5 random repeats). Therefore, the LUT-database-1 included 50 LUTs (2 age-groups × 5 entries × 5 random repeats). The LUT-database-2 consisted of two age-groups and 30 x random repeats with the suitable number of entries determined from the LUT-database-1, yielding a total of 60 LUTs (2 age-groups × 30 random repeats).
The LUT-based-inversion adopted ARTMO (Automated Radiative Transfer Model Operator), a scientific software package created by   [64]. This package can be obtained on the website https://artmotoolbox.com/ (10 June 2021). Two cost functions (i.e., RMSE and Cressie-Read functions) were selected for inversion. The general definition of the cost function is shown in Equation (1): where COST(R) is the formula of R; R is the reflectance value; R Measured is the measured mean reflectance; and R LUT is the simulated reflectance of the model.
The general performance of LUT-based-inversion was evaluated by coefficient of determination (R 2 ), root mean square root (RMSE) and normalized root mean square error (rRMSE): where x i ,x i and x i are the measured values of Cab and Car, the estimated values and the average measured values for number i, respectively; and i represents the sample number.

Remote Sensing Data Processing Results
The visualization and cross-sectional structure diagram of seven sampling trees is shown in Figure 5. The maximum crown width of single trees was mostly located at the middle and lower vertical position of canopies ( Figure 5(a 4 -g 4 )).  Figure 6 shows the box and whisker plots of the pigment content of Ginkgo leaves for different ages, heights, and their combinations. In terms of different ages (Figure 6a), except for Chla which showed a decreasing tendency with increase of age (decreased by 8.71%), Chlb (by 33.60%), Cab (by 2.1%) and Car (by 4.06%) were all increased. In terms of different height layers (Figure 6b), the vertical distribution patterns of different pigments were relatively consistent, with the highest content occurring in the upper vertical layer (7.41-12.20% higher than the middle layer; 3.18-8.00% higher than the lower layer), followed by the lower vertical layer, then the middle vertical layer. Particularly, in terms of the combinations between ages and heights (Figure 6c), the pigment content of the 13-year-old part samples generally increased with height, whereas the pigment content of the 22-year-old part samples generally decreased with height.

The Selection of Optimal LUT Entries for Estimating Cab and Car
The estimation accuracies of Cab and Car content using LUT-database-1 (50 LUTs) are shown in Figure 7. Compared with the RMSE function, the Cressie-Read function had better retrieval performance for Cab and Car content.  (Figure 7d). With the increase in the number entries of LUTs, the R 2 for Cab and Car generally showed a tendency to increase from 1000 to 4000, and then decrease or keep stable from 4000 to 16,000, with more concentrated scatter points for five random repetitions. Therefore, this paper used a LUT with 4000 entries, which had a relatively high R 2 and a relatively small RMSE, combined with the Cressie-Read cost function for subsequent pigments' inversion.

The Selection of Specific Optimal LUTs for Estimating Cab and Car
The horizontal and vertical distributions mapping of pigments within Ginkgo canopies from two strips (i.e., of 13-year-old and 22-year-old plantations) are shown in Figure 8. The results of the best-performing LUT for estimating Cab (R 2 = 0.60 for 13-year-old trees; R 2 = 0.49 for 22-year-old trees) and Car (R 2 = 0.46 for 13-year-old trees; R 2 = 0.36 for 22-year-old trees) are shown in Figure 8(a 1 -d 1 ). The horizontal and vertical distributions of pigments for two strips are shown in Figure 8(a 2 -d 2 ). The spatial distribution of Cab content and Car content within canopies exhibited similar patterns (Figure 8). In the horizontal direction, the pigment content of 13-year-old and 22-year-old ginkgo increased from inside the canopy to outside. In contrast, in the vertical direction, pigment content generally decreased along with the increase in the tree heights (Figure 8(a 2 -d 2 )). Figure 9 shows the stacked histogram of Cab and Car estimations for different heights from the plots in Figure 8. Cab and Car content of the 22-year-old plots were generally higher than those of the 13-year-old plots (Figure 9c), yielding an increase of 5.37-55.67% for Cab content and 3.37-42.88% for Car content at different height levels (Figure 9a,b). The stacked histograms for 13-year-old trees showed that high frequency occurred at the middle range (32-52 µg·cm −2 for Cab; 5.8-9.8 µg·cm −2 for Car; Figure 9a,b). The stacked histograms for the 22-year-old trees showed that high frequency occurred at the high range (34-64 µg·cm −2 for Cab; 7.8-11.8 µg·cm −2 for Car; Figure 9a,b).

Discussion
In this paper, Cab and Car content of ginkgo canopies were well estimated using the LUT-based-inversion approach and fused multi-scale remote sensing data (LiDAR and hyperspectral images). To our best knowledge, this is the first attempt to map the 3D distribution of pigment content within tree canopies by the synthetic use of the PROSAIL model and fused hyperspectral cloud points.

Performance of LUT-Based-Inversion in Estimating Pigment Content
For the estimation of pigment content, the LUT-based-inversion approach (R 2 = 0.36-0.60, Figure 8(a 1 -d 1 )) performed better than SI models (R 2 = 0.25-0.46, Supplementary Figure S3). In particular, the LUT-based-inversion approach can further improve inversion accuracy when applying the optimal entry size of 4000 for two ages with the Cressie-Read cost function (Figure 7). This can be explained by the fact that the PROSAIL model was developed based on radiative transfer principles, which help improve its universality and applicability in biochemical and biophysical estimations for multiple vegetation species and growth environments. Therefore, the LUT-based-inversion approach has a more robust estimation performance for pigments, which can be effectively used to determine the distribution of pigments within Gingko canopies of individual trees or forests.
Despite lower absorption features of Car than Cab, extracting the Car content by adopting the Cab-Car relationship might also lead to bias (overestimation for samples of Car < 6 µg/cm 2 ; underestimation for samples of Car > 10 µg/cm 2 ) in the worse performance of Car inversion (R 2 = 0.46-0.36, Figure 8(b 1 ,d 1 )) to some extent, as compared to the Cab estimation (R 2 = 0.60-0.49, Figure 8(a 1 ,c 1 )).

Potential for Mapping Pigments Content Based on Hyperspectral Fusion Point Cloud Coupling with PROSAIL Model
Numerous studies have proved that the combination of the PROSAIL model and hyperspectral data has the capability to estimate biochemical traits [43,65]. Different from previous studies which usually mapped the horizontal distribution of pigments only, this study was successful in coupling hyperspectral fusion point clouds with the PROSAIL model for mapping horizontal and vertical distribution of pigments. Specifically, based on the fusion remote sensing data, the LUT-based-inversion approach derived from the PROSAIL model achieved relatively preferable accuracies (R 2 = 0.36-0.60, Figure 8(a 1 -d 1 )) for mapping the 3D distribution of pigments [65].
Additionally, the analysis results of pigment distribution ( Figure 9) exhibited that the lower canopy had higher pigment content to some extent. This suggests that some uncertainties occurred when estimating pigments of total canopies using the upper canopy samples only. Thus, coupling the fused remote sensing data and PROSAIL model provides a promising approach for accurately mapping the 3D distribution of pigments.

Factors Affecting Pigment Distribution on Canopy Surfaces
In this study, the effects of age, light and canopy structure on pigment content and their distribution within forest canopies were analyzed.
In terms of age, with the growth of age, the tree growth curve presents an 'S' shaped growth process (i.e., slow-fast-slow) [9]. Therefore, the growth potential and photosynthetic capacity for ages 13-year-old and 22-year-old are different. The results showed that Cab and Car content of 22-year-old Ginkgo trees (middle stand age) was higher than that of 13-year-old (young-aged) (Figures 6 and 8). This may be attributed to the fact that 22-year-old Ginkgo trees are in the middle-aged growth stage, with stronger photosynthesis and more Cab accumulated in leaves.
In terms of light conditions, the content and distribution of Cab and Car in leaves are affected by light conditions, leading to variation in the strength of photosynthetic capacity of leaves [66]. In the horizontal direction, outer canopies received more direct solar radiation than the internal canopies, while the internal canopies are dominated by diffuse radiation [67]. This could be the main reason for the fact that the pigment content increases from inside the canopies to outside.
In terms of the canopy structure, the differences in LAI between various vertical layers may result in shadow effects to a different extent [68,69]. Additionally, the height and vertical position of leaves in the crown are important factors, which affect the pigment content and photosynthetic capacity [70]. The middle canopies have higher pigment content than upper and lower layers. This can be attributed to differences in the leaf expansion and senescence processes between the three height layers. Compared to middle layers, with the expansion of young leaves yielding higher photosynthetic capacity, young leaves might grow and wait for spreading out from the top canopies. In contrast, the leaves of the lower layers can be faster in the senescence process with decreased photosynthetic capacity compared to the other two height layers. Meanwhile, the maximum crown width mainly depends on the middle and lower canopy, resulting in a larger leaf area ratio and a wider light-receiving area in the middle and lower canopy, thus yielding higher content of Cab and other pigments synthesized in leaves. Moreover, the middle and lower leaves of the canopy may accumulate more Cab to obtain more light energy for alleviation of the impact of shadow effects [71,72].
Moreover, light and canopy structures have interaction influences with each other. However, it is complicated to quantify their respective effects on pigment estimation. We found that the estimated results from point-based reflectance (Figure 8) showed inconsistent patterns with the measured data of sampling trees ( Figure 6). This may be because the sampling trees and corresponding leaf samples are limited to some extent.
The change process for leaf biochemical pigments content is complex and continuous, comprehensively controlled and regulated by multiple factors to meet its growth needs. This paper indirectly proves that different height levels influence the pigment content of leaves. More importantly, the estimation of leaf biochemical traits using fused hyperspectral point clouds at different height levels can effectively improve 3D distribution information for pigments within canopies, which is consistent with most previous research results [45,69].

Future Research and Application
Previous studies proved that RTM coupled with machine learning algorithms can partly improve inversion accuracies (mainly modeled with a mean R 2 > 0.60 [73]). Thus, further studies may consider settings with machine learning algorithms (such as Gaussian processes regression) for estimation of pigments to improve accuracy. Furthermore, the PROSAIL model or other RTM can couple with the hyperspectral fusion point cloud to estimate the biochemical content of canopies in different regions, different tree species and different growth conditions in the future.

Conclusions
We attempted to perform multi-source remote sensing data inversion estimation in the PROSAIL model for biochemical pigments of Ginkgo, and fused hyperspectral imagery with LiDAR data to obtain the fused hyperspectral point clouds, which were put into the PROSAIL model to construct the exclusive LUT for Ginkgo. Following the LUT, we explored the 3D distribution of pigments within canopies and the regularities in the horizontal and vertical profiles of pigments with age and height. The fusion of multi-source remote sensing data provides a 3D visualization method for pigment estimation. The main results from this study are summarized as follows: 1.
The pigment content of Ginkgo trees increased from inward to outward in the horizontal direction, but decreased along with increase in height in the vertical direction; 3.
The hierarchical study in optimizing the inversion model once again emphasized its importance.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14030715/s1. They can be found at the Supplementary Materials. Table S1. Summary of flight parameters and UAV-LiDAR and UAV-hyperspectral sensor properties. Figure S1. Comparison of the maximum-minimum value, mean value and mean ± standard deviation of the measured reflectance spectra of the 42 part samples (six parts for seven individual sampling trees).  Figure S2. The pixel-based and point-based maximum-minimum, mean, and mean ± standard deviation reflectance curves for each of the two ages with six part samples.