Detection of Shoot Beetle Stress on Yunnan Pine Forest Using a Coupled LIBERTY2-INFORM Simulation

Abstract

Yunnan pine shoot beetles (PSB), Tomicus yunnanensis and Tomicus minor have spread through southwestern China in the last five years, leading to millions of hectares of forest being damaged. Thus, there is an urgent need to develop an effective approach for accurate early warning and damage assessment of PSB outbreaks. Remote sensing is one of the most efficient methods for this purpose. Despite many studies existing on the mountain pine beetle (MPB), very little work has been undertaken on assessing PSB stress using remote sensing. The objective of this paper was to develop a spectral linear mixing model aided by radiative transfer (RT) and a new Yellow Index (YI) to simulate the reflectance of heterogeneous canopies containing damaged needles and quantitatively inverse their PSB stress. The YI, the fraction of dead needles, is a physically-explicit stress indicator that represents the plot shoots damage ratio (plot SDR). The major steps of this methods include: (1) LIBERTY2 was developed to simulate the reflectance of damaged needles using YI to linearly mix the green needle spectra with the dead needle spectra; (2) LIBERTY2 was coupled with the INFORM model to scale the needle spectra to the canopy scale; and (3) a look-up table (LUT) was created against Sentinel 2 (S2) imagery and inversed leaf chlorophyll content (LCC), green leaf area index (LAI) and plot SDR. The results show that (1) LIBERTY2 effectively simulated the reflectance spectral values on infested needles (mean relative error (MRE) = 1.4–18%), and the YI can indicate the degrees of needles damage; (2) the coupled LIBERTY2-INFORM model is suitable to estimate LAI (R² = 0.73, RMSE = 0.17 m m⁻², NRMSE = 11.41% and the index of agreement (IOA) = 0.92) and LCC (R² = 0.49, RMSE = 56.24 mg m⁻², NRMSE = 25.22% and IOA = 0.72), and is better than the original LIBERTY model (LAI: R² = 0.38, RMSE = 0.43 m m⁻², NRMSE = 28.85% and IOA = 0.68; LCC: R² = 0.34, RMSE = 76.44 mg m⁻², NRMSE = 34.23% and IOA = 0.57); and (3) the inversed YI is positively correlated with the measured plot SDR (R² = 0.40, RMSE = 0.15). We conclude that the LIBERTY2 model improved the reflectance simulation accuracy of both the needles and canopies, making it suitable for assessing PSB stress. The YI has the potential to assess PSB damage.

1. Introduction

The forest ecosystem plays an important role in the global carbon cycle [1,2]. However, forest pests can destroy the stability of forest ecosystems and increase their vulnerability to natural disturbance [3,4]. There are two major types of pests, each attacking different organs of trees: defoliators eat leaves or needles and borers bore into the shoots, bark, or wood. These pests may cause a temporary physiological decline in trees (i.e., reduced moisture) or permanent damage to biochemical and biophysical parameters like leaf chlorophyll content (LCC) and green leaf area index (LAI). These changes result in canopy reflectance changes, making it possible to use Earth observation (EO) data to detect forest pests [5,6]. Remote sensing has become an effective technology for detecting forest damage (i.e., diseases and insects) because of the richness of high temporal, spatial and spectral resolution images [7,8,9,10,11,12,13].

Defoliator damage is relatively easier to detect due to the visible symptoms of forest cover and leaf discoloration. For example, using a time series approach, EO data can effectively characterize defoliator impact [14,15]. However, borer (i.e., beetle) attacks are often hidden in the trunk or shoots. The discoloration symptoms have a lag of several weeks, causing difficulties in detecting them. Despite various remote sensing algorithms having been proposed to assess forest damage, prediction of beetle attacks is still a challenge [16,17,18]. Most empirically-based algorithms link variations in spectral reflectance to levels of pest damage or stages of decline [19]. Chen and Meentemeyer [19] categorized them into four groups: (1) thresholding involves defining one or multiple thresholds to isolate the infested trees from the entire forest using spectral indices (i.e., RGI [20], NDVI [21], NDMI [22,23], DI [24], EWDI [25]); (2) classification involves classifying forest damage based on the distinct symptoms of damaged trees; (3) statistical regression involves correlating image-derived explanatory variables (i.e., spectral bands, spectral ratio indices, and spectral vegetation indices) to field-measured damage indicators (i.e., LAI, LCC, canopy chlorophyll content, defoliation intensity) to assess or predict forest vulnerability to pest attacks; and (4) change detection involves employing multi-temporal EO data to identify the spatial patterns of forest disturbance through time. These empirical methods may successfully detect forest damage under certain conditions; however, they do not explain the mechanism of the spectral response to vegetation stress [26].

Theoretically, the canopy reflectance change as a result of damage is mainly caused by the changes in leaf pigment content and canopy geometrical structure parameters in response to insect attacks [18,26,27]. To increase the detection accuracy of beetle damage, the mechanistic link between vegetation parameters and spectral features has to be clarified [28]. Instead of using empirical simple equations or thresholds, canopy radiative transfer models (RTMs) quantitatively depict the whole process of light transfer within and above the vegetation canopy. In RTMs, top-of canopy reflectance is a function of canopy geometrical structure parameters (e.g., LAI), sun and viewing geometry, and leaf and background soil optical properties. There are four types of canopy RTMs [29,30]: (1) RT models (i.e., SAIL [31]) including plate models, N-flux models, and stochastic models, which are generally applied in continuous and homogeneous vegetation canopies [29]; (2) geometrical optical (GO) models (i.e., 5-SCALE [32]), applied to canopies with distinct architecture particularly sparse forest canopies [33,34]; (3) hybrid GORT models (e.g., Invertible Forest Reflectance Model (INFORM) [35,36]), which combine the advantages of GO models with the hotspot effect and RT models with multiple scattering processes [37]; and (4) computer simulation models (i.e., DART [38]), which use computer graphics to generate the realistic 3-dimensional (3D) structure of the vegetation scene and accurately render the scene [39,40,41]. For these computer simulation models, Monte Carlo ray tracing and radiosity are the most widely used theories. Leaf reflectance and transmittance can be user-defined or simulated by a leaf RTM, which is a function of leaf structure and biochemical parameters (e.g., LCC). Therefore, the coupled model (a canopy RTM and a leaf RTM) can be used to explore the canopy reflectance response of the variation of leaf biochemical and partial canopy biophysical parameters (e.g., LCC and LAI) [42,43]. However, it is important to choose the suitable RTM for a specific forest canopy and application scenery.

