Assimilating Multiresolution Leaf Area Index of Moso Bamboo Forest from MODIS Time Series Data Based on a Hierarchical Bayesian Network Algorithm

The highly accurate multiresolution leaf area index (LAI) is an important parameter for carbon cycle simulation for bamboo forests at different scales. However, current LAI products have discontinuous resolution with 1 km mostly, that makes it difficult to accurately quantify the spatiotemporal evolution of carbon cycle at different resolutions. Thus, this study used MODIS LAI product (MOD15A2) and MODIS reflectance data (MOD09Q1) of Moso bamboo forest (MBF) from 2015, and it adopted a hierarchical Bayesian network (HBN) algorithm coupled with a dynamic LAI model and the PROSAIL model to obtain high-precision LAI data at multiresolution (i.e., 1000, 500, and 250 m). The results showed the LAIs assimilated using the HBN at the three resolutions corresponded with the actual growth trend of the MBF and correlated significantly with the observed LAI with a determination coefficient (R2) value of >0.80. The highest-precision assimilated LAI was obtained at 1000-m resolution with R2 values of 0.91. The LAI assimilated using the HBN algorithm achieved better accuracy than the MODIS LAI with increases in the R2 value of 2.7 times and decreases in the root mean square error of 87.8%. Therefore, the HBN algorithm applied in this study can effectively obtain highly accurate multiresolution LAI time series data for bamboo forest.


Introduction
The leaf area index (LAI) is defined as one-half of the total green leaf area per unit horizontal ground area [1], and is a requirement in a variety of ecological applications [2].The spatially and temporally distributed LAI data is of crucial importance for researches on carbon and water cycling and on the energy exchange of terrestrial vegetation [3,4].In addition, variations in LAI time series data always considered one of the physiological parameters and indicators to reflect the status of vegetation growth [5,6]; it is therefore a significant land surface parameter for quantitative analysis of many physical and biological processes related to vegetation dynamics [7].
Remote sensing technology is one of methods feasible for obtaining LAI data with large spatial coverage [8].Currently, there are many LAI time series products available, e.g., MODIS [9], CYCLOPES [10], GLOBCARBON [11], and ECOCLIMAP [12].However, because of the effects of cloud cover, snow cover, aerosols, sensor failure, and limitations of the inversion methods [13], many satellite-based LAI products are characterized by low accuracy, high noise, and large errors in their spatiotemporal distributions, which constrained their widespread application [14,15].The development of data assimilation methods, which can incorporate heterogeneous and spatiotemporally discontinuous data [16,17], has led to a series of advances in obtaining accurate spatiotemporal continuous LAI data.The adoption of data assimilation techniques such as the ensemble Kalman filter [18,19], dual ensemble Kalman filter [7], and particle filter [20,21] produced LAI time series data with improved precision, which enhanced the simulation accuracy of the carbon cycle of vegetation [22,23].
All of the above studies show the ability to improve the accuracy of individual LAI products, but they merely work at given single resolution because they are based on the LAI product itself, or they change multisource data with different spatial resolutions into a dataset with uniform resolution via resampling.However, the current LAI products suffer from coarse resolution, and products with multiresolution are scarce.Except for a few LAI products with resolution of 500 m (MOD15A2H/MYD15A2H), the resolutions of most products are ≥1 km with large interval and little diversity between them, such as 3, 5, 8, 10 km [24,25].The discontinuous spatial resolutions make a limitation of using LAI products to accurately describe the carbon cycle at relatively continuous spatial resolution, which make it difficult to comprehensively understand the spatiotemporal evolution and control mechanism of carbon cycle [26,27].Consequently, continuous multiresolution LAI data is highly desirable.
To receive high accuracy multiresolution LAI data, Xiao et al. [28] obtained multiscale LAI data using a multiscale Kalman filter with integrated MODIS and Landsat LAI products; Wang et al. [29] utilized a multiresolution tree combined with a Kalman filter to obtain multiscale LAI data by fusing MISR LAI products with different resolutions; and Jiang et al. [30] developed an ensemble multiscale filter to retrieve multiscale LAI data from satellite observations with different spatial resolutions.Compared with the above methods, a hierarchical Bayesian network (HBN) algorithm also has considerable potential in multiresolution data assimilation with the advantage of not being restricted by linear and Gaussian error distribution hypothesis [16,31].It divide a complex data assimilation model into several relatively simpler models using conditional independence theory, transforming the complex posterior probability into a series of relatively simpler conditional probabilities [32][33][34].Thus, conditional but not complete independence is required in HBN which is more in line with the natural situation.Several studies had promoted the idea and support for the application of HBN in the assimilation of multiresolution data [35,36], and it had been used widely in atmospheric [32,[37][38][39][40], environmental [41][42][43][44], and hydrological research [45,46].
Bamboo forest is a special type of forest in subtropical regions of China, which has the characteristics of rapid growth and high yield [47].Many studies on the carbon cycle of bamboo forest ecosystems, undertaken by domestic and international researchers, have shown that bamboo forests have substantial carbon sequestration capability and play an important role in maintaining the regional ecological environment and carbon balance [48][49][50][51].As one of the most important parameters in simulation of the carbon cycle, multiresolution LAI data plays a big part in quantifying and clarifying the evolution of carbon cycle at different scales.This study used MODIS LAI data (MOD15A2) and MODIS reflectance data (MOD09Q1) of Moso bamboo forest (MBF) in Anji County (Zhejiang Province, China) from 2015, and adopted an HBN algorithm coupled with a LAI dynamic model and the PROSAIL model to assimilate highly accurate LAI time series data with multiresolution (i.e., 1000, 500, and 250 m).The assimilated LAI time series data could be used as spatiotemporal data to enhance the accuracy of carbon cycle multiscale simulation of bamboo forest.

Study Area
In this study, the MBF flux measurement site and the area around it within the range of 1 km were used as the study area and the sampling area of observed LAI, respectively.The location of the MBF flux measurement site and observed area in Shanchuan town, Anji County are shown in Figure 1.

Study Area
In this study, the MBF flux measurement site and the area around it within the range of 1 km were used as the study area and the sampling area of observed LAI, respectively.The location of the MBF flux measurement site and observed area in Shanchuan town, Anji County are shown in Figure 1.E).This area has a subtropical monsoon climate with distinct seasons and abundant rainfall.The MBF is distributed widely across Anji County covering an area of 55,367 ha, which accounts for approximately 45% of the total forested area.The area (1 × 1 km) around the flux tower consists of MBF.The average diameter of the bamboo at breast height is 9.3 cm and the canopy height is 12-18 m with a sparse understory of shrubs and herbs [20].

