Short-Term Solar Irradiance Forecasts Using Sky Images and Radiative Transfer Model

In this paper, we propose a novel forecast method which addresses the difficulty in short-term solar irradiance forecasting that arises due to rapidly evolving environmental factors over short time periods. This involves the forecasting of Global Horizontal Irradiance (GHI) that combines prediction sky images with a Radiative Transfer Model (RTM). The prediction images (up to 10 min ahead) are produced by a non-local optical flow method, which is used to calculate the cloud motion for each pixel, with consecutive sky images at 1 min intervals. The Direct Normal Irradiance (DNI) and the diffuse radiation intensity field under clear sky and overcast conditions obtained from the RTM are then mapped to the sky images. Through combining the cloud locations on the prediction image with the corresponding instance of image-based DNI and diffuse radiation intensity fields, the GHI can be quantitatively forecasted for time horizons of 1–10 min ahead. The solar forecasts are evaluated in terms of root mean square error (RMSE) and mean absolute error (MAE) in relation to in-situ measurements and compared to the performance of the persistence model. The results of our experiment show that GHI forecasts using the proposed method perform better than the persistence model.


Introduction
The use of solar power as a renewable energy source has grown substantially over the last few decades.However, the amount of solar irradiance that reaches the Earth's surface can vary significantly over relatively short time periods, resulting in fluctuation and intermittence of power generation.Therefore, it is crucial that we develop new techniques to accurately produce short-term forecasts of solar irradiance at the surface.
Historically, solar irradiance forecasts have relied heavily upon information retrieved via satellites and from Numerical Weather Prediction (NWP) models.Currently, NWP models generate forecasts in the range of 6 h to several days ahead [1,2] while satellite-based methods are generally used for 1 to 6-h period forecasts [3,4].However, both NWP and satellite imagery lack the spatial and temporal resolutions necessary for providing detailed information regarding high frequency fluctuations of solar irradiance.To overcome this issue, ground-based sky imagers have been developed and deployed as a means of capturing local sky images with highly-refined spatial and temporal resolutions for intra-hour solar irradiance forecasts [5].
Cloudiness is the most significant determinant of the amount of solar radiation reaching the Earth's surface [6].The use of sky imagers allows us to accurately identify local cloud cover and subsequently, predict relative atmospheric evolution through the analysis of time series images.This information can be used to forecast changes in solar availability within time frames of 10-15 min for most environments.There are currently two main methods employed when utilizing sky-imagers for solar irradiance forecasting: Global Horizontal Irradiance (GHI) or Direct Normal Irradiance (DNI) which is directly obtained from the empirical relationship with cloud cover [7], and forecast models that integrate cloud images with machine learning techniques.Chow et al. [8] used digitally-processed sky images from Total Sky Imager (TSI) to detect cloud cover.The cloud images, which were used to calculate motion vectors by a cross-correlation method, were propagated forward in time resulting in an intra-hour forecast of GHI.Marquez and Coimbra [9] also made use of TSI for intra-hour DNI forecast, rather than GHI.This work computed the cloud field with Particle Image Velocimetry (PIV), which uses the minimum quadratic difference method to calculate velocity vectors for each segment from successive pairs of sky images.Then, they employed grid cloud fractions in an area of interest to forecast DNI.Bone et al. [10] used a clear-sky DNI model and a cloud fraction prediction algorithm to generate DNI predictions.The clear-sky model parameters were adaptively estimated from identified clear-sky DNI measurements over a moving window.The sector-ladder method, proposed by Quesada-Ruiz et al. [11], that computes 1-D radial cloud velocity vectors on a polar grid was modified to detect ramp events in the DNI predictions.Some researchers believe that it is necessary to use non-linear characteristics of machine learning methods to predict stochastic processes of solar radiation fluctuations.The machine learning methods for solar radiation forecast are various, such as Artificial Neural Network (ANN) [12,13], Support Vector Machine (SVM) [14,15], and k-Nearest Neighbor (k-NN) [16,17].In the work of Chu et al. [13], the expected cloud cover indices were calculated using the method described by Marquez and Coimbra [9].GHI forecasts that were 5, 10, and 15-min ahead of time were produced by ANNs that used the cloud indices and GHI lagged data as inputs.Cheng et al. [15] incorporated history measured irradiance data with features extracted from the all-sky images to perform predictions.The bi-model forecasting, which is trained according to cloud obstruction conditions near the solar disk, was constructed by support vector regressions.Pedro et al. [17] proposed a distribution-free methodology based on the k-Nearest-Neighbors (kNN) algorithm to predict intra-hour GHI and DNI irradiances using ground telemetry and sky images.Although forecast models based on machine learning methods achieve a lower rate of forecast errors, they can only predict GHI or DNI for specific time horizons, rather than producing a continuous solar irradiance forecast.
During broken cloud conditions, when the direct component of solar radiation is nearly unaffected, GHI at ground level often exceeds expected clear-sky values due to the scattering of the cloud in the sky.Additionally, when the sun is occluded, diffuse solar irradiance is the main energy component reaching the surface.Therefore, more accurate retrievals of the direct and diffuse components of solar irradiance will largely benefit GHI forecasts.The main novelty of this study is that Direct Normal Irradiance (DNI) and Diffuse Horizontal Irradiance (DHI) are quantitatively forecast within a 10-min time horizon by mapping the DNI value and downward diffuse radiation intensity fields obtained from Radiative Transfer Model (RTM) into the prediction sky images.The DNI model value and downward diffuse radiation intensity fields under clear sky and overcast conditions are calculated using the modified DISORT Radiative Transfer Model [18].The corresponding prediction sky images are produced by the optical flow method described in Section 3.2.
The data used for solar irradiance forecast is presented in Section 2. Detailed methods are presented in Section 3, including sky image processing to detect cloud, cloud-motion field calculation via the optical flow method, solar irradiance simulation based on the modified DISORT RTM, and solar irradiance forecast through the use of prediction images in combination with corresponding DNI model values and the downward diffuse radiation intensity field.The results and discussion are presented in Section 4, and the conclusions are presented in Section 5.