For Yunnan pine, INFORM was chosen as the canopy RTM and the LIBERTY [44] (Leaf Incorporating Biochemistry Exhibiting Reflectance and Transmittance Yields) model as the leaf RTM. Because the leaf of Yunnan pine has a long needle shape, the LIBERTY model was better than the PROSPECT [45] model. Due to the dry climate and Karst rocky desertification, Yunnan pine trees grow poorly in this area, with large canopy gaps, low canopy height, and a sparse distribution of forest. Therefore, GORT and computer models were best suited. Considering the low time cost, GORT models were chosen. GORT has two types of input parameters: (1) GO parameters, which include sunlit or shaded view fractions of soil or crowns; and (2) RT parameters, which include the brightness of sunlit or shaded soil or crowns. The Yunnan pine crown diameter, density, and tree height were used to estimate the GO parameters. LIBERTY coupled with a canopy RTM was used to simulate the optical properties of sunlit or shaded components (soil and crown). INFORM was a GORT model combining the forest light interaction model (FLIM [46]), SAILH [31], and a suitable leaf model (e.g., PROSPECT or LIBERTY). The INFORM model was suitable to simulate sparse forest because, firstly, it accounts for the effects of shadowing, canopy geometric structure, crown transmittance, the hotspot and the clumping of leaves in crowns on the forest reflectance. Previous studies [36,47] have indicated that the INFORM model achieves a good performance in medium-dense and dense managed conifer stands with low a variation in tree height and open canopy. Secondly, the FLIM model considered the canopy geometrical structure and the SAILH model combined with PROSPECT or LIBERTY implemented a full multiple scattering scheme to simulate discontinuous canopy optical properties. Due to its good documentation and mature software framework for inversion, as well as the embedded LIBERTY [44] model, INFORM was the best choice for our study.

Although the LIBERTY-INFORM model was suitable for Yunnan pine, it still cannot describe the non-uniform distribution of needle chlorophyll, which produces the significant symptoms of heterogeneous discoloring, and red dieback damage of shoots. Pine shoot beetles (PSB) bore into the shoots and block water and nutrient delivery to the needles. This leads to the color of needles above the wormhole of the shoot changing from green to red during the infestation stages, while the color of the needles below the wormhole of the shoot remain green until the whole shoot undergoes a significant physiological decline. PSB thereby cause heterogeneous needles with patchy chlorosis which strongly influences the whole leaf and canopy reflectance due to the non-uniform distribution of chlorophyll [48]. The non-uniform distribution of chlorophyll may impact RTM modeling leaf/canopy reflectance accuracy and the ability to deduce LCC and LAI. However, previous studies have ignored the effect of the non-uniform distribution of chlorophyll on canopy model reflectance for the estimation of biochemical and biophysical parameters in vegetation stress [49,50,51]. In addition, pigment composition such as the ratio between chlorophylls a and b and the relative amount of carotenoids can change in response to certain forms of stress. Many studies [27,52] have indicated that the chlorophyll a/b ratio becomes higher with deeper damage and that the ratio of total chlorophyll (a + b) to total carotenoids in highly damaged leaves was significantly lower than in healthy or slightly damaged leaves. Furthermore, stress pigments such as carotenoids and anthocyanins are strongly responsible for controlling the color of the leaves in response to change in environmental conditions. Leaf optical properties are directly impacted by the composition and distribution of pigments. Therefore, the LIBERTY model should be extended to include the non-uniform distribution of chlorophyll and the changed pigment components for modeling heterogeneous needles. Efficient inversion approaches like look-up-table (LUT [53]) should be tested, adapted, and used to estimate these vegetation parameters from measured spectral data.

In this paper, instead of attempting to quantify the concentration of pigments such as carotenoids and anthocyanins, we aimed to develop an extended LIBERTY-INFORM simulation method (named LIBERTY2-INFORM) to deal with heterogeneous needles and sparse spatial distribution so as to improve the detection accuracy of PSB stress in Yunnan pine forests. Section 2 gives the material and methods. Then, the results and discussion are presented in Section 3 and Section 4, respectively. The conclusions are in the final section.

2. Material and Methods

2.1. Study Site

The study site (25°14′–25°29′N, 100°48′–101°3′E), is located in the Tianfeng Mountains, Yunnan Province, in the southwest of China (Figure 1). The forest is dominated by air-sowing stands of Yunnan pine (Pinus yunnanensis Franch) with an area of about 1000 ha. Since 2010, PSB (Tomicus yunnanensis and Tomicus minor), which bore into both shoots and trunks, have been damaging the Yunnan pines, which causes serious forest dieback (approximately 50% of the forest area to date).

The field campaign was carried out from October to November in 2016. In order to cover different degrees of the damaged forest, 34 plots (30 m by 30 m), representing different levels of damage (11 healthy plots and 23 stressed plots), were set in the study area. The central GPS coordination of each plot was measured using a differential GPS device (GeoExplorer 6000, Trimble Inc., Sunnyvale, CA, USA) with an accuracy of about ±0.75 m. To quantify the tree damage, the shoots damage ratio (SDR) has been defined and measured [54]. The SDR (shown in

Table 1. The parameter statistics for the forest plot variables (n = 34).

Type of Plots	Mean		SD		Max		Min		Range
	HP	SP	HP	SP	HP	SP	HP	SP	HP	SP
CD [m]	2.5	2.3	0.27	0.22	5.5	5.3	0.5	0.5	5	4.8
H [m]	4.68	5.07	1.5	1	11	12	1	1	10	11
CO [%]	28	32	8.9	6.8	42	44	15	20	27	24
SD [ha−1]	990	1342	404	556	2022	2644	500	556	1522	2088
LAI [m2 m−2]	0.78	0.83	0.37	0.29	1.86	1.61	0.43	0.37	1.43	1.24
LCCplot [mg m−2]	422	342	28.21	61.14	428	389	372	205	56	184
SDR	0.03	0.37	0.06	0.31	0.23	1	0	0	0.23	1
Plot SDR	0.00	0.34	0.03	0.13	0.02	0.65	0	0.1	0.02	0.55

Note: SD = Standard deviation, Max = Maximum, Min = Minimum, HP = Healthy Plots, and SP = Stressed Plots.

) is the ratio of damaged shoots to total shoots of a single tree. For each tree, all damaged and healthy shoots were counted to estimate the SDR of a single tree. At each plot, we measured all trees with a diameter at breast height (DBH) greater than 4 cm. By averaging the SDR of the total of all eligible trees in a plot, the plot SDR (shown in Table 1) was calculated by:

$$plotSDR=\frac{1}{N_{\mathrm{trees}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{trees}}}{SDR}_i}$$

(1)

where N_trees is the number of trees in a plot with a DBH greater than 4 cm.

2.2. Measurements of Forest Structural Parameters

The forest structural variables were measured in each plot, including green LAI, canopy closure (CO), stem density (SD), stand average height (H), and stand average crown diameter (CD). CO was sampled in two cross-line of each plot. SD was obtained by counting the number of trees in a plot. H and CD were calculated from the average of all trees (DBH > 4 cm) within each plot. Plot green LAI was averaged using 25 evenly-spaced measurements. Each LAI measurement was conducted by a Li-Cor LAI-2200 plant canopy analyzer with a 270-degree view restrictor. Clumping is simultaneously estimated by the LAI-2200 [55,56]. In the winter season, there is a low density of understory. Summary statistics for the forest plot variables are given in Table 1.

2.3. Measurements of Needle Reflectance and Chlorophyll Content

The feeding behaviors of T. yunnanensis and T. minor are divided into two stages: (1) boring into the shoot; and (2) moving to the trunk. During the first stage, damaged shoots block the water and nutrient transfer to the needles, resulting in localized regions of chlorosis across the needle surface. In one needle, the color may vary from green to yellow. During the second stage, the “red-attack” stage, the needle color turns entirely to red. Therefore, it is interesting to measure the needle reflectance with different color combinations.

According to the common color changing sequence (green, light green, yellow, red) of pine needles, four degrees of needle damage were defined, including healthy (green), slight damage (light green), moderate damage (green-yellow, yellow), and severe damage (yellow-red, red). For each level, 50 samples were collected, sealed, and brought to an indoor laboratory. Two Avantes spectroradiometers, AvaSpec-EDU-VIS (350–850 nm), and AvaSpec-NIR1.7 (900–1700 nm) were combined to measure the needle reflectance indoors from visible (VIS) to near-infrared (NIR) bands. The spectrum resolutions were 1 nm in VIS and 3 nm in NIR, respectively. Reflectance was normalized by the sample spectrum from a white reference panel (98% diffuse reflectance). The process of the spectral measurement was conducted in a dark room. In each sample with the same degree of needle damage, about 30 to 40 number needles were tightly arranged in a grooving device with a clip because one needle is too slim to fully fill the field of view (FOV) of the optical fiber probe. For both the reflectance and transmittance measurements, the AVANTES Halogen light source (AvaLight–DHc) was used, which has a wavelength range from 400 nm to 2500 nm. Reflectance fiber with a wavelength ranging from 200 nm to 2500 nm was used, which was contacted with the needles. Transmission Dip Probes (TDP, 2 m length) were mounted on the light source and used to measure the needle transmittance. To avoid external light interference, the optical fiber probe with a FOV of 25° was installed in a reflection probe holder (RPH-1, https://www.avantes.com/products/fiber-optics/item/260-reflection-probe-holders). The fiber probe holder was attached to the sample grooving device to measure the needle spectra. Each sample (from the tip to base of the needles) was measured five times and averaged to minimize instrument noise.

Correspondingly, the LCC of the needles were measured using the Chlorophyll Content Meter (CCM-300, Opti-Science, Inc., Hudson, NH, USA) based on the fluorescence ratio of red emissions measured at 700 nm and 735 nm. The chlorophyll (Chl) content was calculated using the equation [Chl] = 634F735/F700 + 391 [57]. There was a linear response to Chl content ranging from 41–675 mg m⁻², with a determination coefficient of R² greater than 0.95. The CCM-300 was capable of conducting nondestructive measurement of the chlorophyll of tiny leaves (such as needles) with a leaf clip holder. The CCM-300 was calibrated (R² = 0.94) using 28 independent samples measured by a standard UV-2600 Spectrophotometer and using 95% anhydrous ethanol as dissolving reagent [58].

For each sample, three sections of the tip, center, and base of the needles were measured to calculate the average LCC. The results of the LCC are shown in

Table 2. The statistics of LCC (unit: mg m⁻²) for different degrees of damaged needles and shoots. (50 samples for each class of needles or shoots).

	Statistics	Healthy	Slight	Moderate	Severe
Needles	Max	495	406	303	60
	Min	385	216	29	0
	Mean	433	328	138	19
Shoots	Max	485	398	262	98
	Min	400	256	102	0
	Mean	428	334	174	40

.

Similar to the needle scale, the infested shoots are also divided into healthy, slight, moderate and severe levels of damage (Figure 2). The severity of damage to individual shoots was visually estimated: (a) healthy, for shoots with green needles and no insect attack; (b) light, for shoots with a boring attack and consisting of light or yellowish needles; (c) moderate, with half yellow needles; and (d) severe, for shoots with almost completely red needles. The mean LCC of the four classes of needles was used to determine the LCC of a single shoot. For each class of shoots, 50 shoot samples with varied needle density and length were chosen to calculate the mean LCC of the shoots, which was then defined as the shoot LCC (LCC_shoot). The statistics of LCC_shoot are shown in Table 2.

It is difficult to determine the average LCC for a stressed tree or plot. Traditionally, the total plot LCC is a multiplied value of a mean LCC and LAI [59,60]. This method is easy and quick for estimating the total plot LCC when the canopy leaves are relatively homogenous. However, our key question is around how to determine the mean LCC for a heterogeneous canopy. For Yunnan pines, the trees may have different degrees and numbers of needles or shoots. The heterogeneous crowns lead to difficulty in determining the mean LCC through random sampling of needles or shoots. In this paper, the LCC of a single tree (LCC_tree) is calculated with a weighted average method:

$$LC{C}_{\mathrm{tree}}=\frac{LC{C}_{\mathrm{shoot},\mathrm{healthy}}{n}_{\mathrm{shoot},\mathrm{healty}}+LC{C}_{\mathrm{shoot},\mathrm{slight}}{n}_{\mathrm{shoot},\mathrm{slight}}+LC{C}_{\mathrm{shoot},\mathrm{moderate}}{n}_{\mathrm{shoot},\mathrm{moderate}}+LC{C}_{\mathrm{shoot},\mathrm{severe}}{n}_{\mathrm{shoot},\mathrm{severe}}}{N_{\mathrm{shoots}}}$$

(2)

where LCC_{shoot,healthy}, LCC_shoot,slight, LCC_{shoot,moderate}, and LCC_shoot,severe are the mean chlorophyll content of the healthy, slight, moderate, and severely damaged shoots, respectively. N_shoots is the total number of shoots of a single tree.

Then, the average LCC of a plot with N_trees trees, the LCC_plot (see Table 1), is as follows:

$$LC{C}_{\mathrm{plot}}=\frac{1}{N_{\mathrm{trees}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{trees}}}LC{C}_{\mathrm{tree},\mathrm{i}}}$$

