An Artiﬁcial Intelligence Approach to Predict Gross Primary Productivity in the Forests of South Korea Using Satellite Remote Sensing Data

: Many process-based models for carbon ﬂux predictions have faced a wide range of uncertainty issues. The complex interactions between the atmosphere and the forest ecosystems can lead to uncertainties in the model result. On the other hand, artiﬁcial intelligence (AI) techniques, which are novel methods to resolve complex and nonlinear problems, have shown a possibility for forest ecological applications. This study is the ﬁrst step to present an objective comparison between multiple AI models for the daily forest gross primary productivity (GPP) prediction using satellite remote sensing data. We built the AI models such as support vector machine (SVM), random forest (RF), artiﬁcial neural network (ANN), and deep neural network (DNN) using in-situ observations from an eddy covariance (EC) ﬂux tower and satellite remote sensing data such as albedo, aerosol, temperature, and vegetation index. We focused on the Gwangneung site from the Korea Regional Flux Network (KoFlux) in South Korea, 2006–2015. As a result, the DNN model outperformed the other three models through an intensive hyperparameter optimization, with the correlation coe ﬃ cient (CC) of 0.93 and the mean absolute error (MAE) of 0.68 g m − 2 d − 1 in a 10-fold blind test. We showed that the DNN model also performed well under conditions of cold waves, heavy rain, and an autumnal heatwave. As future work, a comprehensive comparison with the result of process-based models will be necessary using a more extensive EC database from various forest ecosystems.


Introduction
Forests, which cover 30% of the Earth's land area [1], play an essential role in global carbon flux because of their ability to store much more significant amounts of carbon than other terrestrial ecosystems [2]. Terrestrial gross primary productivity (GPP) in the forest ecosystem is the total carbon generated by the photosynthesis of vegetation. GPP refers to the total uptake of carbon dioxide (CO 2 ) during photosynthesis, which is the primary driver of subsequent ecosystem processes, including vegetation growth and yields [3][4][5]. GPP depends on the change in environmental factors, such as light, temperature, air humidity, air turbulence, and CO 2 concentration. It is also affected by  Figure 2 shows a summary of the data used in this study. Input variables for the four AI models (SVM, RF, ANN, and DNN) include downward shortwave radiation (Rs↓), air temperature (Tair), vapor pressure deficit (VPD), and enhanced vegetation index (EVI) obtained from the Moderateresolution Imaging Spectroradiometer (MODIS) products. GPP is the output value that should be adjusted using the in-situ observations from the EC flux tower.   Figure 2 shows a summary of the data used in this study. Input variables for the four AI models (SVM, RF, ANN, and DNN) include downward shortwave radiation (R s ↓), air temperature (T air ), vapor pressure deficit (VPD), and enhanced vegetation index (EVI) obtained from the Moderate-resolution Imaging Spectroradiometer (MODIS) products. GPP is the output value that should be adjusted using the in-situ observations from the EC flux tower.

Data
2.2.1. Overview Figure 2 shows a summary of the data used in this study. Input variables for the four AI models (SVM, RF, ANN, and DNN) include downward shortwave radiation (Rs↓), air temperature (Tair), vapor pressure deficit (VPD), and enhanced vegetation index (EVI) obtained from the Moderateresolution Imaging Spectroradiometer (MODIS) products. GPP is the output value that should be adjusted using the in-situ observations from the EC flux tower. Figure 2. Summary of the data used in this study. Input variables were the satellite data such as downward shortwave radiation (Rs↓), air temperature (Tair), vapor pressure deficit (VPD), and enhanced vegetation index (EVI). Gross primary productivity (GPP), the output value, was adjusted using the in-situ observations from the eddy covariance (EC) flux tower.

Figure 2.
Summary of the data used in this study. Input variables were the satellite data such as downward shortwave radiation (R s ↓), air temperature (T air ), vapor pressure deficit (VPD), and enhanced vegetation index (EVI). Gross primary productivity (GPP), the output value, was adjusted using the in-situ observations from the eddy covariance (EC) flux tower.