Data
The Total Sky Imager (TSI) and Multifilter Rotating Shadowband Radiometer (MFRSR) are located at the Southern Great Plains (SGP) atmospheric observatory in Oklahoma, United States, which was the first field measurement site established by the Atmospheric Radiation Measurement (ARM) Climate Research Facility.
The TSI is located at 36.6060 • N, 97.4850 • W at an altitude of 316 m and is composed of a heated rotating hemispherical mirror with a down-pointing Charge Coupled Device (CCD) camera located above it.The mirror contains a sun tracking shadowband that continuously shields the mirror from direct sunlight in order to protect the camera sensor from the sun's reflection [9].Experimental Images of TSI are taken at 1 min intervals [19].
The MFRSR is a seven-channel radiometer with six passbands 10 nm Full Width Half Maxinum (FWHM) centered near 415, 500, 610, 665, 862, and 940 nm, and it contains an unfiltered silicon pyranometer [20].It is positioned several hundred meters from the TSI with coordinates of 36.6050• N, 97.4850 • W at an altitude of 318 m.It uses an automated shadowbanding technique to measure the total-horizontal, diffuse-horizontal, and direct-normal spectral irradiances through a single optical path [21].In this study, the MFRSR channel at 415 nm, where gaseous absorption is minimal, was selected for solar irradiance forecasting.The sampling time of the DNI coincided with that of the TSI images.We used one-minute averages of DHI obtained from our three-sampling data, each of which was acquired at a 20s interval, as a means of eliminating any abnormal values caused by the cloud edge effect [19].