(3)

2.4. Sentinel 2A Images and Processing

The Sentinel-2 (S2) satellite features a Multi-Spectral Instrument (MSI) with a 290 km field of view (FOV). This MSI provides 13 spectral bands, ranging from the visible and near-infrared (VNIR) to the short-wave infrared (SWIR), with four bands (B2–B4, B8) at 10 m spatial resolution, six bands (B5–B7, B8a, B11, B12) at 20 m spatial resolution, and three bands (B1, B9, B10) at 60 m spatial resolution [61]. The high spatial resolution, wide FOV, and broad spectral coverage of S2 are designed to map land use and land cover (LULC), as well as change detection, geophysical variables (i.e., LCC, LAI) mapping, and risk, such as forest insect detection. S2 imagery has a high quality of geometry and radiation [62]. The Level 1C product is geo-coded with a sub-pixel registration accuracy for users and provides the top of the atmosphere (TOA) reflectance [62]. Images collected on 31 December 2016 (Relative Orbit number: R104) were used in this paper, taken at a time that was almost cloud-free (cloud cover less 1%). The atmosphere correction was conducted using the SEN2COR (http://step.esa.int/main/third-party-plu-gins-2/sen2cor/) tools. We used the bands from B2 to B8a, which are highly sensitive to chlorophyll and LAI variations. The bands at a spatial resolution of 10 m were aggregated to 20 m, resulting in eight bands in VNIR at 20 m.

2.5. The Needle-Reflectance Model LIBERTY

The LIBERTY model was developed to describe the spectral properties of conifer needles [44]. The model employs nine biophysical/biochemical parameters (

Table 3. The ranges and distribution of input parameters to generate LUT.

Type of Variable	Input Parameters	Designation	Unit	Range and Step	Distribution
L	Cell diameter	d	m−6	40–60, 5	Uniform
L	Intercellular airspace	xu	arbitrary	0.03–0.06, 0.005	Uniform
L	Leaf thickness	t	arbitrary	1–10, 1	Uniform
L	Baseline absorption	b	arbitrary	0.0005	-
L	Albino absorption	a	arbitrary	2–4, 0.5	Uniform
L	Leaf Chlorophyll content	LCC	mg m−2	200–450, 5	Uniform
L	Water content	CW	g m−2	100	-
L	Lignin and cellulose content	CL	g m−2	40	-
L	Nitrogen content	CP	g m−2	1	-
I	Yellow index	YI	arbitrary	0–0.5, 0.05	Uniform
C	Single tree of leaf area index	LAIs	m2 m−2	0.1–4.5, 0.5	Uniform
C	LAI of understory	LAIU	m2 m−2	0–1, 0.2	Uniform
C	Stem density	SD	ha−1	500–2500, 50	Uniform
C	Average leaf angle	ALA	deg	30–70, 5	Uniform
C	Tree height	H	m	1–12, 1	Uniform
C	Crown diameter	CD	m	0.5–5.5, 0.5	Uniform
C	Hot spot parameter	hot	m m−1	0.02	-
E	Sun zenith angle	${\theta }_s$	deg	52.50	-
E	Observation zenith angle	${\theta }_o$	deg	7	-
E	Relative Azimuth	phi	deg	0	-
E	Fraction of diffuse radiation	skly	fraction	0.1	-

Type of variable: E = External, C = Canopy parameters, L = Leaf parameters of LIBERTY, and I = Introduced.

) to simulate needle reflectance and transmittance in the VNIR and SWIR wavelengths (400–2500 nm) at 5 nm intervals. The model linearly sums in vivo specific absorption coefficients, scaled by their respective concentrations, to calculate a global absorption coefficient [63]. The global absorption coefficient, in conjunction with the structural parameters of average cell diameter (d), intercellular air gap (xu), and needle thickness (t), were used to calculate needle reflectance and transmittance.

2.5.1. Calibration of Pigment Absorption Coefficients

The LIBERTY model in vivo specific absorption coefficients were measured from jack pine needles during the Boreal Ecosystem Atmosphere Study (BOREAS) campaign [44]. The leaf model (e.g., PROSPECT/LIBERTY) physically describes the relationship between biochemical/structure parameters and spectral reflectance; thus, it is expected to be independent of species and site conditions. However, the specific absorption coefficients of pigment in different species are not completely consistent. Many studies indicated that better accuracy of leaf/needle optical properties and biochemical parameters were estimated after being validated and calibrated with real measurement of leaf/needle reflectance [44,45,63,64,65,66,67,68]. Yunnan pine is a different species to jack pine, thus, the LIBERTY model needs to be calibrated before it can accurately simulate vegetation reflectance spectra at the local scale and inversed to estimate leaf biochemical parameters [64,65]. In order to accurately model needle reflectance, calibration of pigment absorption coefficients should be considered. The Newton–Raphson iterative method was applied to calibrate the LIBERTY model [63] using a randomly selected 40% of healthy needles from field measurements for model calibration. The remainder were used for model validation. The validation accuracy was assessed by the mean relative error (MRE) between the modeled and measured needle reflectance.

2.5.2. Simulation of Heterogeneous Leaf Reflectance

Insect damage can cause heterogeneous leaves with patchy chlorosis which strongly influences the whole leaf reflectance because of the non-uniform distribution of chlorophyll [48]. According to the BARTON study [48], a heterogeneous leaf area was split into two regions with distinct chlorophyll concentrations: (1) the damaged region with the fraction of leaf (f) and chlorophyll concentration C_α; and (2) the healthy region (1 − f) with chlorophyll concentration C_β. The mean leaf chlorophyll concentration C_mean is

$$C_{\mathrm{mean}}=f{C}_{\alpha }+(1-f)C_{\beta }$$

(4)

Assuming that C_α = 0 for the chlorophyll content of the dead leaf area, and f is the dead leaf area fraction, Equation (4) is simplified as

$$C_{\mathrm{mean}}=(1-f)C_{\beta }$$

(5)

Because f represents the degree of needle damage, we use it as a damage indicator called the “Yellow Index” (YI) [69]. Therefore, Equation (5) can be converted into a new form (C_mean = (1 − YI)C_β), in which C_β refers to the chlorophyll content of the healthy region for a heterogeneous needle.

Assuming that two regions of a heterogeneous leaf are scattered or absorbed without interaction, the reflectance R_λ,Cmean is linearly mixed [48] between the dead region reflectance R_λ,Cα=0 and the healthy region reflectance R_λ,Cβ:

$$R_{\lambda ,C_{\mathrm{mean}}}=\mathrm{YI}{R}_{\lambda ,C_{\alpha }=0}+(1-\mathrm{YI})R_{\lambda ,C_{\beta }}$$

(6)

Similarly, the mixed transmittance of a heterogeneous leaf T_λ,Cmean is as follows:

$$T_{\lambda ,C_{\mathrm{mean}}}=\mathrm{YI}{T}_{\lambda ,C_{\alpha }=0}+(1-\mathrm{YI})T_{\lambda ,C_{\beta }}$$

(7)

where T_λ,Cα=0 and T_λ,Cβ are the transmittance of the dead region and healthy region, respectively. The healthy needle spectra (R_λ,Cβ, T_λ,Cβ) can be effectively simulated by the LIBERTY model following calibration with healthy needles. However, under stress conditions, the LIBERTY model performs bad and should be modified to achieve better accuracy [69]. The reason for this is that the LIBERTY model only considers a single pigment absorption coefficient for chlorophyll and carotenoid instead of separating them into individual pigments [44,70]. In fact, chlorosis results from the interaction of multiple pigments of the leaf. Pigments such as carotenoid, brown, and anthocyanins, in particular, cause a yellow or red leaf [71]. To overcome this limitation of LIBERTY, YI was used as a corrected fraction to describe the effect of the mixture of multiple pigments correlated with chlorosis on needle reflectance. In the process of the needle reflectance simulation, YI was independent of the radiative transfer modeling equation of LIBERTY. The measured spectra of red needles (LCC = 0 mg m⁻²) of Yunnan pine was used as the prior-knowledge of the dead needle spectra (R_λ,Cα=0, T_λ,Cα=0); it can be replaced by a generally dead or senescent leaf spectra.

For convenience, the original, calibrated (including the calibrated pigment absorption coefficient), and modified (includes the calibrated pigment absorption coefficient and YI) LIBERTY models are named LIBERTY0, LIBERTY1, and LIBERTY2, respectively. The performance of the infested needle reflectance simulation was evaluated using the mean relative error (MRE) based on measured and simulated needle reflectance.

As described above, there are four levels of LCC: needle, shoot, single tree, and plot. SDR is this concept at the tree level. At the plot scale, the mean SDR of all trees is used to represent the degree of damage. Therefore, the plot LCC is a combined value of three scales (tree, shoot, and needle):

$$LC{C}_{\mathrm{plot}}=\frac{1}{N_{\mathrm{trees}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{trees}}}\left(LC{C}_{\mathrm{shoot},\mathrm{healthy}}\left(1-SD{R}_i\right)+LC{C}_{\mathrm{shoot},\mathrm{damaged}}SD{R}_i\right)}$$

