Estimating Downward Shortwave Solar Radiation on Clear-Sky Days in Heterogeneous Surface Using LM-BP Neural Network

: Downward surface shortwave radiation (DSSR) plays an important role in the energy balance of the earth’s surface. Accurate estimate of DSSR is of great signiﬁcance for the rational and effective use of solar energy. Some parameterization schemes were proposed to estimate DSSR using meteorological measurements given ground-based radiation observation sites are scare and uneven. With the development of remote sensing technique, remotely sensed data can be applied to obtain continuous DSSR in space. Commonly, the spatial resolution of most radiation products is relatively low and cannot meet the needs of certain ﬁelds. Moreover, some retrieval algorithms based on the radiation transfer models are complicated for non-professionals. In this study, a back-propagation (BP) neural network method with Levenberg–Marquardt (LM) algorithm (hereafter referred to as LM-BP) was applied to predict DSSR by building the relationship between measured DSSR and high-resolution remote sensing data from the Advanced Space-borne Thermal Emission Reﬂectance Radiometer (ASTER). The DSSR observations from the four-component radiation sensor installed at the land covered by vegetable, village, maize, orchard, Gobi, sandy desert, desert steppe, and wetland were used to validate the model estimates. The results showed that the estimates of DSSR from LM-BP agreed well with the site measurements, with the root mean square error (RMSE) and the mean bias error (MBE) values of 27.34 W/m 2 and − 1.59 W/m 2 , respectively. This indicates that by combining the LM-BP network model and ASTER images can obtain precise DSSR in heterogenous surface. The DSSR results of this study can provide accurate high-spatial resolution input data for hydrological, evapotranspiration, and crop models.


Introduction
Although the solar radiation energy reaching the ground is extremely less than the total radiation emission energy, it is the main source of the earth's light and heat energy and the main driving force of atmospheric movement [1]. Quantitative accurately downward surface shortwave radiation (DSSR) flux is very important to solar energy applications, climate change research, the surface energy balance, and biological activities [2,3]. The Heihe River Basin (HRB) located in arid regions of northwestern China is one of the most vulnerable areas with climate change and frequent natural disasters [4]. DSSR, as one of the important parameters of climate change and surface energy balance, is very critical to the formulation of the sustainable development strategy for the river basin.
Several methods have been developed to calculate regional DSSR. The algorithms can be divided into the empirical models and the theoretical models. Empirical models generally rely on the linear/non-linear relationship between remote sensing data and

Study Area
The Heihe River Basin is the second-largest inland river basin in northwestern China. It is located in the middle of the Hexi Corridor, roughly between 98°-101°30′ E and 38°-42° N. The climate of the basin is dry, the precipitation is sparse in space and concentrated in time, and the sunshine is sufficient. There is a large temperature difference between day and night. Zhangye Oasis belongs to the middle reaches of the Heihe River Basin (Figure 1), where drought, frost, hot and dry wind, and dust storms are the main disasters that limit the economic and social development of the area. The Heihe Water Allied Telemetry Experimental Research project (HiWATER) experiment was carried out in this area in 2012, including two observation matrices [16]. The large observation matrix is 30 km × 30 km (Figure 1a); the small observation matrix is 5.5 km × 5.5 km (Figure 1b), and the small observation matrix is nested within the large observation matrix. The landcover types of Zhangye Oasis are complex, mainly including maize, vegetables, villages, orchards, deserts steppe, Gobi, sandy deserts, and wetlands. The observation points on each type of underlying surface were set during the HiWATER experiment. Sites 2, 3, 5-16 are set in fields of maize, and the remaining sites were set in the land covered by vegetable (No.1), residential (No.4), orchard (No.17), Gobi (No.18), sandy desert (No. 19), desert steppe (No. 20), and wetland (No.21). The observation items mainly include wind speed/direction, air temperature, humidity, air pressure, four-radiation components, soil temperature and humidity, and surface temperature etc. During the intensive observation period, aerial remote sensing experiments were launched, and a series of ground-based observations such as vegetation coverage, surface albedo, and surface temperature, as well as calculations of relevant eco-hydrological parameters, were carried out.