Eddy Covariance Flux Data
For the ground-truth data at the Gwangneung forest, we accumulated the measurements from an open-path EC and the meteorological system at 40 m above ground. In the EC system, vertical wind speed (ω) was observed by a three-dimensional sonic anemometer (CSAT3, Campbell Scientific, Logan, UT, USA), and the CO 2 concentration (ρc) was measured by an open-path infrared gas analyzer (LI-7500, LI-COR Environmental, Lincoln, NE, USA). These data were collected by a data logger (CR3000, Campbell Scientific, Logan, UT, USA) at the frequency of 10 Hz. Then, the half-hourly vertical CO 2 flux, which is equal to NEE, was calculated from the covariance between the fluctuations in ω and ρc. To derive the NEE from the 10 Hz raw data, we carried out the preprocessing such as air density correction and spike detection. Further, quality control, gap-filling, and nighttime CO 2 flux correction were also conducted according to the standard data-processing protocol of KoFlux [28,29]. The forest GPP was calculated by adding the ecosystem respiration (R eco ) to the NEE. Given that the R eco is mainly governed by the thermal regime, the R eco was estimated using the relationship between the temperature and the CO 2 efflux rate derived from nighttime EC measurements. Such efforts should have improved the quality of the GPP data. In addition, meteorological measurements, such as solar radiation, photosynthetically active radiation, net radiation, air temperature, humidity, and precipitation, were sampled every second, averaged over 30 min, and logged in another data logger (CR3000, Campbell Scientific, Logan, UT, USA). They were used as the ancillary variables for the nighttime correction and gap-filling of NEE, and the estimation of R eco and GPP.

Input Data Processing
MODIS is a satellite sensor that has been developed to monitor Earth's atmosphere, land, and ocean environment with the launch of Terra and Aqua satellites in 1999 and 2002, respectively [30]. The input data for the AI models (R s ↓, T air , VPD, and EVI) were derived from multiple MODIS products (Table 1). Most of the data have a spatial resolution of 1 to 5 km, which is suitable for the area size (2.2 km 2 ) of the forest stand at the Gwangneung site.

MYD07_L2
Air and dew-point temperatures, total ozone, and surface pressure 5 km 1 day R s ↓, T air (4) , and VPD (5) MOD13A2 (Terra) MYD13A2 (Aqua) Enhanced vegetation index 1 km 8 days (2) EVI (6) (1) 16 day composite was assigned to each of the 16 days; (2) 8 day composite was interpolated to daily values; (3) downward shortwave radiation; (4) air temperature; (5) vapor pressure deficit; (6) enhanced vegetation index. R s ↓, a key component for the land surface energy budget, is primarily used as a meteorological input variable for the retrieval of GPP [31]. To calculate R s ↓, we used a 16 day composite of the black-and the white-sky albedos for shortwave radiation extracted from the MCD43B3 product. For the parameterization of R s ↓, we employed the Bird Model [32] that considers both diffuse and direct radiation. The R s ↓ is then expressed as: where I dir is the direct radiation; I di f is the diffuse radiation; a g is the ground albedo; a s is the atmospheric albedo. The 16 day albedo was assigned to each of the 16 days because the forest albedo does not significantly vary during a 16 day period [30,33]. From the daily black-and white-sky albedos, the actual albedo was calculated using the MODIS lookup table for a fraction of the diffuse radiation regarding a given solar zenith angle and a 550 nm aerosol optical depth (AOD) [33]. For the parameterization of AOD, we rearranged the daily AOD images using a statistical gap-filling method [34,35]. Total column precipitable water (TPW), total ozone, and surface pressure were also parameterized according to the Bird Model [33]. The daily R s ↓ calculated by this process showed high accuracy with the correlation coefficient (CC) over 0.9 against the Automated Surface Observing System (ASOS) in South Korea [35]. The air and dew-point temperatures from the lowest atmospheric layer (surface layer) were used in the calculation of VPD by subtracting the actual vapor pressure (e a ) from the saturated vapor pressure (e s ) [36]. The e s is a function of air temperature, and the e a can be derived from the e s and the relative humidity calculated by air and dew-point temperatures.
The EVI can mitigate the saturation problem of the normalized difference vegetation index (NDVI), thus the EVI is commonly used in forest studies because it is more suitable for the representation of vegetation greenness in densely vegetated areas [37][38][39]. The MOD13A2 (Terra) and the MYD13A2 (Aqua) products provide 16 day composites for NDVI and EVI. Because the 16 day periods of the Terra and the Aqua products were crossing each other, we could obtain 8 day composites by combining the Terra and the Aqua products. The 8 day EVI was then converted to daily value using a cubic spline interpolation method because the forest greenness usually shows gradual changes such as the cosine curve during a year.