(8)

$$LC{C}_{\mathrm{plot}}=\frac{1}{N_{\mathrm{trees}}{N}_{\mathrm{shoots}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{trees}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{shoots}}}\left(LC{C}_{\mathrm{needle},\mathrm{healthy}}\left({1-\mathrm{YI}}_{\mathrm{j}}\right)+LC{C}_{\mathrm{needle},\mathrm{damaged}}{\mathrm{YI}}_{\mathrm{j}}\right)}}$$

(9)

where LCC_{shoot,healthy}, LCC_{shoot,damaged}, LCC_{needle,healthy}, and LCC_{needle,damaged} are the mean LCC of the healthy and damaged shoots or needles.

Assuming that all trees in a plot have the same degree of damage and needles or shoots are randomly distributed in the tree crowns with either healthy or dead properties, Equations (8)–(9) can be derived as

$$LC{C}_{\mathrm{plot}}=LC{C}_{\mathrm{shoot},\mathrm{healthy}}\left(1-\mathrm{plot}SDR\right)$$

(10)

$$LC{C}_{\mathrm{plot}}=LC{C}_{\mathrm{needle},\mathrm{healthy}}\left(1-\mathrm{YI}\right)$$

(11)

Based on the assumption above, LCC_{shoot,healthy} is equal to LCC_{needle,healthy}. Although YI is defined at the needle scale, by assuming that needles are randomly distributed, YI is numerically equal to the plot SDR.

2.6. INFORM Model

The INFORM model [36], coupled with LIBERTY0, LIBERTY1, and LIBERTY2, is used to simulate the bidirectional reflectance (BRF) of forest stands between 400 and 2500 nm. The main input parameters of the INFORM model are shown in Table 3.

2.7. Look-Up-Table (LUT) Inversion