TSI Image Processing
Considering that pairs of consecutive TSI images are processed by an optical flow method to compute a cloud motion field, the shadowband areas of TSI images are filled by interpolation, as shown in Figure 1b.
A cloud detection method is then used to identify cloudy regions within the TSI images.Each pixel of the image is transferred into the normalized red to blue ratio, as defined by Yang et al. [22]: where b and r refer to the values of blue and red channels, respectively.The entire image of λ intensities is partitioned into three regions in order to identify clouds, as illustrated in Figure 1c.The first region (red) is represented by the circular area with a radius of 45 pixels and is centered within the circumsolar position.The second region (green) is represented by the circle area of radius 80 pixels, but not including pixels in the first region.The third region (blue) includes any pixels not covered by the first and second regions.Each local threshold value, which segments the clear and cloudy part of the sky, is computed in the individual partitioned region using Ostu's method [23].The selected threshold value within the first region will be determined by whether or not the sun is occluded.Therefore, circularity is adopted to determine the sun's occlusion [24].The circularity is calculated for the circumsolar saturated pixels whose brightness is greater than 245 in the unsigned 8-bit image, using the following Equation (2) [25]: where C is circularity, S refers to the area and L represents the perimeter.Higher values of C imply that the possibility of the sun being occluded is low.When the sun is not blocked by clouds (C ≥ 0.3), the threshold value is less than 0.05 to reduce any adverse effects on cloud detection caused by higher scattering of the direct beam.While the sun is occluded (C < 0.3), the threshold value is greater than 0.08 to better identify various types of cloud.The cloud detection image is shown in Figure 1d.
Lastly, all TSI images and cloud images are projected to a flat rectangular space for the purpose of removing the image distortion caused by the convex mirror.The steps are explained in detail by Marquez and Coimbra [9].The pixel locations of the original image are converted into polar coordinates: where r I is the radial distance from the zenith pole (center of image) to the pixels' coordinate (x I , y I ), and φ is the azimuth angle with respect to the horizontal axis and refers to the image, I. Due to the image distortion, the relationship between r I and the tangent of the zenith angle (θ) is nonlinear.
The relationship between r I and θ is approximated by fitting pairs of the sun's apparent radial distance (r I,sun ) and the corresponding solar zenith angle (θ s ) to obtain a cubic polynomial.The polynomial coefficients, a n , are obtained using the Least Square Method for data computed from images of three selected clear sky days: a n = {−0.0540, 0.0140, 0.8100, −0.0106}.
An example of the projected image is shown in Figure 1e.

Computing Velocity Fields by the Optical Flow Method
The core idea of optical flow is to obtain the motion field between two successive frames of an image sequence while enforcing the following gradient constraint equation under the assumption of brightness constancy: where I represents the pixel values in color or gray scale; I x , I y and I t are the spatial and temporal partial derivatives of I at time t; and u and v denote the velocity magnitude of each pixel's motion vector in the x and y directions, respectively.The equation is an ill-posed problem because one equation with two unknowns, u and v, cannot be solved without additional constraints.In order to solve this problem, Horn and Schunck [26] (abbreviated to HS) introduced the global strategy into the optical flow objective function.They formulated the optical flow as an optimization problem by minimizing the objective function under the assumption of brightness constancy and global smoothness of motions.Black and Anandan [27] (abbreviated to BA) introduced a non-convex robust penalty function into the objective of optical flow.They formulated the data term and the regularization term for smoothness with a series of non-quadratic robust penalty functions.In this research, the optical flow method proposed by Sun et al. [28] was used to compute the motion field of sky images.They formulated a new weighted non-local term, which is added to the classical objective function descended from HS and BA, to prevent over-smoothing across boundaries.The formula by Sun et al. [28] is composed of four terms: the data term, the regularization term for smoothness, the coupling term and the weighted non-local term.Its spatially discrete form is as follows: where u and v are the horizontal and vertical components of the optical flow field to be estimated from images I 1 and I 2 ; û and v denote an auxiliary flow field; i, j indexes a particular image pixel location; u i,j and v i,j are elements of u and v, respectively; λ is a regularization parameter; λ c is a scalar weight; ρ D and ρ s are the data and spatial penalty functions; N i,j is the neighboring set of pixel (i, j) in a possible large area; and W i ,j i,j represents the likelihood of i , j belonging to the same surface as i, j.The weighting is defined by the spatial distance, color-value distance, and the occlusion state as where I(i, j) is the color vector in the lab space that mathematically describes all perceivable colors in the three dimensions (L for lightness and a and b for the color components green-red and blue-yellow); n c is the number of color channels; and the occlusion variable, O(i, j), is calculated as where d(i, j) is the piecewise divergence function defined as in which the flow divergence, div(i, j), is where ∂ ∂x and ∂ ∂y are, respectively, the horizontal and vertical flow derivatives.In practice, the weighting defined by Equation ( 7) in a 15 × 15 neighborhood is used for the boundary regions, while equal weights in a 5 × 5 neighborhood are used to compute the median in the non-boundary regions.The parameter setting is the same as that used by Sun et al. [28], consisting of = 3, λ c = 0.1, In addition to this method, several best practices were also used including a standard, incremental, multi-resolution technique to estimate flow fields with large displacements.In addition, the Rudin-Osher-Fatemi (ROF) structure texture decomposition method was used to preprocess the input image sequence to reduce the influence of lighting changes.Other methods included incremental warping and warping with the spline-based bicubic interpolation to obtain the lowest energy solutions, median filtering of intermediate flow results after each incremental estimation step to remove outliers, and a graduated non-convexity (GNC) scheme to minimize non-convex energies.In practice, a two-stage GNC scheme was used, with the penalty functions for the first and second stages being the quadratic penalty (ρ(x) = x 2 ), and the generalized Charbonnier penalty function (ρ(x) = x 2 + 2 a , ε = 0.001, a = 0.45), respectively.The output of a previous stage served as the initialization to the next stage.An example of a velocity field calculation is shown in Figure 1e.