Remote Sensing Data
The ASTER images were acquired during the HiWATER Project on nine clear-sky days in 2012 (15 June, 24 June, 10 July, 2 August, 11 August, 18 August, 27 August, 3 September, and 12 September) [16]. The surface reflectance of Visible/Near-infrared (VNIR) spectral bands and thermal infrared (TIR) bands were derived from ASTER data by performing radiation correction, atmospheric correction, and geometric rectification. The spatial resolution of VNIR spectral bands is 15-m, and the spatial resolution of TIR bands is 90-m. Land surface temperature (LST) is calculated from the ASTER thermal infrared (TIR) data by referring to Li et al.' research [17]. The land surface albedo (Albedo) is required by shortwave infrared bands and VNIR [18]. Unfortunately, shortwave infrared sensor malfunctioned, and the data was unavailable since April 2008. Therefore, the Albedo is calculated by using the effective three VNIR bands (including green band, red band, and near-infrared band) from ASTER images, as follows: where α 1 , α 2 , and α 3 are the surface reflectance values of the green, red, and near-infrared bands, respectively. The Normalized Difference Vegetation Index (NDVI) was calculated utilizing the relationship between the near-infrared band (ρ NIR ) and the red band (ρ R ) as in Equation (2).
Air temperature (T a ) was provided by Feng et al. [19]. The ratio of the air pressure (press) is the ratio of the local air pressure to the one at sea level, as shown in Equation (3) [12]. DEM comes from http://www.gscloud.cn/.
where z is the altitude, unit m.

Observations
The meteorological observation data in the same period mainly includes air temperature, barometric pressure, surface temperature, and four radiation components, as well as the location information of each station (longitude, latitude, altitude). The detailed information of each site is shown in Table 1. The raw observation data measured at each station is sampled at 10 min. The selected meteorological data has been processed, and the abnormal values have been removed. A simple linear regression method was applied to interpolate the missing values. Additionally, the ground observation data (NDVI, Albedo) corresponding to remote sensing was also manually measured. NDVI was calculated by using the leaf area index (LAI) observations. LAI was measured via Model LAI-2000 plant canopy analyzer. Albedo was measured by FieldSpec Handheld from ASD company. An intercomparison between the Albedo, Press, NDVI, LST, and T a values obtained from the stations with the corresponding remote sensing data was done, the result showed that the root mean square errors (RMSE) of Albedo, Press, NDVI, LST, and T a between the observation values and the satellite values are 0.05, 0.003, 0.002, 1.24 K, and 0.13 K, respectively. The variations of Albedo, NDVI, T a , LST, and Press from June to September in 2012 at each station for different land cover types are shown in Figure 2. The Albedo is generally higher in the land covered by the Gobi, desert steppe, and sandy desert; whilst lower values of Albedo are observed in the maize, orchard, and wetland sites (Figure 2a). The mean Press at wetland station is generally higher than other stations from June to September, and the lowest Press is observed at desert steppe station (Figure 2b). The highest NDVI is observed in maize site (maximum value is greater than 0.7), and the lowest NDVI values are observed at desert steppe (lower than 0.1 for mean value) (Figure 2d). T a and LST are higher in summer and lower in autumn. The highest LST is observed at the sandy desert site, which is generally higher than 320 K throughout the growing season; while the lowest LST at the orchard site is about 300 K (Figure 2e). The T a at each site in the growing season generally shows an increasing trend, followed by a gradual decrease ( Figure 2f). The mean T a at the maize site is 299.07 K and 289.09 K in July and September, respectively. The DSSR shows a decreasing trend for the period of June to September in 2012 ( Figure 2c); there are also obvious differences at different stations for each variable.