LUT is a simple but effective strategy for remote sensing inversion, which can help to alleviate the drawback of the ill-posed problem [72]. The LUT inversion can be divided into two steps: (1) The generation of LUT by RTM; and (2) the selection of the solution corresponding to the given measurement. LUT consists of input parameter combinations and the associated output simulation of canopy reflectance. For the generation of LUT, the range and distribution of the LIBERTY and INFORM model parameters should be defined. Table 3 shows the ranges (minimum, maximum, and a step) of the leaf and canopy variables. The minimum and maximum of single tree LAIs (LAIs = LAI/CO), understory LAI (LAI_u), LCC, H, CD, SD, thee hotspot parameter (hot), and the diffuse light fraction (skly) are defined based on prior knowledge from field measurements and other similar studies [47]. According to the sensitivity analysis of LIBERTY [44,63,64], water, lignin/cellulose, and protein contribute mostly in NIR with wavelengths greater than 1000 nm [44,63]. Therefore, the water content (C_W), lignin and cellulose content (C_L), and nitrogen content (C_P) were held constant. Sun and viewing conditions correspond to the situation of the satellite overpass. A LUT length of 100,000 cases of parameter combinations was chosen using uniform random sampling in accordance with previous research [73,74]. This LUT size proved to be feasible and a larger LUT size may not influence the estimation accuracy according to similar research [55,56,72]. Subsequently, the coupled models LIBERTY0-INFORM, LIBERTY1-INFORM, and LIBERTY2-INFORM were applied to simulate canopy reflectance (ρ_canopy) using the 100,000 cases of parameter combinations. The simulated canopy reflectance with 5 nm intervals was re-sampled to the eight bands of S2A using S2 MSI spectral response functions (SRF) (https://earth.esa.int/web/sentinel/document-library/content/-/article/sentinel-2a-spectral-responses) with the following form [75]:

$${\rho }_{\mathrm{image}}\left(\lambda \right)=\frac{{\displaystyle {\int }_{{\lambda }_{\max }}^{{\lambda }_{\max }}{\rho }_{\mathrm{canopy}}\left({\lambda }_o\right)f\left({\lambda }_o\right)}d{\lambda }_o}{{\displaystyle {\int }_{{\lambda }_{\max }}^{{\lambda }_{\max }}f\left({\lambda }_o\right)}d{\lambda }_o}$$

(12)

where ρ_image (λ) refers to the broadband reflectance of S2A for given band λ, f(λ_o) refers to the SRF of the corresponding wavelength, and λ_max and λ_min refer to the upper and lower limit of each band interval, respectively. Finally, the input parameter combinations and associated output simulation of canopy reflectance were stored in LUT.

In the S2A imagery, we used the average reflectance of 3 × 3 pixels around the plot location to reduce the uncertainty of location error due to the inconsistency between plot dimensions (30 m) and the S2A pixel resolution (20 m). The relative error (RE) (see

Table 4. The statistic of RE (%) between 3 × 3 pixels and the plot location pixel for 34 plots.

	Mean	Max	Min	SD
RE	4.49	13.40	0.08	3.06

) between the average reflectance of 3 × 3 pixel reflectance and plot location pixel reflectance was used to explore pixel variation. Table 4 shows that the difference between the plot location pixel and 3 × 3 pixels was small, at 4.49% of the mean RE. The average reflectance of 3 × 3 pixels can represent the measured reflectance of each plot.

Table 5. The twelve cost functions used in LUT-inversion.

Cost Function	Algorithm
Shannon (1948)	$D[P,Q]=-{\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}(\frac{p({\lambda }_l)+q({\lambda }_l)}{2})}\log (\frac{p({\lambda }_l)+q({\lambda }_l)}{2})+\frac{1}{2}({\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}p({\lambda }_l)\log (}p({\lambda }_l))+{\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}q({\lambda }_l)\log (}q({\lambda }_l)))$
L-divergence lin	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}p(}{\lambda }_l)\ln (p({\lambda }_l)+q({\lambda }_l))\ln (q({\lambda }_l)-(p({\lambda }_l)+q({\lambda }_l)))\ln (\frac{p({\lambda }_l)+q({\lambda }_l)}{2})$
K(x) = log(x) + 1/x	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}\left(\log \left(\frac{q({\lambda }_l)}{p({\lambda }_l)}\right)+\frac{q({\lambda }_j)}{p({\lambda }_j)}\right)}-1$
K(x) = -log(x) + x	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}\left(-\log \left(\frac{q({\lambda }_l)}{p({\lambda }_l)}\right)+\frac{q({\lambda }_j)}{p({\lambda }_j)}\right)}-1$
K(x) = x(log(x)) − x	$D[P,Q]=1+{\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}\frac{q({\lambda }_l)}{p({\lambda }_l)}\left(\log \left(\frac{q({\lambda }_l)}{p({\lambda }_l)}\right)-1\right)}$
Jeffreys-Kullback-Leibler	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}(}p({\lambda }_l)-q({\lambda }_l))(ln(p({\lambda }_l)-q({\lambda }_l)))$
Exponential	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}q({\lambda }_l)(\exp (}-\frac{p({\lambda }_l)-q({\lambda }_l)}{q({\lambda }_l)})-1)$
RMSE	$D[P,Q]=\sqrt{\frac{{\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}(}p({\lambda }_l)-q({\lambda }_l){)}^2}{n}}$
Normal distribution-LSE	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}(}p({\lambda }_l)-q({\lambda }_l){)}^2$
Least absolute error	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}\left|p({\lambda }_l)-q({\lambda }_l)\right|}$
Neyman chi-square	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}\frac{{(p({\lambda }_l)-q({\lambda }_l))}^2}{q({\lambda }_l)}}$
Generalised Hellinger	$D[P,Q]={\displaystyle \sum _{{\lambda }_1=1}^{{\lambda }_n}{(\sqrt{p}-\sqrt{q})}^2}$

shows the 12 different potentially promising cost functions (CFs) for LUT inversion from three broad families (information measures, M-estimates, and the minimum contrast method) [76]. As demonstrated by Rivera [55], these CFs should be competed with one another to find the best inversion solution. The competition is based on D [P, Q], the distance between two functions, where $P=(p({\lambda }_1),\dots ,p({\lambda }_n))$ is the reflectance of satellite data; $Q=(q({\lambda }_1),\dots ,q({\lambda }_n))$ is the LUT simulated reflectance; and λ₁, …, λ_n represent the n spectral bands.

Other than the optimization of cost functions, two regularization options are also applied in LUT-based inversion strategies to further optimize the robustness of the estimates [55,56]. First, a series of Gaussian white noise, from 0–20% with a step 1%, was added to LUT simulate canopy reflectance to account for uncertainties attached to the models and measurements [55,56,73,77]. Second, multiple best solutions were used instead of a single best solution to calculate the estimated values [55,56,73,78]. Several studies have demonstrated that the single best parameter combination corresponding to the smallest distance calculated by a cost function (e.g., RMSE) does not necessarily lead to the best accuracies [43,78]. Different optimization numbers and the mutual effect on these regularization options have been systematically assessed in previous studies [55,56].