Support Vector Machine
SVM is a technique for optimal classification, and the support vector regression (SVR) is a regression model for the optimally classified data groups derived from the SVM. For optimal classification, a maximum margin hyperplane (MMH) can be set up by maximizing the margin between data groups using nonlinear kernel functions such as polynomial and Gaussian functions [40]. We used the e1071 library in R and the gaussian radial basis function for the nonlinear MMH. The cost parameter was set to 1 for outlier management, and the gamma parameter was set to 1/n (n is the data dimension) for the definition of a normal distribution.

Random Forest
RF is an ensemble method that utilizes a number of decision trees derived from random samples. If necessary, a bootstrap process is performed for resampling by taking account of the sample distribution. A bagging (bootstrap aggregating) process creates a final solution by averaging the results from the bootstrapped trees [41]. In our experiment, the number of trees was set to 500, and the number of variables used for splitting tree branches was set to n/2 (n is the number of input variables), using the randomForest library in R.

Artificial Neural Network
ANN emulates a biological neural system. It consists of an input layer for explanatory variables, multiple hidden layers for nonlinear computations, and an output layer to produce a result [42,43]. The sets of weight and bias in the neural network are optimized by minimizing a cost function between true labels and estimated values. We carried out the hyperparameter optimization using the nnet library in R, such as 2 for the number of hidden units; 100 for the maximum iteration; Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm for the adaptive learning rate [43].

Deep Neural Network
The classic ANN might stop at a locally optimized state, and the generic machine learning might cause an overfitting problem. DNN is a technique developed to mitigate these problems by an intensive optimization in the deep and thick hidden layers in a neural network. The rectified linear units (ReLU) activation function can avoid the vanishing gradient phenomenon during the minimization of a loss function. The regularization can make the weighting scheme of a neural network less biased to a specific neuron, and the dropout methods can cope with unexpected cases in a newly given dataset through a handicapped training with the random elimination of part of the neurons [44]. We adopted the AdaDelta optimizer, which can adjust the learning rate of a neural network adaptively using a dynamic update mechanism and thus can train a DNN model more efficiently and precisely [45]. The configuration process of our DNN model is summarized in Figure 3 [46]. We used the h2o library in R for the hyperparameter optimization, such as 500-500-500 for the hidden unit; 300 for the epoch; ReLU for the activation function; AdaDelta for the optimizer; 20% for the dropout ratio. The classic ANN might stop at a locally optimized state, and the generic machine learning might cause an overfitting problem. DNN is a technique developed to mitigate these problems by an intensive optimization in the deep and thick hidden layers in a neural network. The rectified linear units (ReLU) activation function can avoid the vanishing gradient phenomenon during the minimization of a loss function. The regularization can make the weighting scheme of a neural network less biased to a specific neuron, and the dropout methods can cope with unexpected cases in a newly given dataset through a handicapped training with the random elimination of part of the neurons [44]. We adopted the AdaDelta optimizer, which can adjust the learning rate of a neural network adaptively using a dynamic update mechanism and thus can train a DNN model more efficiently and precisely [45]. The configuration process of our DNN model is summarized in Figure  3 [46]. We used the h2o library in R for the hyperparameter optimization, such as 500-500-500 for the hidden unit; 300 for the epoch; ReLU for the activation function; AdaDelta for the optimizer; 20% for the dropout ratio. which was adopted from [46]. L1-regularization uses the Manhattan norm; L2-regularization uses the Euclidean norm, for the adjustment of weights.

