The Ultra-Short-Term Forecasting of Global Horizonal Irradiance Based on Total Sky Images

: Solar photovoltaics (PV) has advanced at an unprecedented rate and the global cumulative installed PV capacity is growing exponentially. However, the ability to forecast PV power remains a key technical challenge due to the variability and uncertainty of solar irradiance resulting from the changes of clouds. Ground-based remote sensing with high temporal and spatial resolution may have potential for solar irradiation forecasting, especially under cloudy conditions. To this end, we established two ultra-short-term forecasting models of global horizonal irradiance (GHI) using Ternary Linear Regression (TLR) and Back Propagation Neural Network (BPN), respectively, based on the observation of a ground-based sky imager (TSI-880, Total Sky Imager) and a radiometer at a PV plant in Dunhuang, China. Sky images taken every 1 min (minute) were processed to determine the distribution of clouds with di ﬀ erent optical depths (thick, thin) for generating a two-dimensional cloud map. To obtain the forecasted cloud map, the Particle Image Velocity (PIV) method was applied to the two consecutive images and the cloud map was advected to the future. Further, di ﬀ erent types of cloud fraction combined with clear sky index derived from the GHI of clear sky conditions were used as the inputs of the two forecasting models. Limited validation on 4 partly cloudy days showed that the average relative root mean square error (rRMSE) of the 4 days ranged from 5% to 36% based on the TLR model and ranged from 12% to 32% based on the BPN model. The forecasting performance of the BPN model was better than the TLR model and the forecasting errors increased with the increase in lead time.


Introduction
Solar photovoltaics (PV) is developing at an unprecedented speed, which is attributed to dramatic cost reductions, technology advancements and government policy support [1][2][3]. According to the IEA (International Renewable Energy Agency) report, global solar PV power generation increased by 22% to 720 TWh in 2019 [4]. However, the ability to forecast PV power remains a key technical challenge for the integration of large-scale PV into the electrical grid due to the variability and uncertainty of solar irradiance, which affects PV power most directly [5][6][7][8][9]. Moreover, accurate forecasting of solar irradiance is crucial to the dispatch and management of power systems [5]. Further ultra-short-term systems for forecasting highly spatially and temporally resolved solar irradiance in a timeframe of 15 min have been put forward to optimize the operation of solar power plants [10][11][12]. In fact, the amount of GHI (global horizontal irradiance) is closely related to atmospheric conditions including cloud, transparency, aerosol concentration and water content [13], in which cloud is the most dominant factor [14,15], data were sampled synchronously with the frequency of 1 min, and simultaneously stored locally and remotely to guarantee the timeliness, integrity and reliability of the data.

Cloud Detection and Retrieve
Feature extraction from cloud images is a key component in the application of the Total Sky Imager [12]. Although TSI-880 provides cloud detection according to the built-in algorithm, the single fixed red-blue-Ratio (RBR) threshold for all images and untreated shielding caused significant error in the results. To be specific, on the one hand, the cloud information near the sun shielded by the shadow band is vital for ultra-short-term irradiance forecasting [34]. On the other hand, the image RBR is greatly affected by atmospheric conditions including transparency, aerosol content and sand dust, even though for the same clear sky image, the RBR values around the sun and near the horizon are larger due to aerosol forward scatter and long optical path [40,41]. Given this, we performed cloud detection based on the RBR threshold to determine the cloud fraction.
The first step of cloud detection is recovering the shielding areas that result from the support arm and shadow band, as shown in Figure 2a. The angle of the shadow band deviates from the true north on the image, which is equal to the solar azimuth (SAA), because projection of the sky on the image is isometric [42], and the area of shadow band is shown in the green area of Figure 2b and interpolated according to the pixels on both sides. The projection of the support arm with a fixed position is shown in the white area of Figure 2b and is restored by horizontal bidirectional interpolation. The restored image is shown in Figure 2c. Additionally, the sky, effectively occupying 260 × 260 pixels, is cropped from the original image and distortion is corrected by the spherical distortion correction model [43]. In the TSI spherical coordinate system, the relationship of the zenith angle, azimuth and scattering angle (marked with Θ) is shown in Figure 3. The small black square indicates a point on the spherical mirror, θ and ∅ are the zenith angle (IZA) and the azimuth angle of the point, respectively; and ∅ are the solar zenith angle (SZA) and SAA, respectively. The IZA is used to represent the zenith angle of the pixel on the image, and the SPA (Sun-Pixel Angle) is the angle between the sun and the pixel, representing the positional relationship between the pixel and the sun on the image. The SPA is approximately equal to Θ due to the long distance between the sun and the earth. From Figure 3, we have:

Cloud Detection and Retrieve
Feature extraction from cloud images is a key component in the application of the Total Sky Imager [12]. Although TSI-880 provides cloud detection according to the built-in algorithm, the single fixed red-blue-Ratio (RBR) threshold for all images and untreated shielding caused significant error in the results. To be specific, on the one hand, the cloud information near the sun shielded by the shadow band is vital for ultra-short-term irradiance forecasting [34]. On the other hand, the image RBR is greatly affected by atmospheric conditions including transparency, aerosol content and sand dust, even though for the same clear sky image, the RBR values around the sun and near the horizon are larger due to aerosol forward scatter and long optical path [40,41]. Given this, we performed cloud detection based on the RBR threshold to determine the cloud fraction.
The first step of cloud detection is recovering the shielding areas that result from the support arm and shadow band, as shown in Figure 2a. The angle of the shadow band deviates from the true north on the image, which is equal to the solar azimuth (SAA), because projection of the sky on the image is isometric [42], and the area of shadow band is shown in the green area of Figure 2b and interpolated according to the pixels on both sides. The projection of the support arm with a fixed position is shown in the white area of Figure 2b and is restored by horizontal bidirectional interpolation. The restored image is shown in Figure 2c. Additionally, the sky, effectively occupying 260 × 260 pixels, is cropped from the original image and distortion is corrected by the spherical distortion correction model [43].
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 17 data were sampled synchronously with the frequency of 1 min, and simultaneously stored locally and remotely to guarantee the timeliness, integrity and reliability of the data.