LM-BP Neural Network
Artificial neural network (ANN) is an information processing technology similar to the human nervous system. It is a mathematical model that performs distributed parallel information processing [20]. ANN mainly processes the input data and changes its own structure by adjusting the weights between neurons based on the complexity of the system [21]. According to the differences in algorithm, there are nearly 40 kinds of ANNs; BP neural network is the most commonly used one. It is a forward network while error back propagation [22]. The LM-BP is the combination of the gradient descent method and the Gauss-Newton method [23], that is, it is easy to calculate fast. Generally expressed by function trainlm. In the forward propagation process, the input mode is processed gradually from the input layer through the hidden layer and turned to the output layer. The state of each layer of neurons only affects the result of the next layer of neurons. If the expected value cannot be obtained at the output layer, it is transferred to backpropagation, and the error signal is returned along the original connection path. By correcting the weight of each neuron, the error signal is minimized. The specific structure of BP is shown in Figure 3.

LM-BP Neural Network
Artificial neural network (ANN) is an information processing technology similar to the human nervous system. It is a mathematical model that performs distributed parallel information processing [20]. ANN mainly processes the input data and changes its own structure by adjusting the weights between neurons based on the complexity of the system [21]. According to the differences in algorithm, there are nearly 40 kinds of ANNs; BP neural network is the most commonly used one. It is a forward network while error back propagation [22]. The LM-BP is the combination of the gradient descent method and the Gauss-Newton method [23], that is, it is easy to calculate fast. Generally expressed by function trainlm. In the forward propagation process, the input mode is processed gradually from the input layer through the hidden layer and turned to the output layer. The state of each layer of neurons only affects the result of the next layer of neurons. If the expected value cannot be obtained at the output layer, it is transferred to backpropagation, and the error signal is returned along the original connection path. By correcting the weight of each neuron, the error signal is minimized. The specific structure of BP is shown in Figure 3.

Implementation of LM-BP Neural Network
When solving problems with LM-BP neural network, a training data set is needed first because the LM-BP network uses supervised learning. The design of the LM-BP network mainly includes several problems such as the network layers including input layer and hidden layer, nodes, the transfer function, the training method, and the setting of training parameters.
(1) Input/output set: There were 546 input/output pairs at 21 sites from June to September in 2012 used to train the LM-BP network, and then, the relationship between DSSR and the input vector was established. The remaining 189 pairs from June to September were employed to validate the performance of the LM-BP neural network. Observations corresponding to remote sensing data were used as inputs for training the LM-BP network to establish their mathematical relationship with measured DSSR in this study. Input vectors include Albedo, Ta, LST, NDVI, and Press.
(2) Transfer function: The transfer function of the BP network must be differentiable. BP model generally adopts Sigmoid function or Linear function as a transfer function. According to whether the output contains negative values, the Sigmoid function can be divided into the Log-Sigmoid function and the Tan-Sigmoid function. Here, Tan-Sigmoid (tansig) was applied as the transfer function the input layer, and the Linear function (purelin) was a transfer function in the output layer.
(3) Initial weight: The BP network uses an iterative method to determine the weight, so an initial value is required. If the initial value is too large or too small, the model performance will be affected. Usually, the initial value is set between (−2.4/F, 2.4/F) or (−3/√F, 3/√F), and F is the number of neurons connected to the weight input terminal.
(4) The number of network layers: This mainly refers to the design of the hidden layer. Generally, multiple hidden layers can bring better performance, but it will lead to too long training time. In this study, a single hidden layer network was applied because many studies have shown it can process non-linear relationship well. The number of neurons in the single hidden layer (M) can be calculated by (4) where a and b are the number of neurons in the output layer and input layer, respectively, and the value of c is a constant between [0, 10].