Solar Irradiance Simulation with RTM
In order to quantitatively predict GHI by calculating the direct and diffuse components of solar irradiance, it is necessary to estimate the cloud optical depth in the scene.We used the RTM to compute the DHI at different cloud optical depths and compared them with retrieved values from the nearby MFRSR, thereby allowing for the appropriate cloud optical depth to be chosen for GHI prediction.

Computing the Zenith and Azimuth Angle of Each Pixel on the Sky Image
The cartesian coordinates with respect to the central position of the image were X = x − x 0 , Y = y 0 − y, where x 0 , y 0 are center coordinates of the image, and x, y are the pixel coordinates on the image measured from the upper-left corner for columns and rows.The radial distance measured in pixels from the center position within the image is defined as r = We assumed that the center of the image was the zenith pole.Referring to the method for determining the relationship between the radial distance (r) and the zenith angle (θ) of each pixel in Section 3.1, the zenith angle of each pixel can be computed from a fit cubic polynomial with the variable r: The azimuth angle of each pixel (φ), which is the clockwise angle starting at zero on the left hand side, can be calculated by:

Locating the Sun's Position on the Image
The solar position can be calculated for a specific location on Earth given the date and time.Using the algorithm provided by Reda and Andreas [29], the zenith angle (θ s ) and azimuth angle (φ s ) of sun at time instance i can be obtained.Then the coordinates of sun position (x s , y s ) on the image are determined by: where r s is the sun's radial distance from the center position on the image, which can be calculated according to the relationship with the solar zenith angle, as described in Section 3.1.