Cloud Detection and Retrieve
Feature extraction from cloud images is a key component in the application of the Total Sky Imager [12]. Although TSI-880 provides cloud detection according to the built-in algorithm, the single fixed red-blue-Ratio (RBR) threshold for all images and untreated shielding caused significant error in the results. To be specific, on the one hand, the cloud information near the sun shielded by the shadow band is vital for ultra-short-term irradiance forecasting [34]. On the other hand, the image RBR is greatly affected by atmospheric conditions including transparency, aerosol content and sand dust, even though for the same clear sky image, the RBR values around the sun and near the horizon are larger due to aerosol forward scatter and long optical path [40,41]. Given this, we performed cloud detection based on the RBR threshold to determine the cloud fraction.
The first step of cloud detection is recovering the shielding areas that result from the support arm and shadow band, as shown in Figure 2a. The angle of the shadow band deviates from the true north on the image, which is equal to the solar azimuth (SAA), because projection of the sky on the image is isometric [42], and the area of shadow band is shown in the green area of Figure 2b and interpolated according to the pixels on both sides. The projection of the support arm with a fixed position is shown in the white area of Figure 2b and is restored by horizontal bidirectional interpolation. The restored image is shown in Figure 2c. Additionally, the sky, effectively occupying 260 × 260 pixels, is cropped from the original image and distortion is corrected by the spherical distortion correction model [43]. In the TSI spherical coordinate system, the relationship of the zenith angle, azimuth and scattering angle (marked with Θ) is shown in Figure 3. The small black square indicates a point on the spherical mirror, θ and ∅ are the zenith angle (IZA) and the azimuth angle of the point, respectively; and ∅ are the solar zenith angle (SZA) and SAA, respectively. The IZA is used to represent the zenith angle of the pixel on the image, and the SPA (Sun-Pixel Angle) is the angle between the sun and the pixel, representing the positional relationship between the pixel and the sun on the image. The SPA is approximately equal to Θ due to the long distance between the sun and the earth. From Figure 3, we have: In the TSI spherical coordinate system, the relationship of the zenith angle, azimuth and scattering angle (marked with Θ) is shown in Figure 3. The small black square indicates a point on the spherical mirror, θ and φ are the zenith angle (IZA) and the azimuth angle of the point, respectively; θ and φ are the solar zenith angle (SZA) and SAA, respectively. The IZA is used to represent the zenith angle of the pixel on the image, and the SPA (Sun-Pixel Angle) is the angle between the sun and the pixel, representing the positional relationship between the pixel and the sun on the image. The SPA Remote Sens. 2020, 12, 3671 4 of 17 is approximately equal to Θ due to the long distance between the sun and the earth. From Figure 3, we have: where r is the distance between pixel and image center on the image, and R is the radius of the image.
where r is the distance between pixel and image center on the image, and R is the radius of the image. The images with the SZA of 58° under clear-sky conditions were selected for representing the relationship of the RBR, IZA and SPA, as shown in Figure 4a. It is obvious that the RBR values decrease with the increase in SPA, while they increase with the increase in IZA. The areas of near horizon (IZA > 75°) and around the sun (SPA < 35°) have larger values of RBR in relation to other areas. Hence, the Clear Sky Library (CSL) was applied to eliminate the dependence of RBR on the areas around the sun and horizon [45]. To be specific, the historical clear images were selected as background images to provide RBR reference for images to be detected, and SZA, IZA, SPA, R, G, and B brightness values of historical clear images were saved as the data fields of the CSL. Additionally, we calculated the standard deviation distribution of the RBR of different clear sky images (Figure 4b) and found that the RBR variation had a large range from 0.1 to 0.18 in the area around the sun (SPA < 35°). Given this, the closest image in the CSL was selected for the images to be detected, to reduce the impact of aerosols concentration and sun position projection on the RBR. The images with the SZA of 58 • under clear-sky conditions were selected for representing the relationship of the RBR, IZA and SPA, as shown in Figure 4a. It is obvious that the RBR values decrease with the increase in SPA, while they increase with the increase in IZA. The areas of near horizon (IZA > 75 • ) and around the sun (SPA < 35 • ) have larger values of RBR in relation to other areas. Hence, the Clear Sky Library (CSL) was applied to eliminate the dependence of RBR on the areas around the sun and horizon [45]. To be specific, the historical clear images were selected as background images to provide RBR reference for images to be detected, and SZA, IZA, SPA, R, G, and B brightness values of historical clear images were saved as the data fields of the CSL. Additionally, we calculated the standard deviation distribution of the RBR of different clear sky images (Figure 4b) and found that the RBR variation had a large range from 0.1 to 0.18 in the area around the sun (SPA < 35 • ). Given this, the closest image in the CSL was selected for the images to be detected, to reduce the impact of aerosols concentration and sun position projection on the RBR.
here r is the distance between pixel and image center on the image, and R is the radius of the imag The images with the SZA of 58° under clear-sky conditions were selected for representing th lationship of the RBR, IZA and SPA, as shown in Figure 4a. It is obvious that the RBR valu crease with the increase in SPA, while they increase with the increase in IZA. The areas of ne rizon (IZA > 75°) and around the sun (SPA < 35°) have larger values of RBR in relation to oth eas. Hence, the Clear Sky Library (CSL) was applied to eliminate the dependence of RBR on th eas around the sun and horizon [45]. To be specific, the historical clear images were selected a ckground images to provide RBR reference for images to be detected, and SZA, IZA, SPA, R, G d B brightness values of historical clear images were saved as the data fields of the CS ditionally, we calculated the standard deviation distribution of the RBR of different clear sk ages ( Figure 4b) and found that the RBR variation had a large range from 0.1 to 0.18 in the are ound the sun (SPA < 35°). Given this, the closest image in the CSL was selected for the images to b tected, to reduce the impact of aerosols concentration and sun position projection on the RBR. As for the thresholds of different sky types, a certain number of images taken on different date ring the observation were selected as training samples. The thresholds shown in Figure 5 we rived from the PDF (probability distribution function) of RBR differences of pixels between th lected images of CSL and the sample images, and the thresholds of clear sky (T_clear) and thic ud (T_thick) were 0.067 and 0.225, respectively. As for the thresholds of different sky types, a certain number of images taken on different dates during the observation were selected as training samples. The thresholds shown in Figure 5 were derived from the PDF (probability distribution function) of RBR differences of pixels between the selected images of CSL and the sample images, and the thresholds of clear sky (T_clear) and thick cloud (T_thick) were 0.067 and 0.225, respectively.
A pixel was identified as thick cloud if the difference of RBR between the pixel and the corresponding pixel in the background image (∆RBR) was larger than 0.225; if the ∆RBR was smaller than 0.067, the pixel was identified as clear sky, otherwise the pixel was identified as thin cloud. An instance of cloud detection is shown in Figure 6e. Moreover, the sun parameter (SP) proposed by Pfiste et al. [16] was introduced to correct the type of the pixels around the sun. SP is a dynamic single threshold defined as the RBR average of the pixels around the sun (Figure 4a, −35 • < SPA <35 • ). If the sun is obscured by thick cloud (Figure 6a), RBR of the pixels around the sun are usually smaller than that of the clear day, leading to misidentification as clear sky ( Figure 6e). Consequently, when the RBR of the pixel around the sun was larger than the SP, the pixel was re-recognized as thick cloud ( Figure 6f). Finally, the cloud fraction of thick cloud ( f thick ) and thin cloud ( f thin ) are the ratios of the numbers of thick cloud and thin cloud to total pixels, respectively.
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 17 instance of cloud detection is shown in Figure 6e. Moreover, the sun parameter (SP) proposed by Pfiste et al. [16] was introduced to correct the type of the pixels around the sun. SP is a dynamic single threshold defined as the RBR average of the pixels around the sun ( Figure. 4a, −35° < SPA <35°). If the sun is obscured by thick cloud (Figure 6a), RBR of the pixels around the sun are usually smaller than that of the clear day, leading to misidentification as clear sky ( Figure 6e). Consequently, when the RBR of the pixel around the sun was larger than the SP, the pixel was re-recognized as thick cloud ( Figure 6f). Finally, the cloud fraction of thick cloud ( ) and thin cloud ( ) are the ratios of the numbers of thick cloud and thin cloud to total pixels, respectively.