Implementation of LM-BP Neural Network
When solving problems with LM-BP neural network, a training data set is needed first because the LM-BP network uses supervised learning. The design of the LM-BP network mainly includes several problems such as the network layers including input layer and hidden layer, nodes, the transfer function, the training method, and the setting of training parameters.
(1) Input/output set: There were 546 input/output pairs at 21 sites from June to September in 2012 used to train the LM-BP network, and then, the relationship between DSSR and the input vector was established. The remaining 189 pairs from June to September were employed to validate the performance of the LM-BP neural network. Observations corresponding to remote sensing data were used as inputs for training the LM-BP network to establish their mathematical relationship with measured DSSR in this study. Input vectors include Albedo, T a , LST, NDVI, and Press.
(2) Transfer function: The transfer function of the BP network must be differentiable. BP model generally adopts Sigmoid function or Linear function as a transfer function. According to whether the output contains negative values, the Sigmoid function can be divided into the Log-Sigmoid function and the Tan-Sigmoid function. Here, Tan-Sigmoid (tansig) was applied as the transfer function the input layer, and the Linear function (purelin) was a transfer function in the output layer.
(3) Initial weight: The BP network uses an iterative method to determine the weight, so an initial value is required. If the initial value is too large or too small, the model performance will be affected. Usually, the initial value is set between (−2.4/F, 2.4/F) or (−3/ √ F, 3/ √ F), and F is the number of neurons connected to the weight input terminal.
(4) The number of network layers: This mainly refers to the design of the hidden layer. Generally, multiple hidden layers can bring better performance, but it will lead to too long training time. In this study, a single hidden layer network was applied because many studies have shown it can process non-linear relationship well. The number of neurons in the single hidden layer (M) can be calculated by (4) where a and b are the number of neurons in the output layer and input layer, respectively, and the value of c is a constant between [0, 10].