The inversed LCC and LAI values were validated using field measurements. The accuracy was assessed by R², the index of agreement (IOA) [79] (IOA is a standardized measure of the degree of model prediction error and varies between 0 and 1), root mean square error (RMSE) and the normalized RMSE (NRMSE) [72]. The inversion performance of INFORM using three versions of LIBERTY models was evaluated by calculating the deviations between the retrieved values and measured values through RMSE and NRMSE. The NRMSE was used to determine the best CFs and optimal regularization options.

3. Results

3.1. Pigment Absorption Coefficient Calibration and Needle Reflectance Simulation

Figure 3 shows the pigment absorption coefficient of the original and calibrated LIBERTY models using the measured needle spectral. It is clear that the values of the calibrated pigment absorption coefficient are lower than the original values in most of the wavelengths. By increasing the YI from 0 to 1 with a fixed chlorophyll content of the healthy region (C_β) of 400 mg m⁻², the corresponding heterogeneous needle reflectance (Figure 4a) was simulated by LIBERTY2. The first derivative is also presented in Figure 4b. From Figure 4, we can see the effect of the chlorotic area fraction on needle reflectance:

(a)

The needle reflectance decreases significantly in green and NIR bands but increases in the red band in response to the increase of the YI (Figure 4a).

(b)

In the first derivative of the heterogeneity of needle reflectance, the red shift in the green peak and blue shift in the red valley are visible due to the damage of pest stress (Figure 4b).

(c)

A slight change in the red edge position (710 nm) is found (Figure 4b).

Figure 5 shows the reflectance simulations using LIBERTY0, LIBERTY1, and LIBERTY2 for different severities of damaged needles. The LIBERTY1 model performs better than the LIBERTY0 model when simulating healthy, slightly, and moderately damaged needle reflectance within 400–800 nm. The better performance of LIBERTY1 over LIBERTY0 in modeling needle reflectance indicates that calibration of the pigment absorption coefficient improved the accuracy. Both the LIBERTY0 and LIBERTY1 models failed in reflectance simulation of severely damage needles; while LIBERTY2 effectively simulated needles with all levels of damage. The result of validation data shows that the accuracy of reflectance simulation has an MRE between 18.1% and 72% for LIBERTY0, between 2.88% and 71.8% for LIBERTY2, and between 1.4% and 18% for LIBERTY2. This result indicates that the linear mixing method is robust in damaged needle reflectance simulation. Furthermore, for each sample of needles, we can obtain the best-simulated reflectance and corresponding YI when inverting LIBERTY2. The result shows that the derived YI (healthy: 0 < YI < 0.1, slight: 0.05 < YI < 0.15, moderate: 0.15 < YI < 0.5, severe: 0.5 < YI < 1) can be used to separate the different degrees of needle damage from the leaf reflectance spectra. If using YI values, the moderate to severe degrees of damage can be distinctly separated from healthy needles or those with a slight degree of damage. However, using the YI, it is still difficult to detect the difference between healthy and slight degrees of damage in the range of 0.05 < YI < 0.1, when the needle damage is too weak.

3.2. Leaf Chlorophyll Content (LCC) and Leaf Area Index (LAI) Retrieval

Based on the LUT strategy and two regularization options, the performances of LCC and LAI retrieval are presented in Figure 6 using the INFORM model coupled with LIBERTY0, LIBERTY1, and LIBERTY2, respectively. Figure 6 shows the accuracies of LCC and LAI retrieval using different CFs and regularization options. It can be seen that different models may achieve the best retrieval with different CFs, which can be explained by different CFs dealing with different types of error distributions [56]. For example, the best CFs are L-divergence for LIBERTY0-INFORM (at 8% multiple solutions and 6% noise), normal distribution-LSE for LIBERTY1-INFORM (at 3% multiple solutions and 5% noise) and least absolute error for LIBERTY2-INFORM (at 1% multiple solutions and 2% noise). The accuracy of LCC retrieval by LIBERTY1-INFORM (R² = 0.44, RMSE = 58.87 mg m⁻², NRMSE = 26.40% and IOA = 0.65) is higher than that of LIBERTY0-INFORM (R² = 0.34, RMSE = 76.44 mg m⁻², NRMSE = 34.23% and IOA = 0.57). The better performance of LIBERTY1-INFORM on LCC retrieval confirms that calibration improved the estimation accuracy. LIBERTY2-INFORM has the highest accuracy (R² = 0.49, RMSE = 56.24 mg m⁻², NRMSE = 25.22% and IOA = 0.72) in LCC retrieval. However, inversed LCC was systematically underestimated above 300 mg m⁻², while being overestimated below 300 mg m⁻².

Similar to LCC retrieval performance, LIBERTY2-INFORM produced better accuracy of LAI (R² = 0.73, RMSE = 0.17 m m⁻², NRMSE = 11.41% and IOA = 0.91) than LIBERTY0-INFORM (R² = 0.38, RMSE = 0.43 m m⁻², NRMSE = 28.85% and IOA = 0.68) and LIBERTY1-INFORM (R² = 0.63, RMSE = 0.23 m m⁻², NRMSE = 15.44% and IOA = 0.86). However, the NRMSE of LAI retrieval is lower than that of LCC.

3.3. Estimation of Plot Shoots Damage Ratio (SDR)

Based on the above results, the LIBERTY2-INFORM model was used to create LUT and inverse LCC, LAI, plot SDR and YI over the 34 plots. The inversed YI is positively related to plot SDR (Figure 7), and the R² was 0.40 and the RMSE was 0.15.

4. Discussion

4.1. The Contribution of Heterogeneous Leaf Reflectance Simulation

In this paper, three different versions of LIBERTY (LIBERTY0, LIBERTY1, and LIBERTY2) were used to model the reflectance of heterogeneous needles. LIBERTY1 performed better than LIBERTY0 (shown in Figure 5) in simulating needle reflectance for needles were healthy to moderately degree damaged. This result confirmed that calibration of the pigment absorption coefficient could help improve the accuracy of modeling needle reflectance, which was consistent with previous research [64]. However, both LIBERTY0 and LIBERTY1 were out-performed for modeling the reflectance of severely damaged needles. There were two major reasons for this: (1) the absorption coefficient was constant in this study, whereas, in reality, the mix of light-absorbing pigments might change in response to declining LCC or localized chlorosis [48]; and (2) a single pigment did not quantify individual pigments, which provide specific services (such as giving leaves their yellow or red color) to leaves and absorb different amounts of radiation [65,70]. Depending on the cause of chlorosis, there may be concurrent change to multiple pigments such as carotenoids and anthocyanin that would alter reflectance. However, the LIBERTY model only considers a total single pigment instead of separating it into specific individual pigments. Additionally, LIBERTY did not consider the anthocyanins, which were closely associated with the color of the needles [71]. Therefore, LIBERTY0 and LIBERTY1 could not accurately model the reflectance of yellow or red needles.