Processing of Satellite Data
The MODIS 8-day 1000-m LAI product (MOD15A2) and the 8-day 250-m reflectance product (MOD09Q1) for Zhejiang Province in 2015 were downloaded from NASA's website (https://ladsweb.nascom.nasa.gov).The MOD15A2 and MOD09A1 datasets were reprojected to the WGS84 coordinate system using MODIS Reprojection Tools software.ENVI v5.3 software was used to determine those pixels corresponding to the locations of the flux measurement site and to extract the pixel values in the range of 1 km; thus, 1 MOD15A1 pixel value and 16 MOD09Q1 pixel values were extracted in each phase.
The locally adjusted cubic-spline capping (LACC) method was applied to remove noise in the MOD15A2 product, following which the smoothed data were used as initial values for the LAI dynamic model.The LACC method has been described in detail in the literature [52].The smoothed LAI time series will be shown in Section 4.1 and used as a comparison for the results.
The MOD09Q1 product includes seven bands; however, this study used only the red band (RED) and the near infrared band (NIR).To reduce the effect of atmospheric factors on reflectance, the Savitzky-Golay smooth and the method based on normalized difference vegetation index (NDVI) developed by Xiao et al. [53], were used to eliminate outliers.This method has been described previously in the literature [18].Take the first pixel as an example, the NDVI time series envelope and the outliers from the MODIS reflectance data are shown in Figure 2. It is seen from the figure that outliers at DOY 49, 73, 249, 361 and others have been eliminated and fluctuations of original NDVI time series data is obviously decreased.Therefore, the quality of RED and NIR reflectance data was improved.MBF.The average diameter of the bamboo at breast height is 9.3 cm and the canopy height is 12-18 m with a sparse understory of shrubs and herbs [20].

Processing of Satellite Data
The MODIS 8-day 1000-m LAI product (MOD15A2) and the 8-day 250-m reflectance product (MOD09Q1) for Zhejiang Province in 2015 were downloaded from NASA's website (https://ladsweb.nascom.nasa.gov).The MOD15A2 and MOD09A1 datasets were reprojected to the WGS84 coordinate system using MODIS Reprojection Tools software.ENVI v5.3 software was used to determine those pixels corresponding to the locations of the flux measurement site and to extract the pixel values in the range of 1 km; thus, 1 MOD15A1 pixel value and 16 MOD09Q1 pixel values were extracted in each phase.
The locally adjusted cubic-spline capping (LACC) method was applied to remove noise in the MOD15A2 product, following which the smoothed data were used as initial values for the LAI dynamic model.The LACC method has been described in detail in the literature [52].The smoothed LAI time series will be shown in Section 4.1 and used as a comparison for the results.
The MOD09Q1 product includes seven bands; however, this study used only the red band (RED) and the near infrared band (NIR).To reduce the effect of atmospheric factors on reflectance, the Savitzky-Golay smooth and the method based on normalized difference vegetation index (NDVI) developed by Xiao et al. [53], were used to eliminate outliers.This method has been described previously in the literature [18].Take the first pixel as an example, the NDVI time series envelope and the outliers from the MODIS reflectance data are shown in Figure 2. It is seen from the figure that outliers at DOY 49, 73, 249, 361 and others have been eliminated and fluctuations of original NDVI time series data is obviously decreased.Therefore, the quality of RED and NIR reflectance data was improved.

Processing of Observed Data
Ground-observed LAI data were collected in study area for each month.There were 11 periods of observation of LAI data in 2015.Observation in February was missing.The LAI was calculated using the LAI (2000G)-LogCI algorithm in the WinSCANOPY v2009a program (Regent Instruments Inc., Québec, QC, Canada), based on canopy images taken in the field using a digital camera with a fisheye lens (Figure 3) [54].The plot design of the observed LAI is shown in Figure 4 [20].Taking the flux measurement site as the center, 5 observation centers (stars in Figure 4) and 20 observation points (dots in Figure 4) were set across the 1 × 1 km area.Based on the coordinates of the observation centers and observation points, the pixels to which they belonged were determined at 1000-, 500-, and 250-