Validations
The daily GPP estimates by the four AI models were validated using the leave-one-year-out (a.k.a. jackknife) method. The calibration was carried out from a 10 year dataset except for one year, and the validation was conducted for the excluded one year. Such cross-validation processes were repeated ten times for each year from 2006 to 2015. Validation statistics such as mean bias error (MBE), mean absolute error (MAE), root mean square error (RMSE), and CC were used.
To compare the performance of the AI predictions with that of the currently operative GPP products, we collected the MODIS and the FluxCom GPP during 2006-2015. MODIS Terra and Aqua products (MOD17A2HGF and MYD17A2HGF) comprise a gap-filled dataset that has less uncertainty than the previous version (MOD17A2 and MYD17A2). The gap-filling process was conducted at the end of each year using a reliable dataset for the fraction of photosynthetically active radiation (FPAR) and LAI [47]. Because these products have a spatial resolution of 500 m and a temporal resolution of 8 days, our daily GPP estimates were aggregated into an 8 day composite for the comparison of performance. For another comparison, we obtained the AI-based FluxCom GPP products with a spatial resolution of 10 km and a temporal resolution of 8 days. However, all the pixels on the Gwangneung site had missing values for the entire period, thus the performance comparison with the FluxCom GPP could not be carried out. which was adopted from [46]. L1-regularization uses the Manhattan norm; L2-regularization uses the Euclidean norm, for the adjustment of weights.

Validations
The daily GPP estimates by the four AI models were validated using the leave-one-year-out (a.k.a. jackknife) method. The calibration was carried out from a 10 year dataset except for one year, and the validation was conducted for the excluded one year. Such cross-validation processes were repeated ten times for each year from 2006 to 2015. Validation statistics such as mean bias error (MBE), mean absolute error (MAE), root mean square error (RMSE), and CC were used.
To compare the performance of the AI predictions with that of the currently operative GPP products, we collected the MODIS and the FluxCom GPP during 2006-2015. MODIS Terra and Aqua products (MOD17A2HGF and MYD17A2HGF) comprise a gap-filled dataset that has less uncertainty than the previous version (MOD17A2 and MYD17A2). The gap-filling process was conducted at the end of each year using a reliable dataset for the fraction of photosynthetically active radiation (FPAR) and LAI [47]. Because these products have a spatial resolution of 500 m and a temporal resolution of 8 days, our daily GPP estimates were aggregated into an 8 day composite for the comparison of performance. For another comparison, we obtained the AI-based FluxCom GPP products with a spatial resolution of 10 km and a temporal resolution of 8 days. However, all the pixels on the Gwangneung site had missing values for the entire period, thus the performance comparison with the FluxCom GPP could not be carried out.  Tables 2 and 3 show the statistics of the daily GPP at the Gwangneung site between 2006 and 2015, including the predictions from SVM, RF, ANN, and DNN as well as the observations from the EC. The mean values of the predicted GPP by the AI models (3.26-3.35 g m −2 d −1 ) were very similar to those of the EC measurements (3.30 g m −2 d −1 ). However, the AI models did not well represent very low or very high GPP values, thus the daily GPP under approximately 0.5 g m −2 d −1 or over about 8.5 g m −2 d −1 was rarely simulated by the AI models. Among the four AI models, the DNN showed the closest predictions in terms of the minimum, maximum, and mean values. Figures 4-7 show the seasonal changes of the predicted daily GPP at the Gwangneung site for the period between 2006 and 2015. The four AI models closely traced the growing season (increasing phase) and the senescence (decreasing phase) of the forest every year. The curve in spring showed an exponential increment, and that of autumn had a relatively longer decreasing period in both AI and EC GPPs. The DNN model followed the GPP peak in May and June better than SVM, RF, and ANN models (Figures 4a, 5a, 6a and 7a). Table 4 shows
Forests 2020, 11, 1000 9 of 17 Figures 4-7 show the seasonal changes of the predicted daily GPP at the Gwangneung site for the period between 2006 and 2015. The four AI models closely traced the growing season (increasing phase) and the senescence (decreasing phase) of the forest every year. The curve in spring showed an exponential increment, and that of autumn had a relatively longer decreasing period in both AI and EC GPPs. The DNN model followed the GPP peak in May and June better than SVM, RF, and ANN models (Figures 4a-7a).           (1) mean bias error; (2) mean absolute error; (3) root mean square error; (4) correlation coefficient.