In this study, the spectral linear mixing method was used to describe the effect of other pigments which were correlated with chlorosis on needle reflectance. Thus, LIBERTY2 improved the simulation of reflectance accuracy for damaged needles. The spectral linear method was a semi-empirical method, it cannot depict the mix of light-absorbing multiple pigments in response to certain forms of stress from a radiative transfer modeling equation. Further, cell structure parameters such as cell diameter (D) and intercellular airspace (xu) would decrease during the beetle attack [27,80]. Because the LIBERTY model requires prior knowledge of D and xu to invert reflectance values [63], accurate cell structure values might improve the accuracy of modeling the reflectance of damaged needles.

4.2. LCC, LAI, and Plot Shoots Damage Ratio Estimation (SDR) Performance

The inversion accuracies of LCC, LAI, and plot SDR were used to compare the capabilities of INFORM coupled with three versions of the LIBERTY model. The LIBERTY2-INFORM has the best accuracy. These results indicate that accurate leaf-scale RT can improve the parameter retrieval performance of the canopy model [64]. The good performance of LAI retrieval is because INFORM effectively considers the sparse canopy structure, understory vegetation, and shadow components [81]. However, the performance of LCC retrieval was worse, although the non-uniform distribution of the chlorophyll impact canopy reflectance has been considered. At the plot scale, LAI can be measured using a suitable optical instrument (e.g., LAI-2000/2200). However, there are no direct approaches or instruments to measure the LCC of a forest plot. Instead, the random sampling method is used to estimate the average value and represent the mean LCC of the plot. The sampling error led to uncertainty in the LCC inversion. In this study, the accuracy of the LCC of the plot depended on the correct classification of the shoot damage class. A more exhaustive classification of shoot damage may improve LCC inversion accuracy.

As for the RT model, 1-dimensional (1D) models are commonly suited to a relatively uniform canopy with a vertical variation of structure but very limited in depicting the complex variation of chlorophyll content in the horizontal structure. INFORM is an extended 1D model, which uses the crown diameter and stem density to describe the horizontal variation of the forest canopy. However, the canopy reflectance simulation was very sensitive to the variations of crown diameter according to the sensitivity analysis of FLIM [82], which added to the inversion error. Hence, the inversion accuracy of LCC is difficult to improve in 1D models. In the future, we will consider using a 3-dimensional (3D) canopy model which can rapidly simulate (e.g., RAPID [83]), to allow the defined 3D scene to invert LCC in PSB stress.

A new index, the YI, can be inversed simultaneously with LCC and LAI. This index theoretically equals plot SDR, which is also evaluated in Figure 7 (R² = 0.40). An unexpected result is the systematic underestimation of YI in damaged plots (plot SDR > 0.2). This may be caused by two major factors: (1) the effect of the clumping of needles or shoots on canopy reflectance; and (2) the overestimation of field SDR because the shoot counting is empirical. As long as a shoot has several yellow needles, it will be counted as a damaged shoot. Conversely, YI is non-zero in healthy plots (plot SDR = 0). This can be explained by the underestimation of LCC (LCC > 300 mg m⁻²) by LIBERTY2-INFORM, influenced by the green understory. The green understory will strongly impact YI effectiveness in the accuracy of predicting plot SDR because of the contribution of the understory to the forest reflectance in a sparse canopy. Therefore, the YI will be more effective in a dense canopy with at least slight damage. Furthermore, the spatial arrangement of damaged needles is another important impact factor. For Yunnan pine, the stressed canopy reflectance is distinctly different from that of a healthy canopy when the damaged needles are at the top part of the crown. In that case, plot SDR can be effectively estimated by the YI. However, when the damaged needles are on the bottom (understory) part of the crown, stressed canopy reflectance is similar with healthy (or slight damage) canopy reflectance because the damaged needles are blocked by healthy needles in a given observation direction. The plot SDR will, therefore, be seriously underestimated when damaged needles are on the bottom (understory) of the crown. A 3D canopy model will be considered in the future to simulate the spatial arrangement of the damaged needles/shoots.

4.3. Role of LUT Setting

In this paper, we used a uniform distribution and fixed ranges to create LUT for both LCC and LAI. However, this setting influences the accuracy of LCC, LAI and the plot SDR estimation. For example, it may cause an error in the variable estimation if the definition of the parameters is too broad compared to the measurements [84]. Because of the external disturbance factors (e.g., pest and pollution), there are difficulties in determining the suitable distribution of parameters. LCC and LAI have different effects on different wavelength regions. LCC takes the dominant role in the visible region; while LAI variations cause spectral variation over the whole VNIR [85,86]. As a result, the band from B2 to B5 are enough to achieve a good LCC retrieval performance; and adding bands may introduce errors. As for LAI, adding bands logically improved the inversion performance [55,56].

5. Conclusions

In this study, the spectral linear mixing method was used in LIBERTY2 which improved the accuracy of modeling heterogeneous needle reflectance. The coupled model LIBERTY2-INFORM was used to quantitatively assess PSB stress in Yunnan pine forests. The YI was introduced in the RT to predict plot SDR representing PSB stress. The YI was used to model the non-uniform distribution of chlorophyll within the needles. The results show significant improvement of needle reflectance simulation under severe damage conditions. Inversion of INFORM-LIBERTY2 using the LUT inversion method leads to the best LAI (R² = 0.73, RMSE = 0.17 m m⁻², NRMSE = 11.41% and IOA = 0.92) and LCC (R² = 0.49, RMSE = 56.24 mg m⁻², NRMSE = 25.22% and IOA = 0.72) retrieval performance for Yunnan pine forests under PSB stress.

Our study implies that a more accurate simulation of leaf optical properties can help enhance the performance of biophysical parameters inversion. The YI is a simple and potential indicator to estimate plot SDR in order to assess PSB outbreaks. To improve the inversion accuracy of the YI, LCC, and LAI, a 3D RT model will be used in the next study.