Processing of Observed Data
Ground-observed LAI data were collected in study area for each month.There were 11 periods of observation of LAI data in 2015.Observation in February was missing.The LAI was calculated using the LAI (2000G)-LogCI algorithm in the WinSCANOPY v2009a program (Regent Instruments Inc., Québec, QC, Canada), based on canopy images taken in the field using a digital camera with a fisheye lens (Figure 3) [54].The plot design of the observed LAI is shown in Figure 4 [20].Taking the flux measurement site as the center, 5 observation centers (stars in Figure 4) and 20 observation points (dots in Figure 4) were set across the 1 × 1 km area.Based on the coordinates of the observation centers and observation points, the pixels to which they belonged were determined at 1000-, 500-, and 250-m resolution, and the average value was taken as the actual observed LAI for the pixels.The order of pixels at 1000-, 500-, and 250-m resolution is shown in flow chart in Figure 5 and the value of observed LAI of corresponding pixels is shown in Table 1.For example, LAI_500_1 indicates the observed LAI corresponding to the first pixel at 500-m resolution.
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 21 m resolution, and the average value was taken as the actual observed LAI for the pixels.The order of pixels at 1000-, 500-, and 250-m resolution is shown in flow chart in Figure 5 and the value of observed LAI of corresponding pixels is shown in Table 1.For example, LAI_500_1 indicates the observed LAI corresponding to the first pixel at 500-m resolution.The leaf bidirectional reflectance (LBR) and canopy bidirectional reflectance (CBR) were measured at a sample culm beside each flux measurement site using a FieldSpec Pro spectroradiometer (Analytical Spectral Devices Inc., Boulder, CO, USA); details concerning the measurement method are available in the literature [18,55].The leaves were transported to the laboratory to measure biochemical parameters for use as inputs to the PROSAIL model.The chlorophyll content (Cab, µg cm −2 ) of the leaves was measured using a spectrophotometer (UV-2102C/PC/PCS, Unico (Shanghai) Instrument Co., Ltd., Shanghai, China) [56].The leaf area (m 2 ) was measured using an AM300 (ADC Bioscientific Ltd., Hertfordshire, UK).The fresh weight and dry matter (after drying for 48 h at 75 °C) were measured using an FA2104 (Shanghai Liangping m resolution, and the average value was taken as the actual observed LAI for the pixels.The order of pixels at 1000-, 500-, and 250-m resolution is shown in flow chart in Figure 5 and the value of observed LAI of corresponding pixels is shown in Table 1.For example, LAI_500_1 indicates the observed LAI corresponding to the first pixel at 500-m resolution.The leaf bidirectional reflectance (LBR) and canopy bidirectional reflectance (CBR) were measured at a sample culm beside each flux measurement site using a FieldSpec Pro spectroradiometer (Analytical Spectral Devices Inc., Boulder, CO, USA); details concerning the measurement method are available in the literature [18,55].The leaves were transported to the laboratory to measure biochemical parameters for use as inputs to the PROSAIL model.The chlorophyll content (Cab, µg cm −2 ) of the leaves was measured using a spectrophotometer (UV-2102C/PC/PCS, Unico (Shanghai) Instrument Co., Ltd., Shanghai, China) [56].The leaf area (m 2 ) was measured using an AM300 (ADC Bioscientific Ltd., Hertfordshire, UK).The fresh weight and dry matter (after drying for 48 h at 75 °C) were measured using an FA2104 (Shanghai Liangping The leaf bidirectional reflectance (LBR) and canopy bidirectional reflectance (CBR) were measured at a sample culm beside each flux measurement site using a FieldSpec Pro spectroradiometer (Analytical Spectral Devices Inc., Boulder, CO, USA); details concerning the measurement method are available in the literature [18,55].The leaves were transported to the laboratory to measure biochemical parameters for use as inputs to the PROSAIL model.The chlorophyll content (Cab, µg cm −2 ) of the leaves was measured using a spectrophotometer (UV-2102C/PC/PCS, Unico (Shanghai) Instrument Co., Ltd., Shanghai, China) [56].The leaf area (m 2 ) was measured using an AM300 (ADC Bioscientific Ltd., Hertfordshire, UK).The fresh weight and dry matter (after drying for 48 h at 75 • C) were measured using an FA2104 (Shanghai Liangping Instrument Co., Ltd., Shanghai, China).The difference between the fresh weight and dry matter was calculated as the water content.

LAI Dynamic Model
In this study, the semiempirical model developed by Dickinson et al. [57] was used as a dynamic model to obtain the simulated LAI.The model and the parameters suited for MBF have been described previously in the literature [18,20,23].The model can be expressed as in the following:

Study Methodology
This study used an HBN algorithm to assimilate multiresolution LAI data.A flow chart of the assimilation process is shown in Figure 5.The steps are as follows: 1.
The 1000-m-resolution MODIS LAI and the 250-m-resolution MODIS reflectance data were processed, which had been carried out in Section 2.2.Smoothed MODIS LAI data were input into the LAI dynamic model to obtain the simulated LAI.

2.
Based on the 250-m-resolution MODIS reflectance data, resolution-specific likelihood inference (RESL) or resolution-specific restricted-likelihood inference (RESREL) method was used to estimate the parameters in the data assimilation.

3.
The LAI data and simulated reflectance data at each scale were obtained using the simulated LAI, parameters in the process model of the HBN, and the PROSAIL model.Combined with MODIS reflectance data, the structure of the HBN at three resolutions (1000, 500, and 250 m) was initialized.4.
An inference of HBN was used to update node (pixel) states by integrating observations with different spatial resolutions.The probability distributions of node states were updated by upward filtering from observations from the finer scale (250 m) to the coarser scale (1000 m).Then, the downward smoothing started from the coarser scale (1000 m) and ended at the finer scale (250 m) to update the posterior probability distributions of all node states according to all of the observations.5.
The Markov chain Monte Carlo (MCMC) sampling method was used to sample the probability distribution.The LAI corresponding to the maximum posterior probability was the assimilation result of the three scales.6.
Because of uncertainty and error in the parameter estimation of the HBN, the above process was repeated 100 times and the average value was used as the final assimilation LAI.The assimilated LAI time series data were compared with the observed LAI.The methodology of this HBN algorithm was implemented within MATLAB.

LAI Dynamic Model
In this study, the semiempirical model developed by Dickinson et al. [57] was used as a dynamic model to obtain the simulated LAI.The model and the parameters suited for MBF have been described previously in the literature [18,20,23].The model can be expressed as in the following: x = (LAI − LAI min )/(LAI max − LAI min ), where LAI t+1 and LAI t represent the LAI at the previous time step (t + 1) and current time step (t), respectively; R(x) is a smoothing function; x is the normalization of the LAI; LAI max and LAI min are the maximum and minimum of the LAI in a year, respectively; L 0 represents the max LAI; L t represents the rate of leaf litter; and L 0 , L t , λ 0 , and c are parameters determined based on experience.In this study, the value of c was set as 0.5 and the initial values of the other parameters were obtained according to observed LAI fitting.

PROSAIL Model
Forest CBR is particularly sensitive to LAI [58].Therefore, the LAI simulated by the dynamic model was used as the input for the PROSAIL model to achieve simulation of bamboo CBR.The PROSAIL model, which is a combination of the PROSPECT5 and 4SAIL models, is an advanced model for simulation of CBR spectra.The PROSPECT5 model simulates LBR [59], and the 4SAIL model simulates CBR [60].The parameters of PROSAIL model were optimized using the measured LBR and CBR, and biochemical parameters of the observed area.The optimized parameters are shown in Table 2.The reader is referred to the literature for specific details of PROSAIL model and optimization methods [54,[59][60][61].

Hierarchical Bayesian Network Algorithm
In this study, the HBN algorithm was used to obtain LAI time series data with different resolutions for MBF.The algorithm mainly includes the steps of construction of the HBN, inference of the HBN, and MCMC sampling.To simplify the model construction, a transition scale using given data with 500-m resolution was constructed between the MODIS LAI data with 1000-m resolution and MODIS reflectance data with 250-m resolution.This represents a dual relationship between the two datasets on adjacent scales, i.e., one parent node corresponds to four child nodes, as shown in Figure 5.

Construction of the HBN
The HBN include a data model, process model, and parameter model [43,44].The data, processes, and parameters are defined as three random variables and establish conditional probability distribution models, such that the data assimilation problem is transformed into updating the posterior probability distribution of the processes and parameters under the condition of the existing data [31,45].

Data model
The data model defines the conditional probability distribution model between the observed data and real processes, related parameters [46].In this study, data with three resolutions gradually increasing from top to bottom were included.The first scale considered LAI data with 1000-m resolution, the second scale had no data, and the third scale comprised reflectance data with 250-m resolution.The data model for each scale was defined based on Equations ( 5)-( 7): Second scaleNo data, No model, ( 6) where Z ch(i,r) and Y ch(i,r) represent the observed data and real process, respectively, for all child nodes of parent node (i, r), i.e., Z ch(i,r

Process model
The process model describes the spatiotemporal distribution and variation of parameters, and it defines the conditional dependency between real processes at different scales [46,62].In this study, the joint posterior probability inference of the real process was transformed into a certain-scale posterior probability inference by defining the conditional dependency between different scales [63].The process model is defined by Equation (8): where Y ch(i,r) represents the real process of node ch(i, r), which obeys the Gaussian distribution is the real process of node (i, r); and w ch(i,r) represents the error of the process model that obeys the Gaussian distribution w ch(i,r) ∼ N 0, σ 2 W ch(i,r) .The process model in Equation ( 8) does not have one-to-one correspondence between Y ch(i,r) and Y (i,r) , w ch(i,r) .Here, Y ch(i,r) is a vector of length m r , whereas Y (i,r) , w ch(i,r) has a total of m r + 1 elements.However, by placing a single linear constraint on the error term w ch(i,r) , a one-to-one correspondence is achieved [64]: To satisfy Equation ( 9), we let Q (i,r) be any m r × (m r − 1) orthonormal matrix with columns that span the space orthogonal to q (i,r) (q (i,r) Q (i,r) = 0 and Q (i,r) Q (i,r) = I).The matrix Q (i,r) can be taken as the Helmert matrix orthogonal to q ch(i,r) [65].Given a Q (i,r) matrix, any w ch(i,r) that satisfies Equation ( 9) can be written as follows: and then Equation ( 8) can be written as follows: where w * (i,r) obeys the Gaussian distribution w * (i,r) ∼ N 0, σ 2 W * (i,r) , and W ch(i,r) = Q (i,r) W * (i,r) Q (i,r) .By defining q (i,r) = a ch(i,r) = a ch 1 (i,r) , a ch 2 (i,r) , • • • , a ch mr (i,r) , a ch j (i,r) represents the area of the j-th child of node (i, r).

Parameter model
The parameter model defines the prior distribution of all parameters [32,62].In this study, the parameters in the data model and the process model included σ, V ch(i,r) , and W * (i,r) .σ was a constant, while V ch(i,r) and W * (i,r) were parameters for Gaussian distribution of v ch(i,r) in Equations ( 5), ( 7) and w * (i,r) in Equation (10).The value of σ was set as 10 −2 ; the values of V ch(i,r) in Equations ( 5) and ( 7) were calculated by covariance of MODIS LAI and reflectance data, respectively.W * (i,r) was estimated based on Q (i,r) in Equation ( 10) using RESL [36] and RESREL [66]; the reader is referred to the literature for specific details [63,64].

Inference of the HBN
The inference of the HBN includes two steps: upward filtering and downward smoothing [67].For upward filtering, observations information is transmitted from the finer-resolution with 250 m (the bottom scale) to the coarser-resolution with 1000 m (the top scale) to predict and update the probability distribution of nodes using finer-scale observations.Conversely, downward smoothing transmits information from the coarser scale (1000 m) to the finer scale (250 m) to update the posterior probability distributions of all node states according to all of the satellite observations [64,68].After the above two steps, each node integrates multiple observations with different spatial resolutions, and the final posterior probability distribution is obtained; the reader is referred to the literature for specific details [63].

Upward filtering
The initial probability distribution of the data at the top scale is defined as Y (i,0) ∼ N(µ 0 , Σ 0 ).Assuming that σ, V ch(i,r) , and W * (i,r) are known, the probability distribution of all nodes is Y (i,r) ∼ N µ (i,r) , Σ (i,r) , where the mean µ (i,r) is the same for all nodes and the variance is Σ (i,r) = W (i,r) Σ pa(i,r) W (i,r) .The steps of upward filtering are as follows [63]: Step 1: Calculate the probability distribution of the data at the bottom scale If there are no observed data at the bottom scale (at the R-1 scale): If there are observed data at the bottom scale (at the R-1 scale): where Y and Z have different meanings at different resolutions and the specific contents are detailed in the data model (as is the same below).
Step 2: Calculate the posterior probability distribution of all nodes Step 3: Predict probability distribution for parent node pa(i, r) whose child node is (i, r) where Z * pa(i,r) represents the observation data of the child nodes that do not contain observation data of parent node pa(i, r).
Step 4: Predict probability distribution of observation of parent node pa(i, r) whose child node is (i, r) Z pa(i,r) |Z * pa(i,r) ∼ N f pa(i,r) , Q pa(i,r) , ( 16) Q pa(i,r) = F pa(i,r) R pa(i,r) F pa(i,r) , (18) where F pa(i,r) represents the PROSAIL model and the specific contents are detailed in the data model.
Step 5: Calculate the posterior probability distribution by correcting the prior probability of the observation of parent node pa(i, r).
Supposing Y pa(i,r) |Z pa(i,r) ∼ N m pa(i,r) , C pa(i,r) , then: If there are no observed data at node pa(i, r): If there are observed data at node pa(i, r): e pa(i,r) = R pa(i,r) F pa(i,r) Q −1 pa(i,r) .
Step 6: Repeat the above steps until at the top scale.

Downward smoothing
The posterior probability distribution at the top scale obtained from upward filtering is used as the initial probability distribution for downward smoothing.The steps of downward smoothing are as follows [63]: Step 1: Obtain the probability distribution of the data at the top scale from upward filtering where Y and Z have different meanings at different resolutions and the specific contents are detailed in the data model (as is the same below).
Step 2: Calculate the probability distribution of all nodes Step 3: Calculate the prior probability distribution of any node at any scale Step 4: Calculate the posterior probability distribution of node (i, r) Y (i,r) |Z (i,0) ∼ N s (i,r) , S (i,r) , ( 30) Step 5: Repeat the above steps until at the bottom scale.

MCMC Sampling
After obtaining the posterior probability distribution using the HBN algorithm, a sampling method was used to obtain the LAI value according to the probability distribution [37,41].Currently, Markov Chain Monte Carlo (MCMC) sampling methods are widely used for such purposes [69,70] and these methods include the Metropolis-Hastings, Gibbs, and Slice sampling [44,71].In this study, the Metropolis-Hastings method was used to sample the LAI time series data at different scales; the reader is referred to the literature for specific details [72][73][74].

Accuracy Assessment
The determination coefficient (R 2 ), root mean square error (RMSE), and absolute bias (aBIAS) were used to evaluate the accuracy of the assimilated LAI.Generally, higher values of R 2 and lower values of both RMSE and aBIAS indicate better accuracy.The RMSE and aBIAS were calculated as follows: where N indicates the number of samples and it is 11 in this study because of 11 available observed LAI; y 0 indicate the values of observed LAI; and y i represent the values of assimilated LAI.

Validation of LAI Assimilation at 1000-m Resolution
The different LAI time series data of MBF at 1000-m resolution are shown in Figure 6.These LAI data include the observed LAI (Field_LAI), original MODIS LAI (MODIS_LAI), smoothed MODIS LAI by LACC (LACC_LAI), and assimilated LAI using the HBN (HBN_LAI_1000) in 2015.
The MODIS_LAI fluctuated considerably (0-6) during the growing season with frequent outliers with low value and errors.Although the LACC_LAI was smoother than MODIS_LAI and the outliers and fluctuations were obviously decreased, it was significantly different from the Field_LAI.The MBF HBN_LAI_1000 showed a trend of slow growth from 3.3 to 4.3 in spring (March-May) before reaching its annual maximum value of 5.6 in summer (June-August).It then decreased gradually from 4.3 to 3.2 in autumn (September-November) with the annual minimum value of 2.3 occurring in winter (December-February).This trend corresponded with the actual trend of MBF LAI.The HBN_LAI_1000 was correlated significantly with the Field_LAI with an R 2 value of 0.91 and RMSE value of 0.27.However, it was slightly higher than the Field_LAI at DOY 200-320 with a range of overestimation of 0.2-0.4.Analysis of the relationships of both the MODIS_LAI and the HBN_LAI_1000 with the Field_LAI revealed the R 2 of the HBN_LAI_1000 and Field_LAI was 2.7 times higher and the RMSE was 87.8% lower compared with the MODIS_LAI and Field_LAI.This indicates the assimilated LAI obtained by the HBN algorithm greatly improves the accuracy and reduces the error of the MODIS_LAI, producing better correspondence with the actual LAI data.Thus, the HBN_LAI_1000 can reflect the dynamic changes of MBF LAI over a long time series.

LAI;
indicate the values of observed LAI; and represent the values of assimilated LAI.

Validation of LAI Assimilation at 1000-m Resolution
The different LAI time series data of MBF at 1000-m resolution are shown in Figure 6.These LAI data include the observed LAI (Field_LAI), original MODIS LAI (MODIS_LAI), smoothed MODIS LAI by LACC (LACC_LAI), and assimilated LAI using the HBN (HBN_LAI_1000) in 2015.The MODIS_LAI fluctuated considerably (0-6) during the growing season with frequent outliers with low value and errors.Although the LACC_LAI was smoother than MODIS_LAI and the outliers and fluctuations were obviously decreased, it was significantly different from the Field_LAI.The MBF HBN_LAI_1000 showed a trend of slow growth from 3.3 to 4.3 in spring (March-May) before reaching its annual maximum value of 5.6 in summer (June-August).It then decreased gradually from 4.3 to 3.2 in autumn (September-November) with the annual minimum value of 2.3 occurring in winter (December-February).This trend corresponded with the actual trend of MBF LAI.The HBN_LAI_1000 was correlated significantly with the Field_LAI with an R 2 value of 0.91 and RMSE

Validation of LAI Assimilation at 500-m Resolution
The time series data of observed LAI (Field_LAI) and assimilated LAI (HBN_LAI_500) of MBF at 500-m resolution in 2015 are shown in Figure 7. Traces 1-4 are the LAI assimilation results of the four pixels with 500-m resolution that correspond to one pixel at 1000-m resolution.
The HBN_LAI_500 of the four pixels showed a trend of growth in spring (3.3-4.1, 3.5-4.0,3.4-4.2,3.3-4.2) before reaching maximum values in summer (6.1, 6.3, 6.3, 6.1), decreasing in autumn (4.9-3.2, 4.5-3.2,4.6-3.2,4.5-3.2),and reaching minimum values in winter (2.4,2.5, 2.4, 2.4).This trend, which corresponded with the actual trend of LAI, was consistent with the HBN_LAI_1000 in MBF; however, compared with the HBN_LAI_1000, more fluctuations were evident in the HBN_LAI_500.The HBN_LAI_500 of the four pixels correlated significantly with the Field_LAI; the R 2 values for pixels 1-4 were 0.89, 0.90, 0.93, and 0.87, respectively with RMSE values of < 0.50 and aBIAS values of < 0.40.However, the HBN_LAI_500 was obviously overestimated in comparison with the Field_LAI during the growing period (DOY 200-250) with overestimations of 1.2, 1.0, 1.1, and 0.8 for pixels 1-4, respectively, similar to the HBN_LAI_1000.Although the accuracy of the MBF HBN_LAI_500 is shown slightly lower than the HBN_LAI_1000, it still has reasonably high accuracy and it corresponds with the actual LAI data.

Validation of LAI Assimilation at 250-m Resolution
The time series data of observed LAI (Field_LAI) and assimilated LAI (HBN_LAI_250) of MBF at 250-m resolution in 2015, and the correlations between them, are shown in Figures 8 and 9. Traces 1-16 are the LAI assimilation results of the 16 pixels with 250-m resolution that correspond to one pixel at 1000-m resolution.
The HBN_LAI_250 of the 16 pixels was consistent with the actual trend of LAI, although considerable fluctuation was evident, especially at DOY 0-200.During the growing period in summer (DOY 200-250), the HBN_LAI_250 of eight pixels (4, 6, 7, 8, 9, 10, 11, and 14) were obviously overestimated in comparison with the Field_LAI with a range of overestimation of 0.4-1.0.For the HBN_LAI_250 and Field_LAI (Figure 9), the R 2 values were >0.80 (the R 2 value for the fourth pixel was 0.92), RMSE values were <0.50, and aBIAS values were <0.40.This indicates that the HBN_LAI_250 for MBF has high accuracy and it is an effective assimilation result of multiresolution LAI using HBN algorithm.

Validation of LAI Assimilation at 500-m Resolution
The time series data of observed LAI (Field_LAI) and assimilated LAI (HBN_LAI_500) of MBF at 500-m resolution in 2015 are shown in Figure 7. Traces 1-4 are the LAI assimilation results of the four pixels with 500-m resolution that correspond to one pixel at 1000-m resolution.

Validation of LAI Assimilation at 250-m Resolution
The time series data of observed LAI (Field_LAI) and assimilated LAI (HBN_LAI_250) of MBF at 250-m resolution in 2015, and the correlations between them, are shown in Figures 8 and 9. Traces 1-16 are the LAI assimilation results of the 16 pixels with 250-m resolution that correspond to one pixel at 1000-m resolution.The HBN_LAI_250 of the 16 pixels was consistent with the actual trend of LAI, although considerable fluctuation was evident, especially at DOY 0-200.During the growing period in summer (DOY 200-250), the HBN_LAI_250 of eight pixels (4, 6, 7, 8, 9, 10, 11, and 14) were obviously overestimated in comparison with the Field_LAI with a range of overestimation of 0.4-1.0.For the HBN_LAI_250 and Field_LAI (Figure 9), the R 2 values were >0.80 (the R 2 value for the fourth pixel was 0.92), RMSE values were <0.50, and aBIAS values were <0.40.This indicates that the HBN_LAI_250 for MBF has high accuracy and it is an effective assimilation result of multiresolution

Discussion and Conclusions
Several researches had achieved the acquisition of multiresolution LAI using multiresolution tree combined with a Kalman filter [28,29], and had validated the precision comparing with high-resolution LAI [30].This study developed an HBN algorithm coupled with a LAI dynamic model and the PROSAIL model to assimilate LAI data with different resolutions (1000, 500, and 250 m) of MBF.The results showed that the multiresolution HBN_LAI corresponded with the actual trend of growth of the MBF and that it was correlated significantly with the Field_LAI (R 2 values reached >0.80).The R 2 value between the HBN_LAI and Field_LAI at 1000-m resolution was 0.91 and the RMSE value was 0.27; at 500-m resolution, the averages of the R 2 values and the RMSE values 0.90 and 0.42; and at 250-m resolution, the averages of the R 2 values and the RMSE values were 0.86 and 0.35.The HBN_LAI at the coarsest scale (1000-m resolution) provided better accuracy, which was consistent with the finding in literature [30].Compared with the MODIS_LAI, the R 2 value of the HBN_LAI_1000 increased by 2.7 times, while the RMSE value decreased by 87.8%.Compared with other researches on improving the precision of MODIS LAI products [7,18,23], this study had achieved better results.Therefore, the HBN algorithm significantly improved the accuracy and reduced the error of the MODIS LAI products, and accurate LAI at different scales could be obtained.
Analysis of the results revealed that the HBN_LAI_1000, HBN_LAI_500, and HBN_LAI_250 obviously overestimated with 7.5%, 20.0%, and 11.3%, respectively, at DOY 217.This overestimation can be explained in terms of the following three aspects.
First, differences between MODIS reflectance and the reflectance simulated by the PROSAIL model will affect the assimilation results [18,22].As there are no MODIS reflectance data available at 1000-and 500-m resolution, the reflectance at 1000-and 500-m resolution is transmitted from the MODIS reflectance at the bottom scale of 250-m resolution, which is different from the simulated reflectance.This will have an impact on the HBN_LAI_1000 and HBN_LAI_500 and overestimation will occur.At 250-m resolution, the difference between MODIS reflectance and simulated reflectance transmitted by the process model has an influence on the HBN_LAI_250.Figure 10 shows the RED and NIR bands of MODIS reflectance and simulated reflectance of 16 pixels at 250-m resolution.The MODIS reflectance data of the 16 pixels has considerable fluctuation and some outliers (hollow circles in Figure 10), while the reflectance simulated by the PROSAIL model is more concentrated and has less fluctuation than MODIS reflectance.In the RED band, the mean value of MODIS reflectance is higher than the simulated reflectance.The maximum MODIS reflectance is up to 1.7 times that of the simulated reflectance.In contrast, in the NIR band, the mean value of simulated reflectance is higher than MODIS reflectance.The minimum MODIS reflectance is 67% of the simulated reflectance.In the hierarchical Bayesian inference, the LAI with 500-and 250-m resolution is inferred from the 1000-m-resolution LAI and the parameter w ch(i,r) in process model, causing an overestimation of LAI at the 500-and 250-m scale.Furthermore, there is considerable difference between the simulated reflectance and MODIS reflectance, and the HBN_LAI_250 is overestimated.Yang et al. [75], Li et al. [18] and Li et al. [22] have all highlighted that error between the simulated reflectance and MODIS reflectance is one of the primary reasons for inaccuracy in LAI estimations.
Second, the uncertainty of the estimated parameter w ch(i,r) in the process model will also cause errors in the assimilation results.In the HBN algorithm, RESL and RESREL used to estimate parameter W * (i,r) by using MODIS reflectance data, and then parameter w ch(i,r) was obtained using Equation (6).In the hierarchical Bayesian inference, the LAI with 500-and 250-m resolution is inferred from the 1000-m-resolution LAI and parameter w ch(i,r) in the process model, uncertainty in the estimated parameter w ch(i,r) causes errors of LAI at 500-and 250-m resolution, which in turn affects the assimilation results.In this study, 100 iterations of the HBN were performed to reduce the errors.Figure 11 shows the mean value and standard deviation (SD) of parameter w ch(i,r) for MBF at the 1000-500-m and 500-250-m scales, respectively.The mean value of w ch(i,r) is maintained between −0.10 and 0.10 at the 1000-500-m scale, except for some nodes with outliers.However, at DOY 230-290, w ch(i,r) is slightly larger than that at other periods, and the SD of w ch(i,r) is very large, indicating that the uncertainty is substantial during this period.The SD for the first, second, and third child node is 6, 5, and 10 times greater, respectively, than the mean value.At the 500-250-m scale, considerable fluctuations are evident in both w ch(i,r) and its SD and the mean value of w ch(i,r) is maintained between −0.20 and 0.20.Similar to the 1000-500-m scale, w ch(i,r) at DOY 230-290 is larger than that at other periods with greater uncertainty, which could be related to the overestimation of the HBN_LAI for MBF at DOY 217.
will occur.At 250-m resolution, the difference between MODIS reflectance and simulated reflectance transmitted by the process model has an influence on the HBN_LAI_250.Figure 10 shows the RED and NIR bands of MODIS reflectance and simulated reflectance of 16 pixels at 250-m resolution.The MODIS reflectance data of the 16 pixels has considerable fluctuation and some outliers (hollow circles in Figure 10), while the reflectance simulated by the PROSAIL model is more concentrated and has less fluctuation than MODIS reflectance.In the RED band, the mean value of MODIS reflectance is higher than the simulated reflectance.The maximum MODIS reflectance is up to 1.7 times that of the simulated reflectance.In contrast, in the NIR band, the mean value of simulated reflectance is higher than MODIS reflectance.The minimum MODIS reflectance is 67% of the simulated reflectance.In the hierarchical Bayesian inference, the LAI with 500-and 250-m resolution is inferred from the 1000-mresolution LAI and the parameter ( , ) in process model, causing an overestimation of LAI at the 500-and 250-m scale.Furthermore, there is considerable difference between the simulated reflectance and MODIS reflectance, and the HBN_LAI_250 is overestimated.Yang et al. [75], Li et al. [18] and Li et al. [22] have all highlighted that error between the simulated reflectance and MODIS reflectance is one of the primary reasons for inaccuracy in LAI estimations.Second, the uncertainty of the estimated parameter ( , ) in the process model will also cause errors in the assimilation results.In the HBN algorithm, RESL and RESREL used to estimate parameter ( , ) * by using MODIS reflectance data, and then parameter ( , ) was obtained using Equation ( 6).In the hierarchical Bayesian inference, the LAI with 500-and 250-m resolution is inferred from the 1000-m-resolution LAI and parameter ( , ) in the process model, uncertainty in the estimated parameter ( , ) causes errors of LAI at 500-and 250-m resolution, which in turn affects the assimilation results.In this study, 100 iterations of the HBN were performed to reduce the errors.Figure 11 shows the mean value and standard deviation (SD) of parameter ( , ) for MBF at the 1000-500-m and 500-250-m scales, respectively.The mean value of ( , ) is maintained between −0.10 and 0.10 at the 1000-500-m scale, except for some nodes with outliers.However, at DOY 230-290, ( , ) is slightly larger than that at other periods, and the SD of ( , ) is very large, indicating that the uncertainty is substantial during this period.The SD for the first, second, and third child node is 6, 5, and 10 times greater, respectively, than the mean value.At the 500-250-m scale, considerable fluctuations are evident in both ( , ) and its SD and the mean value of ( , ) is maintained between −0.20 and 0.20.Similar to the 1000-500-m scale, ( , ) at DOY 230-290 is larger than that at other periods with greater uncertainty, which could be related to the overestimation of the HBN_LAI for MBF at DOY 217.Finally, the overestimation of the HBN_LAI_500 for MBF is most serious because of the lack of original observation data at 500-m resolution.The inference of the HBN and update of the LAI posterior probability at 500-m scale rely on MODIS LAI products with 1000-m resolution and MODIS reflectance data with 250-m resolution, in which the bidirectional accumulation of errors is unavoidable.Therefore, overestimation is more serious in the HBN_LAI_500.To avoid this phenomenon, the application of multisource remote sensing data, e.g., Landsat data, could be considered to provide multiresolution data for the HBN assimilation.Moreover, the complex structure could be taken into account to flexibly adjust the relationship between parent-child nodes in the HBN and to describe irregular shapes of geographical regions [30].
In this study, the HBN algorithm was developed to obtain high precision LAI time series data from MODIS data at different spatial resolutions (1000, 500, and 250 m).The results demonstrated that the assimilated LAI values at each spatial resolution were in good agreement with the actual trend of MBF LAI and the assimilated LAI at the coarsest scale (1000 m) provided best accuracy.The accuracy of assimilated LAI was significantly improved compared with MODIS LAI products.Therefore, the HBN algorithm applied in this study can effectively obtain highly accurate multiresolution LAI time series data for bamboo forest.This study achieved multiresolution LAI time series data using HBN algorithm at MBF flux measurement site.In the near future, we will extend Finally, the overestimation of the HBN_LAI_500 for MBF is most serious because of the lack of original observation data at 500-m resolution.The inference of the HBN and update of the LAI posterior probability at 500-m scale rely on MODIS LAI products with 1000-m resolution and MODIS reflectance data with 250-m resolution, in which the bidirectional accumulation of errors is unavoidable.Therefore, overestimation is more serious in the HBN_LAI_500.To avoid this phenomenon, the application of multisource remote sensing data, e.g., Landsat data, could be considered to provide multiresolution data for the HBN assimilation.Moreover, the complex structure could be taken into account to flexibly adjust the relationship between parent-child nodes in the HBN and to describe irregular shapes of geographical regions [30].

Figure 1 .
Figure 1.(a) Anji county in northwest Zhejiang province and land use types; (b) Shanchuan town in southeast Anji county and Moso bamboo forest distribution; (c) MODIS LAI product for the Shanchuan town at DOY 225 in 2015 and spatial location of MBF flux measurement site and observed area.

Figure 1 .
Figure 1.(a) Anji county in northwest Zhejiang province and land use types; (b) Shanchuan town in southeast Anji county and Moso bamboo forest distribution; (c) MODIS LAI product for the Shanchuan town at DOY 225 in 2015 and spatial location of MBF flux measurement site and observed area.The MBF flux measurement site is located in Anji County in Zhejiang Province (30.46 • N, 119.66 • E).This area has a subtropical monsoon climate with distinct seasons and abundant rainfall.The MBF is distributed widely across Anji County covering an area of 55,367 ha, which accounts for approximately 45% of the total forested area.The area (1 × 1 km) around the flux tower consists of

Figure 3 .
Figure 3. (a) Image of MBF using digital camera with fisheye lens; and (b) analysis of corresponding LAI using LAI (2000G)-LogCI algorithm in the WinSCANOPY v2009a program.

Figure 4 .
Figure 4. Plot design of observed LAI in observed area of 1 × 1 km (Stars represent observation center, and dots indicate observation points).

Figure 3 .
Figure 3. (a) Image of MBF using digital camera with fisheye lens; and (b) analysis of corresponding LAI using LAI (2000G)-LogCI algorithm in the WinSCANOPY v2009a program.

Figure 3 .
Figure 3. (a) Image of MBF using digital camera with fisheye lens; and (b) analysis of corresponding LAI using LAI (2000G)-LogCI algorithm in the WinSCANOPY v2009a program.

Figure 4 .
Figure 4. Plot design of observed LAI in observed area of 1 × 1 km (Stars represent observation center, and dots indicate observation points).

Figure 4 .
Figure 4. Plot design of observed LAI in observed area of 1 × 1 km (Stars represent observation center, and dots indicate observation points).

21 Figure 5 .
Figure 5. Flow chart of process for estimation of LAI using the hierarchical Bayesian network (HBN) algorithm with multiresolution data (LACC respects locally adjusted cubic-spline capping for processing of MODIS LAI; SG respects Savitzky-Golay smooth for processing of MODIS reflectance; MCMC means Markov Chain Monte Carlo sampling method; HBN_LAI_1000, HBN_LAI_500, and HBN_LAI_250 indicate the assimilated LAI by using HBN algorithm at 1000 m, 500m, and 250 m resolution, respectively; and Field_LAI_1000, Field_LAI_500, and Field_LAI_250 indicate observed LAI at 1000 m, 500 m, and 250 m resolution, respectively.).

Figure 5 .
Figure 5. Flow chart of process for estimation of LAI using the hierarchical Bayesian network (HBN) algorithm with multiresolution data (LACC respects locally adjusted cubic-spline capping for processing of MODIS LAI; SG respects Savitzky-Golay smooth for processing of MODIS reflectance; MCMC means Markov Chain Monte Carlo sampling method; HBN_LAI_1000, HBN_LAI_500, and HBN_LAI_250 indicate the assimilated LAI by using HBN algorithm at 1000 m, 500 m, and 250 m resolution, respectively; and Field_LAI_1000, Field_LAI_500, and Field_LAI_250 indicate observed LAI at 1000 m, 500 m, and 250 m resolution, respectively).
r) ; m r is the number of child nodes at scale r; and v ch(i,r) represents the error of the observed data that obeys a Gaussian distribution v ch(i,r) ∼ N 0, σ 2 V ch(i,r) .In this study, at the first scale, Y represents the simulated LAI from the MODIS LAI (MOD15A2) and LAI dynamic model, and Z represents the 1000-m-resolution LAI data with error perturbation.At the third scale, Y represents the LAI with 250-m resolution transmitted from the first scale; the LAI and reflectance data are linked by the PROSAIL model represented by F ch(i,r) , and the simulated reflectance data Z combined with the 250-m-resolution MODIS reflectance data (MOD09Q1) participate in the HBN.

Figure 7 . 4 Figure 7 .
Figure 7. Time series of Field_LAI and HBN_LAI_500 for four pixels of MBF at 500-m resolution in 2015; and comparison between Field_LAI and HBN_LAI_500 of four pixels (N = 11).

Figure 8 .
Figure 8.Time series of Field_LAI and HBN_LAI_250 for sixteen pixels of MBF at 250-m resolution in 2015.

FigureFigure 9 .
Figure series of Field_LAI and HBN_LAI_250 for sixteen pixels of MBF at 250-m resolution in 2015.Remote Sens. 2019, 11, x PEER REVIEW of 21

Figure 10 .
Figure 10.Comparison of MODIS reflectance and canopy reflectance simulated by the PROSAIL model for MBF (The boxplot represents the statistical characteristics of MODIS reflectance; the dot indicates the average of a certain pixel and the hollow circle indicates outlier for the entire year.The red line represents the average reflectance simulated by the PROSAIL model and the region colored light red is the interval of simulated reflectance).

Figure 10 . 21 Figure 11 .
Figure 10.Comparison of MODIS reflectance and canopy reflectance simulated by the PROSAIL model for MBF (The boxplot represents the statistical characteristics of MODIS reflectance; the dot indicates the average of a certain pixel and the hollow circle indicates outlier for the entire year.The red line represents the average reflectance simulated by the PROSAIL model and the region colored light red is the interval of simulated reflectance).Remote Sens. 2019, 11, x FOR PEER REVIEW 17 of 21

Figure 11 .
Figure 11.Mean value and uncertainty of w ch(i,r) of the process model for MBF.

Table 1 .
Observed LAI of MBF in 2015.

Table 2 .
The parameters of PROSAIL model for MBF.