Cloud Velocity and Cloud Map Forecasting
Forecasting of the cloud map is a prerequisite for estimating solar irradiance with total sky images. The MATLAB-based algorithm "mpiv (particle image velocimetry)" developed by Mori and Chang was applied to the two consecutive images for estimating the translational speed of the cloud [46]. The algorithm has been successfully applied to cloud velocity estimation in solar irradiance forecasting with all-sky images, infrared images and satellite images [32,47,48]. To avoid errors caused by cloud detection, the red channel of the raw images was extracted to calculate cloud speed instance of cloud detection is shown in Figure 6e. Moreover, the sun parameter (SP) proposed by Pfiste et al. [16] was introduced to correct the type of the pixels around the sun. SP is a dynamic single threshold defined as the RBR average of the pixels around the sun ( Figure. 4a, −35° < SPA <35°). If the sun is obscured by thick cloud (Figure 6a), RBR of the pixels around the sun are usually smaller than that of the clear day, leading to misidentification as clear sky (Figure 6e). Consequently, when the RBR of the pixel around the sun was larger than the SP, the pixel was re-recognized as thick cloud ( Figure 6f). Finally, the cloud fraction of thick cloud ( ) and thin cloud ( ) are the ratios of the numbers of thick cloud and thin cloud to total pixels, respectively.

Cloud Velocity and Cloud Map Forecasting
Forecasting of the cloud map is a prerequisite for estimating solar irradiance with total sky images. The MATLAB-based algorithm "mpiv (particle image velocimetry)" developed by Mori and Chang was applied to the two consecutive images for estimating the translational speed of the cloud [46]. The algorithm has been successfully applied to cloud velocity estimation in solar irradiance forecasting with all-sky images, infrared images and satellite images [32,47,48]. To avoid errors caused by cloud detection, the red channel of the raw images was extracted to calculate cloud speed