Solar Irradiance Simulation with RTM
Using the modified DISORT Radiative Transfer Model, which combines the delta-fit method with the Nakajima-Tanaka correction procedure, we can accurately and rapidly compute radiances in both the forward and backward directions [18].In addition, the direct and diffuse solar irradiances at the channel of nongaseous absorption wavelength (415 nm) can be simulated under various cloud optical depths with different solar zenith angles [30].In the simulation, we used the Air Force Geophysics Laboratory (AFGL) midlatitude model, the aerosol profile from MODTRAN, 64 streams in the Legendre polynomial expansion, and surface albedos of 0.036 representing normal vegetated surfaces [21].
The RTM was initially utilized to establish a lookup table of the total optical thickness of Rayleigh scattering, gaseous absorption, and aerosol extinction at different solar zenith angles and visibilities.When clouds are broken and the sun is not occluded, the total optical thickness can be determined from direct beam measurements through the Bouguer-Lambert-Beer Law: where I dir is the direct beam transmittance of the solar air mass (A 0 ), which is a function of the solar zenith angle; and τ is the total optical thickness of Rayleigh scattering, gaseous absorption, and aerosol extinction.The total optical thickness at different solar zenith angles calculated from MFRSR measurements is then compared with the lookup table to obtain the visibility at each time instance.The hourly average visibility is then derived from the calculated model parameter values of the direct and diffuse solar irradiances.When the direct solar irradiance is calculated from the model values obtained by RTM, the pixel representing the position of the sun on the image and its four neighbourhoods are used to determine whether the sun is blocked.The direct solar irradiance is calculated by: where F dir represents the computed direct irradiance; DNI clr model denotes the direct irradiance of RTM in the clear sky; and CF sun is the cloud fraction calculated by the pixel representing the sun's position and its four neighbourhoods.Figure 2 shows the calculated direct solar irradiance using the RTM on 7 July 2008.In Figure 2, the red solid line indicates the direct solar irradiance simulated by the RTM in the clear sky (clear sky DNI), the green dash line represents the direct solar irradiance computed by the RTM and the cloud image (Computed DNI), and the gray area is the measured value of the direct solar irradiance (measured DNI).When the diffuse solar irradiance is computed, the simulated diffuse radiation intensity fields under clear sky and overcast conditions with different cloud optical depths are initially obtained from the RTM.In accordance with the zenith and azimuth angles of each pixel on the image calculated in Section 3.3.1, the diffuse radiation intensity fields can then be mapped to the image.The diffuse radiation intensity field under broken cloud conditions can be obtained by: where pixel = clear sky I cld , where pixel = cloud , ( where I, I clr and I cld represent the diffuse radiation intensity under broken clouds, clear sky and overcast conditions, respectively.The cloud image obtained from Section 3.1 is used to determine whether the pixel contains cloud or clear sky. Figure 3a-c show the image-based diffuse radiation intensity fields under clear sky, overcast and broken cloud conditions, respectively, at 19:37 (UTC) on 7 July 2008.Therefore, the total downward diffuse irradiance can be calculated by: where F di f (τ c ) denotes the total downward diffuse irradiance with different cloud optical depths and µ denotes the cosine of the zenith angle of the pixel.Figure 4 shows the total downward diffuse irradiance with different cloud optical depths on 7 July.In Figure 4   In order to select the appropriate cloud optical depths for solar irradiance forecasts, the point sets consisting of the computed values of DHI for different cloud optical depths and the measured value of MFRSR were linearly fitted, respectively.The slope of the fitting line that is closest to 1 means that the difference between the computed value and the measured value is minimal.

Solar Irradiance Forecast
By using the optical flow method with two successive images at 1 min intervals, as described in Section 3.2, the prediction images and the cloud prediction images for time horizons of 1-10 min were initially obtained.Figure 1f,g show the prediction image in two coordinate systems.Figure 1f is the prediction image projected onto rectangular grid coordinates of the sky, and Figure 1g is the prediction image transformed from a rectangular grid to TSI image coordinates.
After that, the DNI and DHI were forecasted by the method described in Section 3.3 with the RTM and corresponding cloud prediction image.However, there are two scenarios that must be accounted for when using the cloud prediction image.One is that the solar irradiance cannot be forecast when the sun position on the image is not within the cloud prediction image.The other occurs when there is an area that cannot be covered by the prediction image, as shown in Figure 1g.The diffuse solar irradiance of this area is estimated by: where denotes the forecast diffuse irradiance of the uncovered area.Therefore, the forecast GHI is calculated by: where F GH I is the forecast value of GHI, and F predict DN I and F predict DH I denote the forecast values of the direct and diffuse solar irradiances from the cloud prediction image, respectively.

Result and Discussion
The experimental dataset comprised daylight hours observed on 20-21 June, and 7 July 2008.These days all had broken cloud conditions with periods of clear sky.However, the cloud fraction throughout the days varied.The cloud cover was least on 20 June and gradually increased on 21 June.The greatest cloud cover occurred on 7 July.A limited number of partly cloudy days is more useful when evaluating short-term forecasts than a prolonged period of many clear sky days and overcast days.Error averages would likely be affected due to the fact that clear sky days and overcast days contain very little forecast error.

Selecting Appropriate Cloud Optical Depths
The diffuse solar irradiance computed with the RTM was divided into two parts, depending on whether or not the sun was occluded, and was then respectively compared to measured DHI values.Figure 5 shows the linear fitting results of scatter points consisting of the DHI computed with RTM and the DHI values measured on 7 July.The anomalous points of the measured value of DHI were caused by the cloud edge effect.When the sun was not blocked, the slope of the fitting line with cloud optical depth of 15 was closest to 1 (as seen in Figure 5a).When the sun was blocked, the slope of the fitting line with cloud optical depth of 25 was closest to 1 (as seen in Figure 5b).Therefore, these two cloud optical depths were used as cloud parameters for solar irradiance forecasting.