Seasonal Characteristics of Validation Statistics
We examined the accuracy of daily GPP predictions according to seasons: spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February). Figure 8 shows the CC between the observed and the predicted GPP by the four AI models. The CC of spring and autumn was quite high (around 0.9), while that of summer was relatively low (approximately 0.7 to 0.8). The accuracy of winter was worse than other seasons: the CC for SVM, RF, and ANN was around 0 to 0.3, while the CC for DNN was about 0.45 ( Figure 8). The CC values in the cold months from November to March were also low. Nevertheless, the DNN model showed an overall concentrated scatterplot along the 1:1 line compared to the other three models that had a tendency of underestimation in spring and summer and a tendency of overestimation in autumn (Figures 9 and 10).  Table 4. Validation statistics of daily gross primary productivity (GPP) at the Gwangneung site for the period between 2006 and 2015, predicted by support vector machine (SVM), random forest (RF), artificial neural network (ANN), and deep neural network (DNN). (1) mean bias error; (2) mean absolute error; (3) root mean square error; (4) correlation coefficient.

Seasonal Characteristics of Validation Statistics
We examined the accuracy of daily GPP predictions according to seasons: spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February). Figure 8 shows the CC between the observed and the predicted GPP by the four AI models. The CC of spring and autumn was quite high (around 0.9), while that of summer was relatively low (approximately 0.7 to 0.8). The accuracy of winter was worse than other seasons: the CC for SVM, RF, and ANN was around 0 to 0.3, while the CC for DNN was about 0.45 ( Figure 8). The CC values in the cold months from November to March were also low. Nevertheless, the DNN model showed an overall concentrated scatterplot along the 1:1 line compared to the other three models that had a tendency of underestimation in spring and summer and a tendency of overestimation in autumn (Figures 9 and 10).   . Figure 10. Comparisons between spring, summer, and autumn for the gross primary productivity (GPP) observed vs. predicted by artificial intelligence (AI) models.   . Comparisons between artificial intelligence (AI) models for observed vs. predicted gross primary productivity (GPP) by season. Figure 10. Comparisons between spring, summer, and autumn for the gross primary productivity (GPP) observed vs. predicted by artificial intelligence (AI) models.   were quite good because of temporal smoothing effects, but the DNN estimates outperformed the Terra and the Aqua products for the Gwangneung site. In terms of the RMSE, our result was 58% better than that of Terra ((1.278−0.537)/1.278) and 55% better than that of Aqua ((1.195−0.537)/1.195). Moreover, the DNN estimates showed an excellent CC of 0.975.

Discussions
In this study, we examined the suitability and the potential of AI models for the forest GPP predictions in South Korea. Among the four AI models, DNN showed the best performance in terms of validation statistics and the applicability to unexpected weather conditions. We discussed the performance of the AI models and explored the advantages and the limitations of our experiment.

