Combining Phenological Camera Photos and MODIS Reflectance Data to Predict GPP Daily Dynamics for Alpine Meadows on the Tibetan Plateau

Gross primary production (GPP) is the overall photosynthetic fixation of carbon per unit space and time. Due to uncertainties resulting from clouds, snow, aerosol, and topography, it is a challenging task to accurately estimate daily GPP. Daily digital photos from a phenological camera record vegetation daily greenness dynamics with little cloud or aerosol disturbance. It can be fused with satellite remote sensing data to improve daily GPP prediction accuracy. In this study, we combine the two types of datasets to improve the estimation accuracy of GPP for alpine meadow on the Tibetan Plateau. To examine the performance of different methods and vegetation indices (VIs), three experiments were designed. First, GPP was estimated with the light use efficiency (LUE) model with the green chromatic coordinate (GCC) from the phenological camera and vegetation index from MODIS, respectively. Second, GPP was estimated with the Backpropagation neural network machine learning algorithm (BNNA) method with GCC from the phenological camera and vegetation index from MODIS, respectively. Finally, GPP was estimated with the BNNA method using GCC and vegetation index as inputs at the same time. Compared with eddy covariance GPP, GPP predicted by the BNNA method with GCC and vegetation indices as inputs at the same time showed the highest accuracy of all the experiments. The results indicated that GCC had a higher accuracy than NDVI and EVI when only one vegetation index data was used in the LUE model or the BNNA method. The R2 of GPP estimated by BNNA and GPP from eddy covariance increased by 0.12 on average, RMSE decreased by 1.13 g C·m−2·day−1 on average, and MAD decreased by 0.87 g C·m−2·day−1 on average compared with GPP estimated by the traditional LUE model and GPP from eddy covariance. This study puts forth a new way to improve the estimation accuracy of GPP on the Tibetan Plateau. With the emergence of a large number of phenological cameras, this method has great potential for use on the Tibetan Plateau, which is heavily affected by clouds and snow.


Introduction
Gross primary production (GPP) refers to the overall photosynthetic fixation of carbon per unit space and time [1]. GPP is extremely important for studying the terrestrial ecosystem carbon cycle and global climate change, because GPP reflects not only the vegetation activity, but also the response of the ecosystem to climate change [2][3][4][5]. Therefore, accurately estimating the temporal and spatial dynamics of GPP on a regional and global scale is of great significance for understanding the terrestrial ecosystem carbon cycle, assessing grain production, predicting future climate change scenarios and making climate-policy decisions [6][7][8]. Although there are many GPP models have been developed in the past few decades, it is still a challenging work to accurately estimate daily GPP.
Monteith first proposed the method of estimating vegetation productivity by using photosynthetic effective radiation absorbed by vegetation and established the basis of light use efficiency model [1]. This method is widely used in regional GPP estimation with remote sensing data [9,10]. However, owing to the influence of atmospheric conditions, cloud and sensor degradation, there were large uncertainties in GPP estimated by remote sensing vegetation index data and the light use efficiency (LUE) model [11], such as the CASA model [12], the global production efficiency model (GLO-PEM) [13], the vegetation photosynthesis model (VPM) [14] or the eddy covariance-light use efficiency (EC-LUE) model [15]. However, it is reported that the satellite remote sensing-based GPP contains many problems caused by snow and cloud contaminants and low temporal resolution. For example, GPP was underestimated at the majority of sites in African ecosystems [16] and has unstable performance at various sites [17,18].
Phenological cameras can provide almost time-continuous images at a small spatial scale; meanwhile, they are rarely disturbed by clouds [19][20][21]. Phenological cameras have the advantages of low cost, small size, convenient installation, and easy handling, and they are widely used to monitor plant dynamics [22]. The phenological camera-based color index was also used to estimate canopy GPP in many studies [23][24][25]. The performance of GPP estimated by color indexes from phenological camera photos was of varying vegetation types, and the estimation effect of GPP on grassland was better than that of evergreen coniferous forests and deciduous broad-leaved forests [11]. Thus, the phenology camera can provide more reliable daily greenness dynamics and can be fused with satellite remote sensing data to improve daily GPP estimate for alpine meadow ecosystems.
In recent years, data-based machine learning has been widely used in remote sensing data fusion and model integration. Many studies have indicated that machine learning works well in vegetation productivity estimation [26][27][28]. Machine learning models find rules by training sample data and use these rules to predict unknown results. GPP estimation based on machine learning algorithms can make full use of existing multisource remote sensing data, and they can provide a long time series and high-precision estimation of GPP. Therefore, machine learning has advantage of combining multisource data by machine learning to estimate GPP. The machine learning method predicted GPP with an RMSE of 1.87 g C·m −2 ·day −1 and an R 2 of 0.71 using meteorological and flux data from 33 AmeriFlux sites between 2000 and 2003 in America [29]. The machine learning method was also used to estimate eight days GPP (at the site level) with superior performance on the European spatial scale [30]. Lee successfully estimated the daily GPP of the forest ecosystem using the machine learning method with MODIS data at site scale, and the RMSE reached 1.03 g C·m −2 ·day −1 [31]. These studies mentioned above indicated that the machine learning methods is an effective and widely used method to correctly predict GPP.
The rapid use of phenology cameras in the Tibetan Plateau provides high-quality data that reflect alpine seasonal changes. However, in what extent it can be used to improve GPP estimation for alpine meadows has not been studied. In view of the advantages and disadvantages of phenology camera digital photos and MODIS data, combing digital camera and satellite remote sensing data is an effective way to improve GPP prediction accuracy on the Tibetan Plateau because digital cameras and satellite remote sensing provide complementary information. Thus, the objective of this study are to: (1) compare the performance of color index from phenological photos and vegetation indices from satellite remote sensing data in GPP estimation for alpine meadows on the Tibetan Plateau; (2) compare the performance of the LUE model and a machine learning method in GPP estimation, and (3) explore to what extent the GPP estimation accuracy can be improved by combining multisource remote sensing data for alpine meadows on the Tibetan Plateau.