Evaluating the Error of Solar Irradiance Simulation with RTM
The metric used to evaluate the simulation error was the Mean Absolute Error (MAE), defined as: where F(t) denotes the measured irradiance at time instance t, F(t) denotes the corresponding computed irradiance, and n denotes the number of simulated samples.Table 1 shows the MAE values of the DNI, DHI and GHI computed by the RTM over these three days.From Table 1, it can be seen that the MAE values of DNI, DHI and GHI all increased with increasing the cloud cover.The main reason for the DNI inconsistency is that the direct solar energy is transmitted through the clouds.The ground measurements in this condition were greater than zero.However, the value calculated by the RTM was zero in the above condition.Another reason for the larger error in DNI on 7 July is the existence of optically thinner clouds during the period from 17.5 to 18 h.The attenuation of direct solar irradiance by thinner clouds during this period resulted in measured DNI values being less than those measured under clear sky conditions (seen in Figure 3).This type of cloud is difficult to detect because it is blue-tinted.Therefore, the measured values are less than computed values under these conditions.There are two reasons for the error in DHI.First, the measured values are abnormally higher than the computed values due to the cloud edge effect.Moreover, the DHI values computed by the RTM are not identical to the measured values because of the change in cloud optical depths caused by the variance in cloud shape and size over time.In addition, the errors between the GHI data simulated by the sky images and the measured data are due to the fact that the TSI and the MFRSR are hundreds of meters apart.When the sun is at the edge of the cloud, there may be a situation where the cloud blocks the direct solar radiation to the MFRSR but not to the sky imager, or vice versa.

Solar Irradiance Forecast Performance
The persistence model is the simplest forecasting model and is often used as a baseline for the performance evaluation of other forecast models.Persistence forecasts are defined as only considering current deviations from clear sky conditions for subsequent forecast time horizons.The persistence forecast of the solar irradiance at the next time step is computed as [7]: where ∆t is the forecast time horizon; F(t) and F clr (t) denote the measured irradiance and the modelled clear sky irradiance at time instance t, respectively; F p (t + ∆t) represents the forecast value of the persistence model at each forecasting horizon; and F clr (t + ∆t) represents the modelled clear sky irradiance at the corresponding forecasting horizon.
The metrics used to evaluate the prediction error for each forecasting horizon were the Root Mean Square Error (RMSE) and MAE, as defined by: where F(t + ∆t) denotes the measured irradiance at time instance (t + ∆t); F(t + ∆t) denotes the corresponding forecast solar irradiance computed by the persistence model or proposed forecast model; and n denotes the number of prediction samples.Table 2 shows the RMSE and MAE of GHI calculated by Equations ( 23) and (24) for the proposed forecast model and the persistence model in the above three days.It is known that the persistence model produces accurate forecasts when the variability of solar irradiance is low.From Table 2, it can be seen that the persistence model errors on these three days showed an increasing trend, which indicates that the fluctuation of the solar irradiance is more frequent with an increase in cloud fraction.During these three days, the number of errors in the proposed forecast model for each forecast horizon were less than the number of errors in the persistence model, which means that the proposed forecast model improves the estimation over random variability.However, as the forecast horizon becomes longer, the errors are larger due to rapid cloud deformation, formation, or evaporation.The GHI forecasting error can be evaluated further by making comparisons with the persistence model error.The forecast skill (s) is computed as: where RMSE f and RMSE p represent the RMSE for the proposed forecast model and persistence model, respectively.Skill values greater than zero imply that the proposed forecast method outperforms the persistence model.Figure 6 shows forecast skill values for the entire evaluation period.In Figure 6, it can be seen that there were higher forecast skill values at forecast horizons of 2-5 min.During these periods, cloud deformation, formation, and evaporation are relatively small so the cloud prediction image could be used for a more accurate quantitative forecast of solar irradiance.The lower predictability of 1 min is due to the fact that the shadowband area of TSI is too wide to cause the loss of information around the sun.For example, when the cloud is in the shadowband area, this cloud may not be obtained for the interpolated image, which could lead to an incorrect cloud prediction image at the next forecasting time.Time-series plots of the measured values, persistence model, and 3-min forecasts with the highest forecast skill are shown in Figure 7. From these figures, it can be found that the forecasts showed considerable improvement in the persistence of predicting the effect of cloud transience on GHI, especially during the period where there was a high frequency of fluctuation.Meanwhile, we also found that when the sun was obstructed by the cloud, the GHI measured values changed with the cloud fraction and cloud optical depth.Therefore, compared to the method of giving a fixed GHI estimated value when the sun is predicted to be blocked, our proposed method can obtain more accurate prediction results through qualitatively calculating the GHI values with the use of the cloud prediction image and RTM.