Results and Discussion
The structure and training results of LM-BP are shown in Table 2. The performance of the trained LM-BP model is satisfactory, with the mean squared error (MSE) values ranging from 0.13 to 0.53 W/m 2 . Albedo, T a , LST, NDVI, and Press are used to estimate DSSR. For the model training, the observations (Albedo, T a , LST, NDVI, and Press) and DSSR at 21 sites are input to LM-BP model. For the model application, 15-m resolution DSSR on clear-sky days in 2012 was generated by utilizing the trained model combined with regional explanatory variables (Albedo, T a , LST, NDVI, and Press).  Figure 4 shows the results of DSSR at ASTER overpass time predicted by LM-BP network model. The 15-m spatial resolution maps thoroughly exhibit its spatial patterns of DSSR over Zhangye Oasis. There were distinct differences in DSSR under the relatively heterogeneous surface. DSSR values in sandy steppe, Gobi, and sandy desert are much larger than the others since Gobi and sandy desert are dry, and clear sky often occurs due to lack of water vapor, and the elevation in sandy steppe is high there, and the optical air mass is rather high. Furthermore, the DSSR estimates showed significant seasonal variations, where DSSR values have the maximum at the end of June and then declined starting in mid-August when the sunshine hours become shorter and the solar zenith angle become greater. Figure 5 shows the scatter plots of estimated and measured DSSR in different landcover types. Significant differences in estimation results were found for heterogeneous surfaces. Overall, the result indicates that the LM-BP network has a better estimation ability. In Figure 5, the estimated value at the maize site (Figure 5c) was closest to the observations; the estimated value at the desert steppe site (Figure 5g), Gobi site (Figure 5e), and sandy desert site (Figure 5f) have relatively poor agreement with the observations. The main reason is that the land-cover type around the corn site is relatively uniform, while the vegetation coverage of the desert steppe site, Gobi site, sandy desert site is sparse and belongs to the mixed surface. Besides, the estimates of the orchard differ greatly from the observations.   Figure 5 shows the scatter plots of estimated and measured DSSR in different landcover types. Significant differences in estimation results were found for heterogeneous surfaces. Overall, the result indicates that the LM-BP network has a better estimation ability. In Figure 5, the estimated value at the maize site (Figure 5c) was closest to the observations; the estimated value at the desert steppe site (Figure 5g), Gobi site (Figure 5e), and sandy desert site (Figure 5f) have relatively poor agreement with the observations. The main reason is that the land-cover type around the corn site is relatively uniform, while the vegetation coverage of the desert steppe site, Gobi site, sandy desert site is sparse and belongs to the mixed surface. Besides, the estimates of the orchard differ greatly from the observations. The statistical difference comparing the DSSR based on LM-BP network versus the observations for these eight sites are list in Table 3. DSSR modeled by LM-BP was prone to have a good agreement with observations from a four-comparison radiometer (R = 0.85). A larger discrepancy between the estimated DSSR and measured DSSR occurred at the orchard site (RMSE = 46.10 W/m 2 ), Gobi site (RMSE = 38.02 W/m 2 ), and desert steppe site (RMSE = 30.34 W/m 2 ) as compared to the other six observation sites. This is because of mixed image pixels including fruit trees, woodland, and grassland in the orchard site according to high-resolution Google Earth images. Gobi Desert site and desert steppe site are also mixed surface types, including bare soil and vegetation, and the other contributing factor might be the uncertainties in the DSSR observations. As a result, mixed pixels cause large errors. In summary, it is practical to utilize LM-BP model to estimate DSSR on heterogeneous surfaces.  Zhang et al. [24] proposed an integrated approach that incorporates a spectral radiation scheme and an atmospheric broadband transmittance model to estimate DSSR. This integrated method can get a good estimate of DSSR, but too many parameters were required and atmospheric remote sensing products are also required. Besides, the spatial resolutions of the atmospheric parameter products are coarse, for example, MODIS aerosol optical thickness (AOT) at a spatial resolution of 10 km. Thus, the trained LM-BP network model shows a great advantage in retrieving DSSR with high-spatial resolution.
The ANN model with multi-layer network topology makes the physical mechanism of solar radiation estimation become a "black box", so which input parameters are appropriate? Undoubtedly, there must be a close physical relationship between the input parameters and the solar radiation. Xu [25] conducted a research on the influence factors of solar radiation in China. The results showed factors have different effects on solar radiation in different times and regions. In this study, we selected Albedo, NDVI, LST, T a, and Press as input parameters by referring to the input parameters of the empirical model and previous studies [1,6,26]. A sensitivity analysis of the input parameters was carried out. The sensitivity coefficient was obtained based on the results of the sensitivity analysis. Figure 6 shows the inter-comparison of these five variables' weights. Each variable has different effects on DSSR for each site. For vegetable site, albedo has the strongest correlation with DSSR; Press has the weakest correlation with DSSR. For village, maize, desert steppe, and wetland site, LST has a greater influence on DSSR. For the orchard and Gobi site, a high relationship between Press and DSSR is found. NDVI has the strongest correlation with DSSR at the sandy desert site. According to the sensitivity analysis, the five variables (Albedo, LST, T a , NDVI, and Press) selected in this study are closely related to DSSR in Zhangye Oasis. The deviations between the modeled DSSR and the observations measured using four-component radiation are to some extent involved the uncertainties in input variables of Albedo, LST, Ta, NDVI, and Press. Remote sensing pixel scaled data contain noise [27,28]. These uncertainty inputs affect the accuracy of the LM-BP network. Besides, the scale mismatch between the pixel resolution (15-m spatial resolution) of LM-BP outputs and the four-component radiation observations can cause the discrepancies. Another cause of bias is the representativeness of the sample, particularly in areas with sparse vegetation coverage. In these regions, the representativeness of the sample is relatively poor, which affects the final result of the model.
The advantage of this research is that high spatial resolution DSSR can be obtained based on artificial intelligence. Previously, artificial intelligence to simulate solar radiation was concentrated on a point scale due to a lack of regional-scale data. Remote sensing technology provides a continuous signal from space that offers us an opportunity to obtain regional-scale solar radiation. Additionally, LM-BP network model is a simple and efficient method as compared to the simple empirical methods. The performance of the multiple linear regression (MLR) model and LM-BP neural network approach was compared in Feng et al.'s research [11]. As shown in Figure 7, the ANN-based method can obtain solar radiation with high accuracy (R = 0.96, RMSE = 1.34 MJ·m −2 , MBE = 0.15 MJ·m −2 ), and good stability. The deviations between the modeled DSSR and the observations measured using four-component radiation are to some extent involved the uncertainties in input variables of Albedo, LST, T a , NDVI, and Press. Remote sensing pixel scaled data contain noise [27,28]. These uncertainty inputs affect the accuracy of the LM-BP network. Besides, the scale mismatch between the pixel resolution (15-m spatial resolution) of LM-BP outputs and the four-component radiation observations can cause the discrepancies. Another cause of bias is the representativeness of the sample, particularly in areas with sparse vegetation coverage. In these regions, the representativeness of the sample is relatively poor, which affects the final result of the model.
The advantage of this research is that high spatial resolution DSSR can be obtained based on artificial intelligence. Previously, artificial intelligence to simulate solar radiation was concentrated on a point scale due to a lack of regional-scale data. Remote sensing technology provides a continuous signal from space that offers us an opportunity to obtain regional-scale solar radiation. Additionally, LM-BP network model is a simple and efficient  [11]. As shown in Figure 7, the ANN-based method can obtain solar radiation with high accuracy (R = 0.96, RMSE = 1.34 MJ·m −2 , MBE = 0.15 MJ·m −2 ), and good stability. The scope of this study was limited in terms of temporal scale. Solar photovoltaic power generation requires long-term solar radiation data with daily scale and even higher time resolution as reference data, while this paper estimated nine DSSR images on clearsky day at satellite overpass time. In the past ten years, especially geostationary satellite technology, which provides long-term ground-level observations at the minute level, for instance, Himawari-8 series [29], can get a scene image every ten minutes; FY-4 series [30] can get a scene image every four minutes. In future studies, these remote sensing data could be used to obtain more time-resolved long-term surface solar radiation.