Study Area
The Tibetan Plateau is the largest and highest plateau on Earth [32]. The ecosystems of the Tibetan Plateau are sensitive and fragile, and they are sensitive to global climate change. Alpine meadows are a dominant vegetation type on the Tibetan Plateau. The A'rou (AR) eddy covariance (EC) observation station ( Figure 1) is located in the eastern part of the Tibetan Plateau. The land surface of the AR station is alpine meadow. The longitude and latitude of the AR station are 100°27′52.9″ E and 38°02′39.8″ N and the altitude is 3033 m. The mean annual temperature and annual precipitation are 0.7 °C and 400 mm, respectively. We divided the seasons on the Tibetan Plateau, with spring as May, summer as June to July, and autumn as August to September [33,34].

Meteorological Data
The meteorological data from 2017 to 2018 observed at the AR station were used as input to simulate GPP, including global solar radiation (RG, w/m 2 ), saturated vapor pressure deficit (VPD, Pa), air temperature (T_2m, air, °C) at a height of 2 m above ground level, and photosynthetically active radiation (PAR, MJ/d/m 2 ). These variables are observed every 10 min, and daily RG, VPD, PAR, and minimum temperature (Tmin) were calculated from raw 10 min data in this study ( Figure 2).

Remote Sensing Data
MOD09GA is the daily MODIS surface reflectance product, which was atmospherically corrected [35]. The band reflectance data of MOD09GA can be obtained from Google Earth Engine (https://earthengine.google.com/) and its spatial resolution and temporal resolution are 250 m and daily, respectively [36]. The 3 × 3 pixels around the AR station were extracted from MOD09GA data ( Figure 3). Daily NDVI and EVI were calculated from the reflectance data using the following formulas: where ρRED, ρNIR, and ρBLUE are the band reflectance of the red, near-infrared, and blue band, respectively. C is the adjustment parameter, and the value of C is set to 1.

GPP from Eddy Covariance
GPP data from AR eddy covariance flux site ( Figure 1) on the Tibetan Plateau were used to evaluate the GPPs derived from remote sensing data and phenological camera photos. Eddy covariance raw flux data were processed rigidly, including despike, time lag correction, coordinate rotation, frequency response correction, and WPL correction [37]. After gap-filling and partition, daily GPP was calculated based on ecosystem respiration (Re) and net ecosystem exchange (NEE) data by using Equations (3) and (4) [38]: where εmax is the maximum LUE in g C/MJ APAR, GEEmax is the maximum gross ecosystem exchange in g C/m 2 /s [39].

Phenological Camera Photos and Color Index Extraction
Phenological cameras (NV201E digital cameras, Haoxin Netscape Network Technology Co., Ltd, Shenzhen, China) were installed at the AR station on 15th April, 2017. They continually record the surface status by taking photos almost vertically and toward the south to avoid the tower shadow. These phenological cameras took photos with a resolution of 1600 × 1200 and can specify the shooting time and frequency ( Figure 4). During the observation period, photos were saved in JPEG format.
These phenological cameras were operated in automatic exposure mode with the white balance adjustment in "automatic mode". In 2017 year, the times the photos were taken were 7:30, 7:45, 8:00, 11:30, 11:45, 12:00, 17:30, 17:45, and 18:00. In 2018 year, the time was adjusted and the times were 8:30, 8:45, 9:00, 11:30, 11:45, 12:00, 17:30, 17:45, and 18:00 respectively. We deleted the photos when the light was too strong or too weak from the photo set, so as to get the photos of real recorded features under normal weather conditions. The green chromatic coordinate (GCC) was calculated to reflect the canopy greenness, which was widely used to monitor canopy development and was strongly related to GPP [40,41]. GCC was calculated using following formulas: where DN is the digital number and R, G, and B were the pixel values of red, green, and blue channels, respectively. GCC1 and GCC2 are both green chromatic coordinate (GCC), but the denominators in the calculation Equations (5) and (6) were different. To ensure the reliability of GCC data, we set some conditions on the quality of the photos as following: (1) Make sure that there were no underexposed photos. (2) No shadows of the the flux tower appeared in the photos. (3) No obvious animals and people appeared in the photos. We extracted digital RGB data from the ROI (400 × 600) of the images and calculated GCC1 and GCC2 for all photos. We found that the R 2 of ROI_GCC1 and GPP from eddy covariance were higher than the R 2 of ROI_GCC2 and GPP from eddy covariance ( Figure 5), so ROI_GCC1 was used for all GPP estimations in this study.

GPP Model Based on VIs and GCC
The LUE model was used in study to estimate GPP, as follows: where ε is light use efficiency in photosynthesis (gC/MJ) and PAR (MJ/d/m 2 ) is photosynthetically active radiation. PAR = n × RG. n is the proportion of PAR in the global total solar radiation and was not a fixed value. For calculations, n = 0.45 was used. FPAR is the fraction of absorbed photosynthetically active radiation by vegetation. In this study, FPAR was estimated from the normalized color index (CInorm) of GCC, NDVI and EVI by Equations (9)- (11). α is the coefficient constant and is assumed to be 1.0 [14]. The value of εmax was set to 1.6 g C·m −2 ·day −1 ·MJ −1 according to literature for alpine meadow [39]. f(Tmin) and f(VPD) are simple linear piecewise slope functions of daily minimum temperature (Tmin) and daily average VPD [42]. f(Tmin) and f(VPD) reflect the environmental stress on vegetation photosynthesis and are expressed as: Tmin < Tmin Tmin − Tmin Tmin − Tmin , Tmin < Tmin < Tmin 1, Tmin > Tmin where VPDmax, VPDmin, Tminmax, and Tminmin are model parameters and the values were from the default parameters for grass in the BPLUT table [43].
The GCC values of RGB color photos of the AR station from 2017 to 2018 were calculated in Python 3.5, and the CInorm index was calculated with the GCC values (Equation (14)). The CImax and CImin were the maximum and minimum value of GCC in each year, respectively. The GPPs based on CInorm and the remote sensing VIs was conducted by Equations (15)-(17), and the equations are given by:

BP Neural Network Algorithm
The BP neural network is a multilayer feedforward network trained by error back propagation. The BP neural network algorithm (BNNA) contains an input layer, hidden layer, and output layer, and each layer includes a specific number of neurons. We chose BNNA for this study because it had the following four capabilities: 1. BNNA has a strong nonlinear mapping ability; 2. BNNA has a highly self-learning and self-adaptive ability; 3. BNNA has the ability to apply learning results to new knowledge; 4. BNNA has a certain fault tolerance. We used Keras based on TensorFlow high-level Application Programming Interface (Google Artificial Intelligence Team Google Brain, Mountain view, Santa Clara County, CA, USA) and BNNA to achieve machine learning fusion of meteorological data, MOD09GA VIs, and phenological cameras GCC to predict GPP. First, the input datasets were normalized and mapped to the range [0,1], and then BNNA was trained with data in 2017. Training was stopped when BNNA reached the iteration stop condition (Terror = 1 × 10 −7 ) and number of iterations (epochs = 100). The data in 2018 were used to validate the BNNA to test the performance of established model. R 2 , RMSE, Mean Absolute Deviation (MAD) indices between the GPP from eddy covariance and simulated GPP were calculated to evaluate the accuracy during training and validation periods. The structure and key parameters of BNNA are reported in Figure 6 and Table 1. The process of machine learning GPP modeling consisted of four groups of experiments (Table 2).  Too many or too few iterations of the training model in the BNNA will cause uncertainty in training errors. Only repeated experiments can sum up a relatively minimum number of trainings and obtain better training results. Finally, we set n = 100. If the activation function was not set properly, the BP neural network will not converge. We compared the convergence effects of different activation functions and found that the effects of the "tanh" activation was better than others for the training. The learning rate was generally set to 0.1 for BNNA.

Model Output Accuracy Evaluation
We calculated the R 2 , RMSE, and MAD in this study to evaluate the performances of the GPPs (GPP_NDVI, GPP_EVI, GPP_CInorm) estimated by the remote sensing GPP model and the GPPs (GPP_train1, GPP_train2, GPP_train3, GPP_train4, GPP_valid1, GPP_valid2, GPP_valid3, GPP_valid4) estimated by the BNNA (Figure 7). P and O, represent the GPP simulation values and the GPP from eddy covariance, respectively, and i and n represent the days of the year and the total number of days in the year, respectively. The R 2 , RMSE, and MAD are expressed as follows:

Simulating GPP with GCC
The seasonal GPP dynamics estimated from GCC are shown in Figure 8 and Figure A1. The total growing season GPP was 790.48 g C·m −2 ·day −1 in 2017 and 809.24 g C·m −2 ·day −1 in 2018, respectively. Generally, GPP_CInorm showed the same seasonal dynamics with the GPP from eddy covariance in 2017 and 2018. The R 2 , MAD, and RMSE of GPP_CInorm and the GPP from eddy covariance in the growing seasons of 2017 were 0.68, 2.11 g C·m −2 ·day −1 , and 2.47 g C·m −2 ·day −1 , respectively. In the year 2018, R 2 , MAD, and RMSE between GPP_CInorm and GPP from eddy covariance were 0.78, 1.63 g C·m −2 ·day −1 , and 1.94 g C·m −2 ·day −1 , respectively. There was a systematic deviation between GPP from GCC and GPP from EC in spring and autumn, and these phenomena appeared both in 2017 and 2018. In 2017, the GPP_CInorm had the highest accuracy in autumn, followed by summer, and spring was the worst in the LUE model (Table 3).  In 2018, the GPP_CInorm had the highest accuracy in summer, followed by autumn, and spring was the worst in the LUE model (Table 3). More information about the changes of GPP from eddy covariance and GPP_CInorm in 2017 and 2018 can be found in the appendix (Appendix A, Figure A1).

Simulating GPP with MOD09GA VIs
The comparison between GPP from MODIS NDVI and GPP from eddy covariance in 2017 and 2018 at the AR station are shown in Figure 9a,b. In the growing season of 2017, the R 2 , RMSE, and MAD between GPP_NDVI and GPP from eddy covariance were 0.36, 3.30 g C·m −2 ·day −1 , and 2.46 g C·m −2 ·day −1 respectively. In 2017 growing season, the GPP_NDVI had the highest accuracy in autumn, followed by spring and summer ( Table 3). The R 2 , RMSE, and MAD between GPP_NDVI and GPP from eddy covariance were 0.46, 3.29 g C·m −2 ·day −1 , and 2.45 g C·m −2 ·day −1 in the growing season of 2018, respectively. In 2018, the GPP_NDVI had the highest accuracy in autumn, followed by summer and spring (Table 3). The comparison between GPP from MODIS EVI and GPP from eddy covariance in 2017 and 2018 at the AR station are shown in Figure 9c,d. In the 2017 growing season, the R 2 , RMSE, and MAD of GPP_EVI and GPP from eddy covariance were 0.53, 2.84 g C·m −2 ·day −1 , and 2.19 g C·m −2 ·day −1 , respectively. The LUE model with EVI as input (GPP_EVI) had the highest accuracy in spring, followed by autumn and summer (Table 3). During growing season in 2018, the R 2 , RMSE, and MAD of GPP_EVI and GPP from eddy covariance were 0.59, 2.85 g C·m −2 ·day −1 , and 2.22 g C·m −2 ·day −1 , respectively. The LUE model with EVI as input (GPP_EVI) performed best in autumn, followed by summer and spring (Table 3).
GPPs estimated from MODIS NDVI (GPP_NDVI) and EVI (GPP_EVI) were able to successfully capture the season dynamics of GPP from eddy covariance (Figure 9a-d). However, throughout the growing season, daily GPP_NDVI/GPP_EVI was dramatically underestimated at some dates comparing with GPP from eddy covariance, this was mainly caused by cloud or snow contamination of the satellite remote sensing data. The number of cloud contaminated data during the growing season in 2017 was 56, and 49 in 2018. The cumulative MAD caused by GPP_NDVI and GPP_EVI were 3.794 g C·m −2 ·day −1 -3.849 g C·m −2 ·day −1 and 3.151 g C·m −2 ·day −1 -3.193 g C·m −2 ·day −1 , respectively. More information about the changes of GPP from eddy covariance and GPP_NDVI and GPP_EVI in two years (2017 and 2018) can be found in the appendix (Appendix A, Figure A2).

Simulating GPP with multiSource Data Using Machine Learning
In machine learning GPP estimation experiment, data in 2017 were used for training, and data in 2018 were used for validation. Four experiments were designed to test the contribution of different input datasets. The time period of the four GPP training experiments was from 26th April to 30th September, 2017. The time period of the four GPP validation experiments was from 26th April to 30th September, 2018. Figures 10a-d and 11a-  During the growing season, autumn had the highest training accuracy, followed by summer and spring corresponding to GPP_train1. For the GPP_train2, spring had the highest training accuracy, followed by autumn and summer. In the GPP_train3 experiment, spring had the highest training accuracy, followed by summer and autumn. In the GPP_train4 experiment, autumn had the highest training accuracy, followed by summer and spring (Table 4). In the four validation experiments (Figure 11), it is found that autumn had the highest validation accuracy for all of these experiments (GPP_valid1, GPP_valid2, GPP_valid3 and GPP_valid4), followed by summer and spring.

Comparing the Performance of Different GPP Estimation Methods
The statistical indicators of GPPs estimated by MOD09GA and GCC according to the LUE model and BNNA are shown in Table 5. The seasonal patterns between the GPP from eddy covariance and GPP from machine learning method were generally consistent, except for some large fluctuations between May and June. We compared the R 2 , RMSE, and MAD of the training results of GPP and found that the training accuracy of the GPP_train4 was the highest followed by GPP_train1, GPP_train3, and GPP_train2. The training accuracy of the GPP_train4 indicated that the training results can be improved by combining multisource data compared with the GPP_train1, GPP_train2, and GPP_train3 in machine learning. In addition, we found that the validation accuracy of GPP_valid4 was the highest followed by GPP_valid1, GPP_valid3, and GPP_valid2. Therefore, the effect of GPP_train4 and GPP_valid4 was the best of machine learning. Among the LUE models forced by different vegetation indices, the model accuracy from high to low were GPP_CInorm, GPP_EVI, and GPP_NDVI. When compared with the LUE models, the R 2 , MAD, and RMSE of machine learning were clearly better than those of the LUE models. Therefore, the machine learning method was an effective approach and can significantly improve the simulation accuracy of daily GPP as compared to the LUE model.

GPP Estimation Based on GCC
The GCC can be used to quantitatively describe the vegetation canopy greenness in alpine meadows on the Tibetan Plateau, and the CInorm derived from GCC can be used to estimate the vegetation GPP. This was consistent with the finding that the GCC can quantify the total GPP of the grassland and deciduous broad-leaved forest [11]. In this study, we found that GPP_CInorm has the highest accuracy, followed by GPP_EVI, and GPP_NDVI was the worst among them. This consequence illustrates that long-term phenological photos can replace remote sensing data thereby obtaining higher-precision regional GPP. In addition, we found that there was no deviation in the GCC peak value obtained from the AR station, but GCC had systematic deviation in spring and autumn. However, previous studies reported that the GCC peak date in growing season is much different from observed GPP in forests and grasslands [21]. The GCC peak value of the phenology photo of the TKY site (deciduous broad-leaved forest) and the TKC site (evergreen coniferous forest) in Japan appeared about 20-30 days before the GPP peaks [26]. GCC and GPP of the Tongyu (meadow steppe) site in the Songnen Plain of China and the Maodeng Ranch (typical steppe) site in Inner Mongolia showed consistent and noticeable seasonal changes. GCC and GPP had large values during the growing season, and the values during the dormant period were very close to zero. The peak of GCC of meadow steppe lagged behind the peak of GPP by about ten days, and the peak of GCC of typical steppe lagged behind by GPP by about 20 days [44].
The advantages of GPP estimated by GCC include high continuity of phenology photos and that phenology photos can be obtained multiple times a day with less affects from clouds. In addition, GCC will be affected by some factors when estimating GPP, such as the exposure rate of the photo, the spatial scale of the photo, the height angle of the sun, the presence of shadows or animals in the photo, temperature, and photosynthesis. These factors will affect the values of GCC and thus GPP estimation. To our knowledge, there have been few studies on the use of phenological photos (GCC) to estimate GPP on the Tibetan Plateau. Therefore, the performance of GCC in GPP estimation for alpine meadow on the Tibetan Plateau was tested in this study.

GPP Estimation Based on Remote Sensing Vegetation Indices
The remote sensing data had great advantages of wide coverage and high spatiotemporal resolution, which can be used to estimate GPP for regional and even global over a long period of time. However, the complexity of cloud cover, snow, moisture, temperature, and vegetation growth will all affect the accuracy of GPP estimation by remote sensing [45]. We used the remote sensing vegetation index MOD09GA NDVI and EVI to estimate the GPPs of the AR station on an annual scale. For NDVI, we found that R 2 was 0.50 and RMSE was 2.78 g C m −2 day −1 in 2017; R 2 was 0.57 and RMSE was 2.76 g C m −2 day −1 in 2018. For EVI, R 2 was 0.85 and RMSE was 1.65 g C m −2 day −1 in 2017; R 2 was 0.87 and RMSE was 1.56 g C m −2 day −1 in 2018 (Appendix A, Figure A2). Li estimated the GPPs of AR station in 2009 and reported that R 2 of LUE models were up to 0.8 at AR station at annual scale [46]. This is consistent with the performance of the LUE model in this study. This study found that at the AR station, using EVI to estimate GPP was better than using NDVI. This was also reported by a previous study [47]. Compared with NDVI, the blue band was used to correct the background effects when calculating EVI. This difference may be one of the reasons that EVI had a better correlation with GPP than NDVI [48]. These findings were consistent with conclusions that the remote sensing GPP model can accurately estimate GPP in alpine meadows.
The GPP_NDVI and GPP_EVI estimated by the LUE model based on the MOD09GA vegetation indices fluctuated greatly during the growing seasons, which were lower than the GPP from eddy covariance (Figure 9). GPP_NDVI and GPP_EVI had many cloud contaminated data points close to zero during the growing season, which may be caused by some noise points in the MOD09GA product. For GPP_NDVI, the number of cloud-contaminated data in the growing season of 2017 year was 56, R 2 of GPP_NDVI and GPP from eddy covariance was 0.22, the number of without cloud contaminated data was 92, and the R 2 of GPP_NDVI and GPP from eddy covariance was 0.61. The number of cloud contaminated data in the growing season of 2018 year was 49, the R 2 of GPP_NDVI and GPP from eddy covariance was 0.26, the number of without cloud contaminated data was 99, and the R 2 of GPP_NDVI and GPP from eddy covariance was 0.60. For GPP_EVI, the number of cloud contaminated data in the growing season of 2017 year was 48, R 2 of GPP_EVI and GPP from eddy covariance was 0.40, the number of without cloud contaminated data was 100, and the R 2 of GPP_EVI and GPP from eddy covariance was 0.64. The number of cloud contaminated data in the growing season of 2018 year was 46, the R 2 of GPP_EVI and GPP from eddy covariance was 0.43, the number of without cloud contaminated data was 102, and the R 2 of GPP_EVI and GPP from eddy covariance was 0.71. For the non-growth seasons, because all the alpine grasses were withered and photosynthesis was weak, GPP_NDVI and GPP_EVI were very small and close to zero, which was consistent with the GPP from eddy covariance (Appendix A, Figure A2).

The Impact of εmax for GPP Estimation
The maximum light energy utilization rate (εmax) of the remote sensing GPP model is an important parameter and has an important impact on the accuracy of GPP modeling. Turner revealed that the MODIS GPP algorithm may cause εmax to be underestimated [49]. MODIS GPP products in the midwestern U.S. croplands were underestimated due to the inaccuracy of εmax estimation. Therefore, the LUE values measured in the field should be used consistently in large-scale GPP modeling [50]. Lin conducted a monthly parameterization of the one-leaf LUE model and the twoleaf LUE model of eighteen typical plant function types. The results indicated that the εmax of most plant function types change with the season and had a significant positive correlation with the leaf area index and environmental temperature changes [17]. Wei developed a data-driven RFR-GPP model using random forest regression method, which showed LUE varied largely with latitude [51]. Zheng examined different sources of uncertainty in GPP simulated by LUE models and found that model structures (e.g., FPAR, Ws, and Ts) influenced both model parameters and model performance [52]. The εmax in alpine meadows at different altitudes on the Tibetan Plateau estimated to be 1.83 gC/MJ and 1.33 gC/MJ [47,53]. The εmax of this study was equal to 1.6 gC/MJ in alpine meadows [39]. εmax also vary with spatial and even in different year. Using a constant εmax for a certain biome may cause errors in LUE-oriented GPP modeling [54]. Therefore, obtaining reasonable εmax is vital to LUE models.

BNNA Performance Evaluation for GPP Estimation
In recent years, machine learning methods have been widely used in GPP estimation for its flexibility in application. In this work, it was found that machine learning algorithm (BNNA) performed well in alpine meadow GPP prediction. The machine learning algorithm (BNNA) was more accurate than the LUE models. The advantages of using machine learning method to fuse multisource data to estimate GPP were that machine learning methods can easily fuse data from different sources together and can combine the advantages of different data for GPP estimation. However, the main shortages were the lack of mechanism process, the need for a large number of training samples, and the results of training in different regions cannot be applied to other regions. The implication of this study was to provide a new idea to improve the GPP estimation accuracy on the Tibetan Plateau. With the large number of phenological cameras, the method will have good application potential on the Tibetan Plateau, which is greatly affected by clouds and snow.

Conclusions
The GPP_CInorm has the highest accuracy, followed by GPP_EVI, and GPP_NDVI was the worst in the LUE model. However, it should be noted that there is a problem with the seasonal dynamics of CInorm results, manifested as a significant deviation existing in spring and autumn. In several sets of machine learning experiments, the accuracy of GPP was found to be the highest when both the CInorm based on phenological camera photo and MODIS vegetation indices (NDVI and EVI) are input at the same time, indicating that the two types of data can complement each other's defects and improve estimation accuracy of GPP. Compared with the GPP accuracy estimated by EVI and NDVI in the LUE model, the GPP accuracy estimated by combining multisource data is increased by 27% and 59%, respectively. These results illustrate that the overall accuracy of machine learning is significantly higher than the traditional LUE model. However, compared with the GPP accuracy estimated by CInorm in the LUE model, it is reduced by 3%, which is caused by the lower quality of GPP data from eddy covariance and the shorter time series of GCC data. Due to the late installation of phenological camera resulting in the quantity of photos available now is very small. Thus, we can get longer time series GCC data later, the purpose is to reduce the impact of the amount of data on the accuracy of machine learning estimation GPP.