Performance of AI Models
The results showed a tendency of strong agreement between the GPPs predicted by the AI models and by the EC measurements, explaining approximately 90% of the temporal changes in daily GPP. The phenology of forest vegetation is primarily divided into growth and senescence phases. The predicted GPP at the start of the growing season (increasing phase of each year) and the end of the growing season (decreasing phase of each year) consistently followed the EC measurements. Understanding the seasonal development on daily GPP and environmental factors such as R s ↓, T air , VPD, and EVI is critical to the applicability of the nonlinear relationships between the input variables and the GPP predictions. Each phase follows a different phenological development, leading to different relationships between the input variables and the GPP. The AI models consistently demonstrated a reasonable correlation in the estimation of daily GPP (Table 4). Indeed, AI techniques in recent years have been considered a viable option for the retrieval of terrestrial carbon fluxes on a spatially continuous grid. Our experiment also showed the possibility of the AI-based GPP prediction, particularly for the local-scale applications for the countries with complex terrains such as South Korea. The integration of AI techniques with satellite images can be an alternative to studies of the interactions between atmosphere and forest ecosystems. Previous studies showed that the AI applications combined with satellite images were effective for estimation of GPP, above-ground biomass, ecosystem respiration, and latent heat flux [24,25,48,49].
Moreover, a deep learning approach such as DNN in our experiment can be thought of as a better method than the classical AI models such as RF, SVM, and ANN. The DNN presented a superior performance to the other three models in terms of the validation statistics and the representation of seasonal patterns. The overall accuracy decreased in winter because the site is a warm-temperate deciduous forest in which winter photosynthesis is almost inactive. However, the DNN could simulate the winter GPP values much better than the other three AI models. It is because the DNN is a method developed to resolve local minima and overfitting problems of the classical AI models. The DNN model was more optimized in a deep and thick network using the ReLU activation function, the AdaDelta optimizer, and the appropriate outlier management by regularization and dropout techniques. Moreover, we demonstrated that the DNN model performed well despite the multiple extreme weather events such as cold waves, heavy rain, and an autumnal heatwave in 2011 in South Korea. The comparison with MODIS products also showed that our DNN model was locally optimized for the forests in South Korea on a daily basis, while the parameterization of the MODIS GPP was not specifically suitable for South Korea. Such an outcome of the DNN model was partly because of the spatiotemporally appropriate arrangement of the satellite data. The 16 day albedo and the daily EVI interpolated from 8 day data were also suitable for use in the estimation of daily GPP, which is in accordance with previous studies [50,51].

Limitations and Future Work
We made sure that the AI models such as DNN were very effective in resolving the nonlinear issues for GPP prediction. The DNN model provides a useful method to overcome the limitations of linear and process-based models. However, a limited volume of data can lead to degradation of accuracy since it may not sufficiently represent the various aspects of GPP changes [24,43]. A large database from various locations of EC measurements will be required for a more advanced DNN model. A longer time series to cover various weather events will also be necessary to cope with unexpected outlier cases. By using such big data for GPP prediction, we can produce an alternative GPP map from the improved DNN model. Moreover, a spatially continuous retrieval of forest GPP over South Korea can be achieved after the improved DNN model is combined with a higher-quality dataset. For example, the advanced versions of the MODIS products, such as gap-filled snow-free daily albedo (MCD43GF), gap-filled daily AOD with the multi-angle implementation of atmospheric correction (MAIAC) (MCD19A2), and gap-filled daily surface radiation (MCD18A1), are now available. Additionally, the meteorological reanalysis such as Local Data Assimilation and Prediction System (LDAPS) of South Korea can be used for the local-scale land surface variables such as sensible/latent/soil heat flux, temperature, humidity, wind, and precipitation. The use of meteorological reanalysis will also help avoid the problem of the cloud-covered pixels in satellite images.

Conclusions
Many process-based models for carbon flux predictions have faced a wide range of uncertainty issues in the model result. The complex interactions between atmosphere and forest ecosystems are the main reason for the uncertainty. However, the setting of too many initial parameters might lead to disturbing accuracy improvement. On the other hand, artificial intelligence techniques, which are novel methods to resolve complex and nonlinear problems, have shown a possibility for forest ecological applications. This study is the first step to present an objective comparison between multiple AI models for the daily forest GPP prediction using satellite remote sensing data. In particular, the DNN model outperformed the other three models through an intensive hyperparameter optimization in a deep learning network. We showed that the DNN model also performed well under conditions of cold waves, heavy rain, and an autumnal heatwave. In the experiment, our DNN model showed better accuracy than the operative MODIS GPP product, presumably because of the localized optimization for the forests in South Korea. A spatially continuous retrieval of forest GPP can be achieved after the integration of more improved DNN models and higher-quality datasets from satellite remote sensing and meteorological reanalysis. As future work, a comprehensive comparison with the result of process-based models will be necessary using a more extensive EC database from various forest ecosystems. Such comparisons can produce a synergy between the process-based models that have advantages for the understanding of the mechanism and the AI models that have the ability to solve nonlinear problems. In this sense, more efforts for the AI-based modeling of forest GPP will also be required.