Conclusions
In this work, LM-BP network combined with ASTER images has been developed to estimate the DSSR in Zhangye Oasis. The advantage of this study is the combination of artificial intelligence method with remote sensing images. The performances of the model have been evaluated through sufficient observations representing land cover types.
The LM-BP network model generally provides better accuracy in predicting DSSR, with the RMSE values ranging from 10.74 to 46.10 W/m 2 , and the MBE from −41.44 W/m 2 to 35.40 W/m 2 . The MPE value is between 1.11% and 4.38%, respectively. The most optimal and worst RMSE, MBE, MPE values were found for the maize site and orchard site. The LM-BP network model overall underestimated DSSR values (MBE = −1.59 W/m 2 ) according to station observations made at several surface types in the Zhangye Oasis. DSSR showed a trend of first increasing and then decreasing in the growing season.
The LM-BP network model has the ability to obtain DSSR at the heterogeneous surface in Zhangye Oasis. The DSSR results of this study are of vital importance for the local surface energy budget and effective use of solar energy resources. Certainly, the combination of artificial intelligence and remote sensing images should be further applied and tested at other regions in different climatic conditions for extending its applicability.  The scope of this study was limited in terms of temporal scale. Solar photovoltaic power generation requires long-term solar radiation data with daily scale and even higher time resolution as reference data, while this paper estimated nine DSSR images on clearsky day at satellite overpass time. In the past ten years, especially geostationary satellite technology, which provides long-term ground-level observations at the minute level, for instance, Himawari-8 series [29], can get a scene image every ten minutes; FY-4 series [30] can get a scene image every four minutes. In future studies, these remote sensing data could be used to obtain more time-resolved long-term surface solar radiation.

Conclusions
In this work, LM-BP network combined with ASTER images has been developed to estimate the DSSR in Zhangye Oasis. The advantage of this study is the combination of artificial intelligence method with remote sensing images. The performances of the model have been evaluated through sufficient observations representing land cover types.
The LM-BP network model generally provides better accuracy in predicting DSSR, with the RMSE values ranging from 10.74 to 46.10 W/m 2 , and the MBE from −41.44 W/m 2 to 35.40 W/m 2 . The MPE value is between 1.11% and 4.38%, respectively. The most optimal and worst RMSE, MBE, MPE values were found for the maize site and orchard site. The LM-BP network model overall underestimated DSSR values (MBE = −1.59 W/m 2 ) according to station observations made at several surface types in the Zhangye Oasis. DSSR showed a trend of first increasing and then decreasing in the growing season.
The LM-BP network model has the ability to obtain DSSR at the heterogeneous surface in Zhangye Oasis. The DSSR results of this study are of vital importance for the local surface energy budget and effective use of solar energy resources. Certainly, the combination of artificial intelligence and remote sensing images should be further applied and tested at other regions in different climatic conditions for extending its applicability.