Conclusions
This paper aimed to use the Radiative Transfer Model and its corresponding TSI prediction image to quantitatively forecast DNI and DHI for time horizons shorter than 10 min.The main processes were the following: (1) The prediction images were produced by image processing, which includes cloud detection, spatial transformation of TSI images, and cloud motion field calculation.The non-local optical flow method was able to capture the independent velocity vector of each pixel across successive images.(2) The direct and diffuse solar irradiances simulated by the RTM were compared with in-situ measurements to determine model parameters-visibility and cloud optical depth.One of the key steps was to produce the image-based diffuse radiation intensity field under broken cloud conditions.(3) The direct and diffuse solar irradiances were forecasted from DNI model values and the diffuse radiation intensity field with certain cloud optical depths obtained from the RTM and corresponding cloud prediction images.
Three partly cloudy days with different cloud cover fractions were chosen to demonstrate the methodology.The RMSE and MAE errors of the presented forecast model were less than those of the persistence model, validating the use of prediction images combined with the RTM for 1-10 min, short-term forecasting horizons.
At the same time, there are some limitations to the proposed method and the existing instruments: (1) the proposed method may generate large errors in the presence of optically thin cloud layers due to their blue-tinting effect; (2) all clouds in images containing the same optical depth values as in the RTM may cause a certain number of DHI errors; (3) the wide shadowband of TSI and the distance between TSI and MFRST could cause the prediction errors.However, all in all, these situations open new possibilities for optimizing the current method and designing a new instrument integrating the sky imager and ground-based radiometer in future work.

Figure 1 .
Figure 1.Main image processing steps: (a) Original 8-bit image; (b) interpolation image; (c) mask used for partitioning image into three regions to obtain local cloud classification thresholds; (d) cloud detection image; (e) image projected onto rectangular grid and velocity field computed by the optical flow method; (f) prediction image projected onto rectangular grid; (g) prediction image transformed from a rectangular grid to TSI image coordinates.The red region indicates the area that the prediction image cannot cover.

Figure 2 .
Figure 2. The simulation of direct solar irradiance on 7 July 2008.
, the black line indicates the measured value of the diffuse solar irradiance (measured DHI), the red solid line indicates the diffuse solar irradiance simulated by the RTM in the clear sky (clear sky DHI), the other solid lines of different colors indicate the diffuse solar irradiance simulated by the RTM under overcast conditions with different cloud optical depths (overcast DHI), and the short dash lines of different colors represent the diffuse solar irradiance computed from the cloud image and model value with different cloud optical depths (computed DHI).

Figure 4 .
Figure 4.The simulation of diffuse solar irradiance on 7 July 2008.
CF represents the cloud fraction computed by the cloud prediction image; F cld model and F clr model denote the diffuse irradiance of RTM under overcast and clear sky conditions, respectively; F cld predict and F clr predict denote the diffuse irradiance in the range of the prediction image under overcast and clear sky conditions, respectively; and F uncover DH I

Figure 5 .
Figure 5.The linear fitting results of scatter points for selecting the appropriate cloud optical depth: (a) the sun is not occluded; and (b) the sun is occluded.

Figure 6 .
Figure 6.The forecast skill over the persistence of the proposed forecast method for the three days.

Table 2 .
RMSE and MAE in W/m2-nm computed for forecasting horizons of 1-10 min.