Cloud Velocity and Cloud Map Forecasting
Forecasting of the cloud map is a prerequisite for estimating solar irradiance with total sky images. The MATLAB-based algorithm "mpiv (particle image velocimetry)" developed by Mori and Chang Remote Sens. 2020, 12, 3671 6 of 17 was applied to the two consecutive images for estimating the translational speed of the cloud [46]. The algorithm has been successfully applied to cloud velocity estimation in solar irradiance forecasting with all-sky images, infrared images and satellite images [32,47,48]. To avoid errors caused by cloud detection, the red channel of the raw images was extracted to calculate cloud speed because "mpiv" can only be used for two-dimensional images, and the red channel of the image has higher contrast to clear sky and cloud than the blue and green channels. The interval of the two consecutive images was 1 min, and the search sub-window was defined as an area of 32 × 32 pixels. Two iterative recursive checks were performed using the MQD (Minimum Quadric Differences) algorithm to calculate the velocity vector field. Median filtering was used to remove discrete error vectors, and finally linear interpolation was used to correct the error vectors.
The cloud velocities were assumed to have the same direction and magnitude due to the limited sight of the TSI. For the clear sky area, the calculated velocities are usually zero or close to zero. Therefore, the calculated velocity vector fields were divided into two categories with the K-means algorithm and the larger value of cluster centers was selected as the representative velocity vector of the cloud. The red and blue scatters represent the velocity fields of the cloud and clear sky, respectively, and the black cross indicates the cluster center in Figure 7. The velocity vector of the cloud was u = 10.3, v = 7.9 (pixels/min); u and v represent the velocity of the horizontal and vertical direction, respectively.
Sens. 2019, 11, x FOR PEER REVIEW 6 r contrast to clear sky and cloud than the blue and green channels. The interval of the cutive images was 1 min, and the search sub-window was defined as an area of 32 × 32 pi terative recursive checks were performed using the MQD (Minimum Quadric Differen thm to calculate the velocity vector field. Median filtering was used to remove discrete e s, and finally linear interpolation was used to correct the error vectors. he cloud velocities were assumed to have the same direction and magnitude due to the lim of the TSI. For the clear sky area, the calculated velocities are usually zero or close to z fore, the calculated velocity vector fields were divided into two categories with the K-m thm and the larger value of cluster centers was selected as the representative velocity vect oud. The red and blue scatters represent the velocity fields of the cloud and clear tively, and the black cross indicates the cluster center in Figure 7. The velocity vector o was u = 10.3, v = 7.9 (pixels/min); u and v represent the velocity of the horizontal and ver ion, respectively. o obtain the forecasted cloud map, the retrieved cloud map was advected to the future presentative velocity vector. However, the information at the boundary of images is likely t dictable due to the limited observation view of TSI. Here, the boundary of the forecasted c as filled with the boundary information of the image 1 min ago. From Figure 8, we found xels at the edge were more error-forecasted because these areas use the boundary informa previous moment, and the pixels of thin clouds are also more error-forecasted because the s are difficult to be accurately detected in the raw image, especially in the area around the ixels of clear sky are likely to be detected as thin clouds due to the forward scattering of aero aze. Therefore, the accuracy of cloud forecasting is limited by the TSI observation view detection errors. To obtain the forecasted cloud map, the retrieved cloud map was advected to the future with the representative velocity vector. However, the information at the boundary of images is likely to be unpredictable due to the limited observation view of TSI. Here, the boundary of the forecasted cloud map was filled with the boundary information of the image 1 min ago. From Figure 8, we found that the pixels at the edge were more error-forecasted because these areas use the boundary information of the previous moment, and the pixels of thin clouds are also more error-forecasted because the thin clouds are difficult to be accurately detected in the raw image, especially in the area around the sun. The pixels of clear sky are likely to be detected as thin clouds due to the forward scattering of aerosols and haze. Therefore, the accuracy of cloud forecasting is limited by the TSI observation view and cloud detection errors.
of the previous moment, and the pixels of thin clouds are also more error-forecasted because the thin clouds are difficult to be accurately detected in the raw image, especially in the area around the sun. The pixels of clear sky are likely to be detected as thin clouds due to the forward scattering of aerosols and haze. Therefore, the accuracy of cloud forecasting is limited by the TSI observation view and cloud detection errors.

The Clear Sky Model and Clear Sky Index
Extraterrestrial solar radiation received by the upper boundary of the atmospheric layer changes with the Earth-Sun distance. For the observation site located at latitude ϕ and longitude λ, the extraterrestrial irradiance (I ets ) at a certain moment can be calculated by the solar constant and the position of the sun at the moment: S c is the solar constant (1367 W/m 2 ) [49], d is the DOY (day of year) of the year, 1 January is 1 and 31 December is 365. Diurnal variation of I ets is a cosine function of the zenithal angle [50]. The GHI exhibits uniform attenuation of I ets under conditions of clear sky, while it shows random and irregular fluctuations under conditions of broken sky with active clouds. Consequently, many irradiance forecasting models were developed from clear-sky radiation models. The existing clear-sky models such as Solis, Kasten and Ineichen typically developed from atmospheric radiative transfer models can simulate the GHI well, but usually require additional meteorological factors such as aerosol optical depth, ozone content, water vapor content and turbidity [51]. In this study, the clear-sky model was established between the I ets and the GHI of clear days during the observation though the least square method.
where I clk is the GHI under clear sky, a n are fitting coefficients. The three-dimensional radiation transmission process in the clouds and the cloud radiative effect of clouds are very complicated; it is difficult to quantify the forcing effect of clouds on radiation. The variation of GHI is assumed to be caused by changes in cloud fraction and type (thick cloud and thin cloud) due to the small observation scope. The clear sky index, defined as the ratio of GHI to I ets , is used to reflect the attenuation degree of entire atmosphere to I ets . Whereas in the studies of radiation forecasting, the clear sky index (k t ) is usually defined as the ratio of GHI to clear sky GHI [31], meaning that only the effect of the cloud is considered while the influence of other radiative forcing Remote Sens. 2020, 12, 3671 8 of 17 factors, including atmospheric molecules and aerosols, is excluded. Additionally, the attenuation degree of radiation by clouds with different optical thicknesses varies considerably. Thick clouds with low-altitude can reduce DNI to almost zero in a few seconds [33], while thin clouds with high-altitude may reduce GHI by 10% [14,52]. Hence, we established GHI estimation models based on the principle of the attenuation degree of cloud to irradiance varying with the distribution proportions of the different types of sky (thick cloud, thin cloud, and clear).

Two Forecast Models of GHI
Based on the retrieved and forecasted cloud fraction and clear sky index, the two rolling forecasting models of GHI using ternary linear regression (TLR) and BP neural network were established. Rolling forecasting refers to training the previous model with n values observed from time t-n to t, and the model is used for the forecasting of future time.
The TLR model constructed the relationships of k t and cloud fraction as shown in Formula (5), describing the attenuation proportion of three sky types (thick cloud, thin cloud and clear sky) 5 min before forecast time, for correcting the clear sky index dynamically.
where f thick , f thin , f clear are the fractions of thick cloud, thin cloud and clear sky. The subscript t − 5_t indicates 5 min from t − 5 to t. The regression coefficients obtained by formula 5 combined with forecasted cloud fraction were applied to estimate the k t and GHI of the short-term future.
where, ∆t represents the forecasting time scale (taken as 0~5 min in this paper). The subscript t + ∆t represents the forecasting time.
The BP (Back Propagation) neural network is a multi-layer feed-forward neural network trained by an error back-propagation algorithm which can obtain complex non-linear processing capability through complex mapping of simple non-linear processing elements. Its basic principle is to use the gradient descent method to minimize the mean square error between the actual output and the expected output of the network. It consists of the input layer, the hidden layer, and the output layer. The hidden layer consists of several neural units and each neural unit has a randomly assigned threshold and assigns a weight to the input layer and output layer. The calculation process includes the forward propagation of the signal and the backward propagation of the error. The error output is calculated in the direction from input to output, and the weight and threshold are adjusted in the direction from output to input. In this paper, a BP neural network is used to establish a ultra-short-term GHI forecast model. The training process of the prediction model is simply expressed as Equation (10): ANN u is an untrained BP neural network, L is the number of hidden layers, N is the number of neurons in each hidden layer. I represents inputs including f thick(t−10_t) , f thin(t−10_t) , I ets , and I clk , and T represents the target. ANN B is the trained BP neural network. The forecasting model is given by: where I n represents new inputs that do not contain the samples involved in training. Forecasting accuracy and operating efficiency are important indicators for evaluating the performance of the ultra-short-term forecasting model due to the demand of online rolling learning. This paper selected 3, 20, and 10 for L, N and the number of samples, respectively, by training the N, L, the number of iterations and sample error and considering the accuracy of forecasting results and the efficiency of running time.

Evaluation of Forecast Models
Pearson correlation coefficient (R) given by formula 10 is used to evaluate the correlation between the forecasted and measured GHI: RMSE (Root Mean Square Error) and rRMSE (relative RMSE) are used to indicate the deviation degree of the predicted and measured irradiance, and smaller values for RMSE and rRMSE indicate higher forecasting accuracy and smaller deviation between the forecasted and observed GHI.
MAE (Mean Absolute Error) is the average of absolute error of the predicted and measured GHI, which can better reflect the prediction error. MAE and rMAE (relative Mean Absolute Error) can be defined by: MBE (Mean Bias Error) is the average value of the predicted and measured GHI errors; MBE and rMBE (relative Mean Bias Error) are given by:

Results of GHI Forecasting
Clear or overcast days will always result in perfect forecasting, making them unsuitable for reflecting the ground irradiance response to changes in cloud fraction and evaluating the performance of forecasting models with the sky imager; hence, these days are eliminated. The 4 representative days with partial clouds were selected from November 2017, including 5th November, 7th November, 8th November, and 11th November.
Although nowcasting of GHI is not affected by the cloud motion, it can be used as a verification method for the accuracy of the cloud detection and the significance of the ternary linear regression equation. Hence, the nowcasting of GHI obtained from formula 5 was used for verifying the rationality of the TLR model in this paper. The results of GHI nowcasting for 5 November 2017 and 7 November 2017 are shown in Figure 9a, and Figure 10a showed that the correlation between forecasted GHI and observed GHI is almost 1, indicating that all the nowcasting GHI can match the observed GHI. Additionally, from Figures 9 and 10b,c, we also found that the forecasted GHI could adequately match the randomness and uncertainty of the observed GHI. In other words, the ternary linear regression equation based on the three different sky types has a significant regression effect, and can consequently be applied for further forecasting.  The scatter diagrams shown in Figure 11 are the result of GHI forecasting ranging from 1 to 5 min on 5 November 2017 and the various errors are given in Table 1. The points dispersion increases with the increase in time scale, with a validation of r = 0.999, 0.983, 0.967, 0.948 for 1 min, 2 min, 3  The scatter diagrams shown in Figure 11 are the result of GHI forecasting ranging from 1 to 5 min on 5 November 2017 and the various errors are given in Table 1. The points dispersion increases with the increase in time scale, with a validation of r = 0.999, 0.983, 0.967, 0.948 for 1 min, 2 min, 3 min, 5 min, respectively. From Table 1, we also found that the forecasting errors increase with the The scatter diagrams shown in Figure 11 are the result of GHI forecasting ranging from 1 to 5 min on 5 November 2017 and the various errors are given in Table 1. The points dispersion increases with the increase in time scale, with a validation of r = 0.999, 0.983, 0.967, 0.948 for 1 min, 2 min, 3 min, 5 min, respectively. From Table 1, we also found that the forecasting errors increase with the increase in the forecast time. The various errors during the 4 days were given by Tables 1-4. As expected from the previous studies, higher lead times resulted in higher errors. The average RMSE of the 4 selected days for 1 min, 2 min, 3 min, and 5 min forecasting are 39 (MBE and rMBE) showed negative values, indicating that more sky clear pixels were incorrectly assigned as cloud.   Table 5. It is obvious that the errors of RMSE, MAE and MBE of the two kinds of forecasting all increase with the increase in lead time, and the correlation coefficient decreases with the increase in lead time.
By taking the neural network dynamic rolling learning method, the comparisons between the forecasted and measured GHI on 5 November and 11 November 2017 are given by Figures 12 and  13. The points with larger deviation correspond to the moment when the radiation changes sharply. Bias values (MBE and rMBE) also showing negative values may be related to excessive boundary clouds.   Table 5. It is obvious that the errors of RMSE, MAE and MBE of the two kinds of forecasting all increase with the increase in lead time, and the correlation coefficient decreases with the increase in lead time.
Although the cloud fraction and clear sky index satisfy the multiple linear regression relation in a short time (several minutes), when using this linear relation to statistically extrapolate the future irradiance, the forecast performance is not as good as the neural network forecast model and the error is larger. The reason for this is that, under broken sky conditions with rapid changes of cloud fraction and speed, the nonlinear BPN-based model can reflect these changes more effectively. This is because neural networks have the ability to learn and can handle complex non-linear problems. Table 1. Forecasting error comparison between the ternary linear regression (TRL) model and the Although the cloud fraction and clear sky index satisfy the multiple linear regression relation in a short time (several minutes), when using this linear relation to statistically extrapolate the future irradiance, the forecast performance is not as good as the neural network forecast model and the error is larger. The reason for this is that, under broken sky conditions with rapid changes of cloud fraction and speed, the nonlinear BPN-based model can reflect these changes more effectively. This is because neural networks have the ability to learn and can handle complex non-linear problems.

Discussion
The two GHI forecasting models developed in this paper are mainly based on the ground observational characteristics of clouds from the Dunhuang region of China, considering the attenuation effect of clouds with different thickness and proportion on GHI. To apply them to other areas, cloud detection and forecasting technologies are universal and necessary adjustment based on observations for the clear sky model are needed. Meanwhile, the radiation transmission in cloudy atmosphere is an extremely complicated process and the interaction mechanism between radiation and clouds remains poorly understood. In fact, clouds with different heights have diverse attenuation effects on radiation. However, cloud base height is not available for a single TSI site. Unified processing of clouds of different heights is likely to increase the error of forecasting.
Moreover, in this study, the temporal-spatial range is limited by the observation range of the ground sky imager. The boundary information of forecasted cloud images comes from the original images and cannot represent the sky situation in the future. When the clouds move fast, the forecasting error of the cloud image increases and the forecasting time scale reduces.
Overall, the average rRMSE of the 4 selected days ranged from 5% to 36% based on the TLR model and ranged from 12% to 32% based on the BPN model. In the previous study, the forecasting with four all-sky cameras, which were claimed to be state of the art, reduced the rRMSE of GHI from 30.9% to 23.5% [11]; however, the computational efficiency and cost were also much higher.
Dunhuang is located in northwestern China; the weather in winter is usually dry. There is almost no typical synoptic process during the observation period, broken days are few and the optical thickness of cloud is small. On the one hand, thin high-altitude cloud is hard to be accurately detected in the sky images, and the error of cloud detection and retrieval limited the accuracy of irradiance forecasting. On the other hand, thin clouds are prone to move and dissipate quickly, making forecasting more difficult and the time scale smaller.
In view of the shortcomings in the current work, the prospects for future work are proposed: The localized network of TSI observation will be deployed to retrieve the cloud base height and extract the cloud information and join the multiple images from different sites at the same time, achieving the monitoring of cloud motion in a larger area and improving spatial-temporal scale and accuracy of refined forecasting of solar radiation.

Conclusions
Two rolling forecasting models of GHI were established in this study based on the observation of a TSI and radiometer at a PV plant in Dunhuang, China. The main technical work included cloud detection and forecasting. To our knowledge, this is first time that such models have been developed and compared. Although the complicated effect of clouds on radiation is simplified in the models, especially in the TLR model, by including the objective observational characteristics of clouds into the surface radiation based on clear sky radiation, we simulated the GHI under the partly cloudy sky conditions. The two models allow us to effectively model the decreasing or increasing trend and degree of local GHI from the deployment of the ground observation of clouds; this is of significant relevance to the scheduling, dispatching and regulation of PV plants' power.
For cloud detection, the area around the sun and near-horizontal area are difficult to distinguish accurately; the error of thin cloud, especially in the areas around the sun and the shading areas, is greater than thick cloud. The threshold of clear sky and thick cloud is 0.067 and 0.225, respectively, in this study.
The accuracy of cloud map prediction is limited by the accuracy of cloud recognition. Because of the limited view of TSI, there will be clouds moving outside of the image, causing errors in the forecasted cloud map, through advection transport, to increase with the increase in lead time, and causing the average rRMSE of the 4 selected days to range from 5% to 36% based on the TLR model and from 12% to 32% based on the BPN model.