Towards the intrahour forecasting of direct normal irradiance using sky-imaging data.

Increasing power plant efficiency through improved operation is key in the development of Concentrating Solar Power (CSP) technologies. To this end, one of the most challenging topics remains accurately forecasting the solar resource at a short-term horizon. Indeed, in CSP plants, production is directly impacted by both the availability and variability of the solar resource and, more specifically, by Direct Normal Irradiance (DNI). The present paper deals with a new approach to the intrahour forecasting (the forecast horizon [Formula: see text] is up to [Formula: see text] ahead) of DNI, taking advantage of the fact that this quantity can be split into two terms, i.e. clear-sky DNI and the clear sky index. Clear-sky DNI is forecasted from DNI measurements, using an empirical model (Ineichen and Perez, 2002) combined with a persistence of atmospheric turbidity. Moreover, in the framework of the CSPIMP (Concentrating Solar Power plant efficiency IMProvement) research project, PROMES-CNRS has developed a sky imager able to provide High Dynamic Range (HDR) images. So, regarding the clear-sky index, it is forecasted from sky-imaging data, using an Adaptive Network-based Fuzzy Inference System (ANFIS). A hybrid algorithm that takes inspiration from the classification algorithm proposed by Ghonima et al. (2012) when clear-sky anisotropy is known and from the hybrid thresholding algorithm proposed by Li et al. (2011) in the opposite case has been developed to the detection of clouds. Performance is evaluated via a comparative study in which persistence models - either a persistence of DNI or a persistence of the clear-sky index - are included. Preliminary results highlight that the proposed approach has the potential to outperform these models (both persistence models achieve similar performance) in terms of forecasting accuracy: over the test data used, RMSE (the Root Mean Square Error) is reduced of about [Formula: see text], with [Formula: see text], and [Formula: see text], with [Formula: see text].

Introduction

Due to the ability of a sky imager to provide a detailed cartography of the sky, forecasting the solar resource using such a device is gaining interest at a fast pace. However, for now, only a few studies are dedicated to forecasting Direct Normal Irradiance (DNI) from sky-imaging data. DNI (I), which can be defined as the direct irradiance received on a plane normal to the Sun, can be split into two terms, i.e. clear-sky DNI () and the clear-sky index () (1):

Clear-sky DNI is about the solar power received at ground level per unit of area, at a specific location, when there is no cloud between the Sun and the observer, i.e. the solar field. This is key information for CSP plant operators since it is directly related to the upper limit of the available solar energy. For its part, the clear-sky index is derived from the attenuation of DNI caused by sun-blocking clouds. It spans from 0, when the path to the observer is obstructed by a thick cloud, to 1, when there is no cloud between the Sun and the observer.

These few studies have been conducted by a group of researchers belonging to the University of California in San Diego (UCSD). First, Marquez and Coimbra proposed in 2013 a deterministic DNI forecasting model based on the cloud cover estimated in a region of the sky with potentially sun-blocking clouds (Marquez and Coimbra, 2013). In this model, clear-sky DNI is supposed to be constant () during the time interval of interest (i.e. from 10 a.m. to 2 p.m.). The forecast horizon is up to ahead. Cloud pixels are identified by using the hybrid algorithm developed by Li et al. (2011) whereas cloud motion is estimated using the PIV (Particle Image Velocimetry) method (Marquez and Coimbra, 2013). Clearly, performance is degraded by not taking changes in into consideration. Following this first study, a new DNI forecasting model has been developed by Chu et al. (2013). The authors decided on a forecast horizon up to ahead. is estimated using a eighth-order polynomial of the cosine of the Solar/Zenith Angle (SZA). Note that a third-order polynomial has proven to do the trick. Both the cloud map and cloud motion are estimated in the same way Marquez and Coimbra have done (Marquez and Coimbra, 2013). In addition, a hybrid forecasting model that makes use of two artificial neural networks has been developed: the first neural network is used in case of high variability whereas the second one is used when variability is low. The authors highlight that two artificial neural networks, instead of only one network dealing with all types of variability, allow DNI to be forecasted with much better accuracy. The cloud cover in regions of the sky having great potential for affecting ground DNI values as well as last-minute values of DNI are serving as model inputs. Only data that agrees with the following condition are considered: , that is to say . As a result, DNI is forecasted when the Sun is near the zenith. In that case, perspective effects are lower. A short time afterwards, a novel method for cloud tracking that makes use of TSI images to forecast DNI up to ahead is presented (Quesada-Ruiz et al., 2014). Cloud pixels are identified by using the hybrid algorithm developed by Li et al. (2011) and is computed from an eighth-order polynomial of the cosine of the solar zenith angle. The authors introduce both a sector method to detect the motion of potentially sun-blocking clouds and an adjustable-ladder method, which is based on a size-adjustable set of grid elements that focus on sky regions of greatest potential for affecting ground DNI values. The model performs better than when cloud motion is estimated using the PIV method, showing improved performance in DNI forecasting under various sky conditions. However, persistence is still the best in case of clear sky. In the last study (Chu et al., 2015a), a generic off-the-shelf, high-resolution fisheye dome network camera is used to capture sky images. Regarding clear-sky DNI, the model proposed in 2002 by Ineichen and Perez is used (Ineichen and Perez, 2002). A Smart Adaptive Cloud Identification (SACI) system (Chu et al., 2015b) is used to conduct cloud detection and a grid-cloud-fraction method is applied to sky images in order to obtain numerical information regarding potentially sun-blocking clouds. So, a multi-ANN platform with cloud information serving as exogenous input allows DNI ramps to be forecasted over the next . Over-fitting is prevented by using cross validation. A Genetic Algorithm (GA) is used to find the optimal ANN topology and identify the platform's parameters. Results show that the proposed forecasting platform outperforms the reference persistence model.

Despite being promising, many of the models to the short-term forecasting of DNI that make use of sky-imaging data do not provide very accurate forecasts. This is especially the consequence of the highly-complex behavior of clouds and, as a result, estimating cloud motion locally is still a bottleneck. One can also highlight that Standard Dynamic Range (SDR) images do not provide quality information in both the circumsolar area and the darkest parts of the sky. Another conclusion can be drawn from the studies conducted in the last few years: using machine learning techniques seems to be an option to consider in the forecasting of DNI (in particular, in the clear-sky index forecasting) but questions are still pending, in particular regarding model inputs. Overall, recent achievements indicate that there is room for improvement. That is why PROMES-CNRS has decided to develop its own system, i.e. a sky imager able to provide High Dynamic Range (HDR) images of the sky (the PROMES Sky Imager, or PSI), as well as a set of new algorithms to the intrahour forecasting of DNI. Forecasts of DNI are given by the combination of clear-sky DNI () forecasts and clear-sky index () forecasts. Clear-sky DNI is forecasted from DNI measurements (made by a rotating shadowband irradiometer), using an empirical model combined with a persistence of atmospheric turbidity, taking advantage of the fact that changes in this quantity are relatively small throughout the day in comparison to changes in DNI, even when the sky is free of clouds. For its part, the clear-sky index is forecasted from sky-imaging data, using an Adaptive Network-based Fuzzy Inference System (ANFIS). Such a system synergizes the two techniques by combining the human-like reasoning style of fuzzy systems with the connectionist structure of artificial neural networks. The fractional cloud cover is used as main input. Preliminary results highlight that the proposed approach has the potential to outperform persistence models (either a persistence of DNI or a persistence of the clear-sky index) in terms of forecasting accuracy.

The paper is organized as follows: Section 2 is about the complete methodology of DNI forecasting. First, this section surveys the state of the art of clear-sky DNI models. Then the approach used to forecast clear-sky DNI in real time is briefly presented. The PSI is described in the sequel, as well as the algorithm used for HDR image generation. Regarding cloud detection, a hybrid algorithm has been developed. The section ends with a comparative study made on various approaches and models, including persistence models, to the forecasting of DNI. In Section 3, the results are interpreted in terms of forecasting accuracy. The paper ends with a conclusion and an outlook to future work (Section 4).

Methodology

The complete methodology of DNI forecasting can be found in this section of the paper.

Clear-sky DNI forecasting

The experimental data used to develop and validate the approach to the forecasting of clear-sky DNI are first (briefly) described. This approach combines the well-known Ineichen and Perez's model (Ineichen and Perez, 2002) with a persistence of atmospheric turbidity, taking advantage of the fact that changes in this quantity are relatively small throughout the day in comparison to changes in DNI, even when the sky is free of clouds (Chauvin et al., 2016). The section of the paper also focuses on a state of the art about clear-sky DNI models.

Experimental data

The database is derived from measurements realized in Perpignan, at the PROMES-CNRS laboratory, using a Rotating Shadowband Irradiometer (RSI) (Figure 1) (CSP Services, 2015). The station is located in southern France, approximately 20 km west of the Mediterranean Sea. Winter is mild and summer is hot and dry. In addition, there is a lot of wind, often resulting in a cloudless sky. Typical uncertainties of the RSI are about . Other data like the air mass or solar angles have been computed. Overall, we have used 80 days of data, collected in Perpignan from July 30 to October 17, 2015.

Figure 1

The Rotating Shadowband Irradiometer (RSI) installed at PROMES-CNRS, in Perpignan (southern France).

The main difficulty in forecasting clear-sky DNI comes from the discontinuities in the signal caused by clouds and nightfall. To overcome this problem, clear-sky data (in that case, there is no disturbance to DNI at all) have been selected using an off-line detection algorithm, i.e. a wavelet-based multi-resolution analysis (Nou et al., 2016). Such a selection will contribute to the better understanding of the transient behavior of atmospheric turbidity during the day. Taking advantage of the local tendencies of atmospheric turbidity is expected to improve forecasting accuracy, in comparison to what is obtained with a persistence of DNI.

State of the art of clear-sky DNI models

Regarding clear-sky DNI, models are divided into two categories: radiative transfer models and empirical models (Engerer and Mills, 2015; Gueymard, 2012a; Gueymard and Ruiz-Arias, 2015). In radiative transfer models, clear-sky DNI is usually derived from the Beer-Lambert law and can be expressed as follows (2): where is the extraterrestrial solar irradiance, m is the relative optical air mass and is the total optical thickness of a cloudless atmosphere.

The total optical thickness of a cloudless atmosphere is expressed as the sum of the broadband transmittances due to Rayleigh scattering (), uniformly mixed gases absorption (), ozone absorption (), nitrogen dioxide absorption (), water vapor absorption (), and aerosol extinction (). It is key in accurately assessing the state of the atmosphere and, as a consequence, the amount of solar energy reaching the ground (Davies and McKay, 1982; Gueymard, 2008; Gueymard and Myers, 2008).

Among all the radiative models one can find in the literature, REST2 (Gueymard, 2008) has proven to assess clear-sky DNI with unsurpassed accuracy. Highly complex phenomena that involve interactions between atmospheric particles (e.g. aerosols and water vapor) and sunbeams are taken into account in that model. In addition, REST2 includes two spectral bands with distinct transmission and scattering properties. Although radiative transfer models often produce better clear-sky DNI estimates than empirical models, they need input data that might be not available at any time. Indeed, aerosol optical depth data are required but happen to be difficult to obtain, as well as rarely available (Gueymard, 2012b). As a consequence, radiative transfer models are not suitable for real-time applications.

Regarding empirical models, several levels of complexity can be distinguished. Usually, the simplest clear-sky DNI models are based on the solar zenith angle and empirical correlations derived from experimental measurements (Daneshyar, 1978; Meinel and Meinel, 1976; Paltridge and Proctor, 1976). Clear-sky DNI can also be computed from an eighth-order polynomial of the cosine of the solar zenith angle (Chu et al., 2013; Quesada-Ruiz et al., 2014) or taking into consideration the altitude of the site, in addition to the solar zenith angle (Laue, 1970). Accuracy of these models is negatively impacted by the lack of information about the state of the atmosphere. So, several studies have been dedicated these past few years to developing more efficient models based on atmospheric turbidity, which can be determined from broadband beam radiation measurements (Chow et al., 2011; Ineichen and Perez, 2002; Yang et al., 2014). However, in these models, information about the state of the atmosphere is not real-time information and is generally derived from mean values or approximations provided by solar energy services (SoDa, 2015). Because atmospheric turbidity can be quite variable over the days, accuracy of the models can be unsatisfactory, especially at high solar zenith angles.

Ineichen and Perez's model

We have decided on the empirical model proposed by Ineichen and Perez (2002) (3): where b is a function of the altitude of the considered site, is the extraterrestrial solar irradiance, m is the relative optical air mass, and is the atmospheric turbidity coefficient revised by Ineichen and Perez.

As it is well known, is relatively stable throughout the day and is less dependent on air mass than , i.e. the Kasten-reviewed Linke turbidity coefficient (Eltbaakh et al., 2012; Grenier et al., 1994; Kasten, 1988; Kasten and Young, 1989). is defined as the power per square meter reaching the top of the Earth's atmosphere on a surface normal to the beam. Due to the elliptical orbit of the Earth, is not constant and varies through the year as follows (4): where is the solar constant defined as the average extraterrestrial solar irradiance at mean Sun-Earth distance ( = 149 597 871 km). is set to (Gueymard, 2012a).

Regarding the computation of the position of the Sun, several algorithms can be found in the scientific literature, with specific computation times, accuracies, and validity dates. Among the best known algorithms, one can cite the Sun Position Algorithm, called SPA (Reda and Andreas, 2003, Reda and Andreas, 2004), developed at the National Renewable Energy Laboratory (NREL). It is considered as the reference algorithm. However, its computation cost is pretty high. This can be unacceptable to real-time applications and, as a result, alternative (faster) algorithms have been developed these past few years (Blanc and Wald, 2012; Blanco-Muriel et al., 2001; Grena, 2008; Michalsky, 1988). Among these algorithms, SG (for Solar Geometry) has been developed in the framework of the European Solar Radiation Atlas (Rigollier et al., 2000; Wald, 2007). In the present work, the ratio comes from the conservation of flux and is given by the SG2 algorithm, a refined version of SG, which is accurate and fast enough (Blanc and Wald, 2012). It allows the Sun-Earth distance to be computed at each time step.

The relative optical air mass (m) is given by the ratio of the optical path length of the solar beam through the atmosphere to the optical path through a standard atmosphere at sea level with the Sun at the zenith. For its accuracy at all zenith angles and because computations are very fast, the formulation proposed by Kasten and Young (Kasten and Young, 1989) has been selected (5): where SZA is the Solar/Zenith Angle.

Regarding atmospheric turbidity, yearly, monthly or daily mean values of are commonly used as model inputs. Due to the calculation method used, such approaches do not always match with real-time applications. In addition, because a new DNI measurement has to be linked to the situation of the sky, a real-time detection algorithm for clear-sky situations has been developed. It makes use of the last-known clear-sky situation and requires the maximum opacification speed of the atmosphere to be evaluated (Nou et al., 2016, 2017). Clear-sky DNI is then forecasted using the Ineichen and Perez's model, combined with a persistence of atmospheric turbidity. This takes advantage of the fact that changes in this quantity are relatively small throughout the day in comparison to changes in DNI, even when the sky is free of clouds.

Approach to the forecasting of clear-sky DNI

The approach relies on a persistence of atmospheric turbidity, taking advantage of the fact that changes in this quantity are relatively small throughout the day in comparison to changes in DNI, even when the sky is free of clouds. Since clear sky defines the nominal operating conditions of Concentrated Solar Power (CSP) plants, forecasts of clear-sky DNI are key information for power plant operators tasked with the management of those plants. Of course, forecasting clear-sky DNI is also needed to forecast DNI for all sky conditions. Forecasts are given by (6): where , with the forecast horizon, b is a function of the altitude of the considered site, is the extraterrestrial solar irradiance, m is the relative optical air mass, and is the last-known value of atmospheric turbidity given by the real-time detection algorithm developed by Nou et al. (2016, 2017).

Overall, for the available data (a total of 80 days), the Mean Absolute Error (MAE) is lower than whereas the Root Mean Square Error (RMSE) is lower than , regardless of the forecast horizon (). Overall, performance is excellent. Both MAE and RMSE becomes lower than if is regularly updated. Three representative examples of clear-sky DNI forecasts (in Perpignan, southern France) are given by Figure 2, with . Taking a look at this figure, one can note the ability of the model to be accurate in forecasting clear-sky DNI, under different sky conditions. In the first day (October 1, 2015), DNI shows a high variability (the sky is partly cloudy all day). In the second day (October 10, 2015), DNI shows a low variability (for a large part of the day, the sky is clear). Finally, in the last considered day (October 15, 2015), the sky is clear at the beginning of the day and becomes totally cloudy from about 15:00 ( from that hour).

Figure 2

30-min forecasts of clear-sky DNI in Perpignan, southern France (red lines). Blue lines are for DNI measurements. Top: October 1, 2015. Middle: October 10, 2015. Bottom: October 15, 2015.

PROMES sky imager

The PROMES Sky Imager (PSI) consists in a 4-megapixel color camera equipped with a fisheye lens and protected by a waterproof enclosure (Figure 3). Sky images are collected every 20 s at a resolution of pixels, with 10 bits per channel. This is a 2nd-generation prototype that is designed to outperform the first camera installed at PROMES-CNRS since mid-2013. Among the different improvements, we find: a temperature controller, a stronger neutral density filter, a faster acquisition frame rate (22 fps), a higher bit depth and, more generally, a full access to camera settings (Table 1).

Figure 3

The PROMES Sky Imager (PSI) installed at PROMES-CNRS, in Perpignan (southern France). It consists in a 4-megapixel color camera equipped with a fisheye lens and protected by a waterproof enclosure.

Table 1

Characteristics of the PROMES sky imager.

Camera	MV1-D2048C-96-G2
Sensor	CMOS CMV4000
Resolution	2048 × 2048 pixels
Pixel size	5.5 μm
Color depth	10 bits (“raw” format)
Format	1″
Dynamic range	≃60 dB (no HDR)
Frame rate	22 fps
Exposure time	0.024 ms–349 ms
Operating temperature	0 °C to 50 °C
Objective	FE185CO86HA-1
Filter	Edmund Optics ND2.0
Enclosure	Turtle S IP67
Temperature controller	Yes [0 °C ; 50 °C]
Dimensions (L × W × H)	160 × 185 × 300 mm
Cost	≃5000€

A geometric angular calibration of the camera has been performed in order to get the relationship between a given pixel point and its projection onto the unit sphere. We used the OcamCalib toolbox (Scaramuzza et al., 2006) to calibrate the camera. The camera orientation has also been corrected comparing the real Sun position with its position as detected on the image. Once the camera is calibrated and its orientation corrected, it is possible to get the position of the Sun on the image even during cloudy days and get the angle between the Sun and each pixel on the image.

The system is not equipped with a solar occulting device to reduce the light intensity reaching the sensor. Indeed, although this device improves the sky visibility by reducing pixel saturation, it occults the circumsolar area, which provides vital information about the very short-term solar irradiance fluctuations. To reduce the considerable amount of incident radiance on the CMOS (Complementary Metal Oxide Semiconductor) sensor, a neutral density filter of transmittance has been installed. The exposure Q of a sensor is defined as follows (Miyamoto, 1964) (7): where is the exposure time and is the incident radiance received by the sensor.

The incident radiance received by the sensor can be approximated as follows (8): where L is the scene radiance, is the filter transmittance, D is the diameter of the aperture, n and are the indices of refraction of object and image space, respectively, S is the surface of the sensor and Ω is the solid angle.

Now, in order to limit disturbances due to diffraction, the PSI has a constant (maximal) aperture. The diffraction pattern resulting from a uniformly-illuminated circular aperture has a bright region in the center, known as the Airy disk, which together with the series of concentric bright rings around is called the Airy pattern. For its part, the starburst effect is essentially due to using a small aperture and how light is entering the lens and wrapping around the aperture blades as it travels to the sensor. The smaller the aperture, the more pronounced the effect. As a result, cannot be changed and reducing the exposure time is the only way to attenuate the sensor's exposure (7).

HDR imaging

The dynamic range between direct sunlight and dim skylight () is more than 120 dB whereas, for a single exposure, the camera's dynamic range is typically around 60 dB. Images are thus called Low Dynamic Range (LDR) images. is formulated as follows (9): where is the maximum radiance emitted by the Sun (typical value is ) and is the minimum radiance emitted by a thick cloud (typical value is ).

As a result, one picture is not enough to fully capture the range of illuminations produced by the Sun and by clouds. The large saturated areas on LDR images induce a loss of information in the circumsolar area, but also errors in cloud detection and cloud motion estimation when clouds are close to the Sun. Using HDR-imaging techniques, the dynamic range can be increased, allowing details to be kept in both the circumsolar area and the darkest parts of the sky.

Merging process

HDR imaging consists in merging several images, taken at different exposure times, such as every pixel is correctly exposed in at least one image (Figure 4). The merging process is as follows: first, for each image, over- and under-exposed pixels are removed. Indeed, pixels having a low value are too noisy, whereas pixels near saturation suffer from the non-linearity of the camera response function. Then, a linear regression analysis is performed on the remaining pixels in order to map their values according to a reference frame. Finally, a weighted average of valid pixel values is used to generate the high dynamic range image of the sky. Each of the just-mentioned steps is described in the sequel.

Figure 4

Study of the linearity of the sensor's response: comparison between pixels' value in a LDR (Low Dynamic Range) image taken at two different exposure times. Correctly-exposed pixels shared by two LDR images taken at exposure times and are highlighted in red in Figure 4(c).

Acquisition of a sequence of LDR images

Since the maximum frame rate of the PROMES sky imager is 29 fps, it is possible to use as much as 29 LDR images to generate an HDR image. In this work, the number of LDR images has been set to 24, for a total acquisition time around 0.8 s. As a key point, the total acquisition time has to be short enough for the scene radiance to stay constant, i.e. for the clouds' motion between each LDR image to be negligible.

First, it is necessary to define both the minimum exposure time and maximum exposure time , so that each pixel is correctly exposed on at least one of the LDR images, regardless of the weather condition. Since the Sun is saturated even with the camera's minimum exposure time, one can easily deduce that must be equal to this value: . The maximum exposure time must be defined so that each pixel is correctly exposed for the darkest overcast skies. Tests showed that this is the case for . The theoretical dynamic range of the HDR image is thus (10): where et are the dynamic ranges of HDR and LDR images, respectively.

In this configuration, is close to . In order to maintain a good overlap between two successive LDR images and since the camera model is linear, each intermediate exposure time should be chosen to be proportional to the previous one. However, generating an HDR image using 24 LDR images taken at 24 different exposure times is time-consuming (an image is taken every 20 s). For this reason, it has been decided to use batches of LDR images taken at the same exposure time. That way, computation time is limited while preserving a high signal-to-noise ratio (due to the averaging in each batch). In the end, batches of four images are acquired, leading to the following intermediate exposure times (11): where is constant and equal to the ratio between two successive exposure times and .

The exposure times , , , , and used to generate the HDR image are listed in Table 2.

Table 2

Exposure times used for HDR image generation.

t_1	t_2	t_3	t_4	t_5	t_6
25.8 μs	93 μs	334 μs	1024 μs	4333 μs	15600 μs

Suppression of under- and overexposed pixels

Camera model

When the sensor's response is linear (this hypothesis is verified in the sequel), the value of a pixel j can be expressed as follows (Granados et al., 2010) (12): where and are the collected charge (in electrons) and the value of pixel j for an exposure time , respectively, g is the overall camera gain factor, is the readout noise and is the round-off operator corresponding to quantization.

The charge collected by pixel j is linked to the incident radiance by (13): where is the incident radiance on pixel j, is the pixel's gain and is the dark current.

Noises introduced during an image acquisition (Aguerrebere et al., 2014) can be classified into temporal noises, that lead to differences in pixel value between images after all acquisition parameters are left untouched (photon shot noise, dark current, readout noise), and spatial noises, that lead to measurement discrepancies between pixels exposed to the same light intensity (photo-response non-uniformity and dark current non-uniformity). Spatial noises are generally calibrated by camera manufacturers. Consequently, they are neglected in the present work. That is to say, coefficients and in equation (13) are supposed to be pixel-independent: and . Results obtained in the sequel will validate this hypothesis. From equations (12) and (13), the incident radiance on a pixel j can be written as follows (14):

Note that the approximation in equation (13) is imposed by temporal noises and quantization. Supposing that the radiance stays constant during the total acquisition time of LDR images, it can be deduced that , the value of pixel j for an exposure time , can be expressed as a function of as follows (15):

As a result, if the sensor's response is linear, then there is a linear relationship between pixels' value in images taken at different exposure times. Reciprocally, if a linear relationship appears between pixels' value of two images, then the sensor's response can be considered as linear.

Linearity of the sensor's response

Pixels' value for two successive LDR images is plotted in Figure 4. It can be seen that there exists a zone, shared by two successive LDR images, for which linearity is verified. It can also be noticed that this linearity is lost either when pixels are saturated, or when the Signal-to-Noise Ratio (SNR) becomes too low. Since these ‘non-linearity zones’ must be suppressed before combining the LDR images, they must be characterized. Saturation occurs when a pixel is fully charged. As the measured signal is encoded in 10 bits, the saturation value is equal to , where DN stands for Digital Number. However, because of the several noise sources introduced during image acquisition, the non-linearity upper limit is reached before pixel's saturation value. One can note that, due to thermal and photon shot noises, this upper limit is all the more lower that the exposure time is long (see, e.g., Figure 4(d), where the linearity is lost around 900 DN). After performing tests on several sequences of LDR images taken at the longest exposure time ( and ), the upper limit above which the linearity of the sensor's response is not verified anymore has been set to the following value (16):

Pixels having a value higher than are thus suppressed from the process of HDR image generation. SNR ratio decreases when pixels are underexposed, because in that case readout noise is not negligible anymore. Most underexposed pixels can be found on images taken at the shortest exposure times (see, e.g., Figure 4(a)). Tests have been performed on several sequences of LDR images taken at the shortest exposure times ( and ), in order to determine the lower limit under which the linearity of the sensor's response is not verified anymore. This lower limit has been set to (17):

Pixels having a value lower than are also suppressed from the process of HDR image generation. Figure 5 shows correctly-exposed pixels from two sequences of LDR images. Exposure times and are not plotted, because, in these images, correctly-exposed pixels are confined in the circumsolar area (typically, less than 10° from the Sun's center). This is also why the scatter plot in Figure 4(a) is thinner than those in Figure 4(b) and Figure 4(d).

Figure 5

Under- and overexposed pixels (in white) and correctly exposed pixels in two sequences of LDR images (PROMES sky imager).

Correction of remaining pixels' value

Since LDR images are now stripped down from under- and overexposed pixels, the sensor's response linearity is verified. So, the next step in combining LDR images is to adjust pixels' value so that it corresponds to a ‘reference exposure time’ , defined as the maximum exposure time . As a result (18):

Using equation (15), a function , linking (the value of pixel j at an exposure time ) to an estimated value (the value of pixel j that would be obtained at exposure time instead of ), is defined. This function is given by (19): where and are coefficients.

Correcting LDR images then consists in calculating for each pixel j and each exposure time . To do that, coefficients and have to be known. Theoretically, can be calculated using equation (11) and then can be estimated as follows (20): where denotes the value obtained by averaging on all pixels j.

In practice, however, equation (20) is not usable because an image taken at exposure time sometimes does not share any correctly-exposed pixels with the reference image (Figure 5). Nevertheless, chosen values of , and make sure that two successive LDR images in the sequence share a sufficient number of correctly-exposed pixels. LDR images are corrected using the 2-step algorithm described in the sequel:

Step 1. Using correctly-exposed pixels shared by two successive LDR images, estimates of and are obtained by linear regression.

Step 2. Noticing that (21): where links to (19).

Estimates of , denoted , are obtained as follows, from the estimated value of , denoted (22): while coefficients are estimated recursively using (23):

Even though , i.e. the ratio between two successive exposure times, is known (11), it is estimated by linear regression in Step 1. The reason is simple: it allows to compensate for uncertainties on the actual exposure time used by the camera, and also to mitigate noise fluctuations due to the fluctuation of the sensor's temperature.

Combination of LDR images

Lastly, the value of a pixel j on the HDR image, for reference exposure time , is obtained by averaging estimated values of each LDR image (24):

A diagram depicting the whole HDR generation algorithm is given in Figure 6. Regarding the contribution of each LDR image to the HDR image, for a particular line of pixels that goes through the Sun, one can notice the overlapping of each LDR image and the high dynamic range of the resulting HDR image (Figure 7). Notably, the sky imaging system has proven to be able to measure information as close as 0.9° from Sun's center.

Figure 6

Algorithm used to generate an HDR image from a sequence of LDR images.

Figure 7

Contribution of each corrected LDR image to the HDR image. Pixels whose value is plotted belong to a line going through the Sun.

HDR images generated using the two sequences of LDR images presented in Figure 5 are given in Figure 8. It should be stressed that the dynamic range of HDR images is too high for the images to be printed on a screen or on paper. A tone mapping technique is thus used to convert HDR images into images that can be visualized (Ledda et al., 2005; Reinhard et al., 2010). This technique requires a compression of the dynamic and some details are lost in the process: in particular, the Sun is smaller in the HDR image than it appears in images of Figure 8.

Figure 8

HDR images obtained from the two sequences of LDR images plotted in Figure 5 (PROMES sky imager).

Computation time of the HDR generation algorithm

The mean computation time of each step of the algorithm is given in Table 3. In this setup, the camera can deliver an HDR image every 3 s (acquisition and pre-processing of LDR images is 1.67 s whereas generation of the HDR image is 1.14 s), which leaves 17 s for other signal and image processing algorithms (keep in mind that all the processing must be performed in less than 20 s).

Table 3

Mean computation time of HDR generation algorithm (using an Intel Xeon E3-1220).

Step		Resulting number of images	Computation time
Acquisition and pre-processing of LDR images	Acquisition of LDR images	24	0.82 s
	Downloading of LDR images	24	0.53 s
	Averaging LDR images in each batch	6	0.32 s
			Sub-total: 1.67 s
Generation of the HDR image	Suppression of under- and overexposed pixels	6	0.24 s
	Computation of δˆ_i+1,i and βˆ_i+1,i	6	0.38 s
	Computation of δˆ_r,i and βˆ_r,i	6	Negligible
	Combination of LDR images	1	Negligible
	Other steps (e.g. storage)	1	0.29 s
			Sub-total: 1.14 s
			Total: 2.81 s

Cloud detection

Various techniques can be used in the detection of clouds on sky images. The most popular ones are fixed or adaptive thresholding techniques that make use of features computed from the RGB (Red/Green/Blue) components of an image. The features one can use to emphasize cloud patterns on sky images are based on the physical and visual properties of the sky. These techniques are widely used because of the simplicity they offer and the limited amount of computation time needed to perform pixel identification. However, such techniques suffer from the anisotropy of the clear-sky background, leading to false detections. Indeed, both the circumsolar area and the Sun are usually detected as clouds, whereas thin clouds are often identified as clear-sky pixels. Therefore, removing the clear-sky anisotropy from sky images would significantly improve performance of cloud detection algorithms.

Clear-sky intensity distribution

Modelling the clear-sky intensity distribution using the PSI allows clear-sky images to be generated. These images can then be used to remove the clear-sky background anisotropy and so improve the real-time detection of clouds significantly (see Section 2.4.3 for a description of the algorithm developed in the framework of the CSPIMP research project). Of course, this detection is key in the short-term forecasting of DNI.

The formulation proposed by Chauvin et al. is especially designed for improving performance of existing models in the circumsolar area (Chauvin et al., 2015). It relies on experimental observations of the radiance profile near the Sun (25): where is the solid angle subtended by pixel p, PZA is the Pixel/Zenith Angle, and SPA is the Sun/Pixel Angle. K, , , , , , and are coefficients estimated in real time from a selection of clear-sky pixels.

When tested over more than 2200 clear-sky images, corresponding to a solar zenith angle spanning from 24° to 85°, the new formulation outperforms a standard approach based on the All-Weather model (Perez et al., 1993) by 15% on the whole sky and more than 20% in the circumsolar area. Overall, it showed a very good ability to reproduce the clear-sky anisotropy (Chauvin et al., 2015).

Fractional cloud cover

In case of a perfect image segmentation, the optical depth of the sky, for pixel p, hereafter noted as , is equal to zero () if that pixel is a “clear-sky” pixel, equal to one () if that pixel is a “thick-cloud” (opaque) pixel, and spans between zero and one () if that pixel is a “thin-cloud” pixel, respectively. Let us define as the set of pixels verifying the following condition (26):

As a first example, is about the set of pixels defined by angles these pixels make with the Sun less than 10°. As another example, is about the set of pixels defined by angles these pixels make with the zenith less than 90° (i.e. pixels above the horizon). So, the fractional cloud cover () can be formulated as follows (27): where is the solid angle subtended by pixel p, PZA is the Pixel/Zenith Angle, and is the maximum angle of aperture for which is estimated.

Correcting by the solid angle is mandatory because the PSI uses an equidistant objective and, as a result, the solid angle subtended by each pixel p in the image is not always the same. Finally, the fractional cloud cover can be defined as the ratio of the sum of the solid angles subtended by clouds to the solid angle 2π sr of the whole sky (28):

If the sky (the PSI's field of view) is free of clouds, then . On the other hand, is for a totally cloudy sky and is for a partly-cloudy sky.

Hybrid algorithm

The algorithm developed for cloud detection (Figure 9) takes inspiration from the classification algorithm proposed by Ghonima et al. (2012) when clear-sky anisotropy is known and the hybrid thresholding algorithm proposed by Li et al. (2011) in the opposite case. First, if a clear-sky RGB image () can be synthetized (case 1), a fixed threshold () is applied to the feature, formulated as follows (29):

Figure 9

The hybrid algorithm developed for cloud detection. As stated in Section 2.4.3, this algorithm takes inspiration from the classification algorithm proposed by Ghonima et al. (2012) when clear-sky anisotropy is known and the hybrid thresholding algorithm proposed by Li et al. (2011) in the opposite case.

In that case, using several thresholds for different parts of the image is no longer necessary because most of the clear-sky background anisotropy is removed (Section 2.4.1). is a fixed threshold because changes in luminance in clear-sky parts of that image are mostly removed thanks to the real-time correction of anisotropy. If that clear-sky RGB image () cannot be synthetized and if the sky is considered to be partly cloudy in the previous iteration (in that case, , with the fractional cloud cover, as defined in Section 2.4.2), an adaptive threshold, calculated using the Otsu algorithm, is applied to the NRBR (Normalized Red/Blue Ratio) feature (case 2). Such feature allows the visual contrast between sky and clouds to be improved. In addition, it is known to be robust to noise (30):

In all other cases, a fixed threshold () is applied to that feature (case 3). When the clear-sky RGB image () is not available, performance is the same as using the hybrid thresholding algorithm proposed by Li et al. (2011). Overall, performance of the proposed cloud-detection algorithm is very good. However, one can remark some light reflections on the RGB images, leading to false detections. Many solutions to this problem can be found in the literature.

Figure 10 shows two examples of RGB images, clear-sky RGB images (), and associated cloud maps (). From a cloud map, the fractional cloud cover () can be computed using equation (28). As highlighted by Figure 11, the sky is clear in Perpignan about 35% of the time whereas it is totally cloudy 28% of the time (from a total of 660 000 images). This is key information for operators in charge of managing a CSP plant. Indeed, this can help planning maintenance procedures as well as investments. Keep in mind that Perpignan is located in southern France. It experiences a Mediterranean climate similar to much of the Mediterranean coastline of France. Winter is mild and summer is hot and dry. In addition, there is a lot of wind, often resulting in a cloudless sky.

Figure 10

Two examples of RGB images, clear-sky RGB images (RGB), and associated cloud maps (C).

Figure 11

The fractional cloud cover (i.e. the ratio of the sum of the solid angles subtended by clouds to the solid angle 2π sr of the whole sky) in Perpignan (southern France) (from a total of 660 000 images).

DNI forecasting models

In this section of the paper, the DNI forecasting models (, with the forecast horizon) included in the comparative study are presented. The two persistence models assume that DNI and the clear-sky index persist over the forecast horizon, respectively. The third model included in the comparative study is a deterministic model based on the fractional cloud cover (). As previously stated, can be defined as the ratio of the sum of the solid angles subtended by clouds to the solid angle 2π sr of the whole sky. In the last model, an Adaptive Network-based Fuzzy Inference System (ANFIS) (Jang, 1993), with serving as input, is used to forecast over the next . Because ANFIS is a widely-used and very powerful neuro-fuzzy architecture, it has been chosen for forecasting the clear-sky index ().

Persistence of DNI

The model that relies on a persistence of direct normal irradiance is the simplest model included in the comparative study. It allows forecasting DNI for all types of horizons and does not need any knowledge about the atmosphere. Having access to real-time measurements of DNI (I) is the only requirement (31): where , with the forecast horizon, and I is the measured DNI.

Persistence models achieve very good performance in case of stationary weather conditions, i.e. for very cloudy days or, under clear-sky conditions, when the Sun is near the zenith. In opposition, persistence models are not suitable for high-variability signals. Changes in DNI are in part the result of the apparent movement of the Sun. As a result, assuming that DNI persists over the forecast horizon removes much of the effect of diurnal solar variation. Of course, this impacts forecasting accuracy significantly.

Persistence of the clear-sky index

The apparent movement of the Sun being known, one can try to take its effect on direct normal irradiance into consideration. Given that DNI can be split into two terms, i.e. clear-sky DNI and the clear-sky index (1), forecasting DNI is about forecasting both terms (32): where , with the forecast horizon, is the clear-sky index and is clear-sky DNI.

So, as another option than a persistence of DNI, one can assume that the clear-sky index persists over the forecast horizon and as a result (33): where , with the forecast horizon, I is the measured DNI and is clear-sky DNI.

Combining equations (32) and (33) leads to the following expression (34): where , with the forecast horizon, I is the measured DNI and is clear-sky DNI.

Regarding the forecasting of clear-sky DNI, the approach described in Section 2.1 is used. Keep in mind that such an approach combines the Ineichen and Perez's model with a persistence of atmospheric turbidity (35): where , with the forecast horizon, b is a function of the altitude of the site, is the extraterrestrial solar irradiance, m is the relative optical air mass, and is the last-known value of atmospheric turbidity (Nou et al., 2016, 2017).

So, the DNI forecasting model based on a persistence of the clear-sky index can be formulated as follows (36): where , with the forecast horizon, I is the measured DNI, and is the last-known value of atmospheric turbidity (Nou et al., 2016, 2017).

Real-time access to DNI measurements as well as (good) atmospheric turbidity () estimates are required. One can remark that using that forecasting model (36) is the same at using the model that relies on a persistence of DNI (31) but correcting the forecasts by a factor allowing the apparent movement of the Sun to be taken into consideration. Such a methodological approach has proven to be efficient in case of very cloudy or clear skies, regardless of SZA (the Solar/Zenith Angle).

Deterministic model based on

The present model aims at taking advantage of the fact that, from a statistical point of view, there is a correlation between the fractional cloud cover () at time t and the clear-sky index () estimated a few minutes later (at time ). As a result, can serve as model input (37): where , with the forecast horizon, and is the fractional cloud cover.

From equation (1), one can formulate as follows the deterministic model based on the fractional cloud cover (38): where , with the forecast horizon, is the fractional cloud cover and is clear-sky DNI.

Regarding the forecasting of clear-sky DNI, the approach described in Section 2.1 is used (39): where , with the forecast horizon, b is a function of the altitude of the site, is the extraterrestrial solar irradiance, m is the relative optical air mass, and is the last-known value of atmospheric turbidity (Nou et al., 2016, 2017).

ANFIS-based model

The last DNI forecasting model, on which the present paper focuses, makes use of an Adaptive Network-based Fuzzy Inference System (ANFIS) (Jang, 1993). First, the principles of neuro-fuzzy modelling are briefly presented. Then, the variables used as model inputs and the tuning of the system (through a supervised learning process) are described. The proposed model allows the clear-sky index to be forecasted over the next .

Neuro-fuzzy systems synergize the two techniques by combining the human-like reasoning style of fuzzy systems (through the use of fuzzy sets and a linguistic model consisting in a set of if-then fuzzy rules) with the connectionist structure of artificial neural networks (Lin and Lee, 1996). A fuzzy inference system is composed of five functional blocks: a collection of fuzzy if-then rules, a data set defining the membership functions of the fuzzy sets used to design the rules, a decision-making unit allowing the inference operations on the rules to be performed, a fuzzification interface to transform the crisp inputs into degrees of match with linguistic values, and a defuzzification interface to transform the fuzzy results of the inference into a crisp output.

ANFIS maps inputs through input membership functions and associated parameters and then through output membership functions and associated parameters to output, is based on the Takagi-Sugeno inference system and can be used to interpret an input(s)/output map. The parameters associated with the membership functions changes through the learning process. Adjustment is facilitated by a gradient vector. It provides a measure of how well the fuzzy inference system maps the input/output data for a given set of parameters. When the gradient vector is obtained, any of several optimization routines can be applied in order to adjust the above-mentioned parameters and reduce error measure.

The Adaptive Network-based Fuzzy Inference System (ANFIS) uses as inputs the clear-sky index () at time t, the fractional cloud cover () at time t, and . It can be expressed as follows (40):

Note that we have used a training set composed of 40 days, a validation set composed of 20 days, and 20 consecutive days to test the model. In addition, we have decided on Gaussian membership functions and split the universes of discourse of , , and with two fuzzy sets only (using more fuzzy sets, forecasting accuracy is more or less the same). The learning process has allowed 8 fuzzy rules to be designed, for a total of 56 parameters in the system. From (40), one can formulate the ANFIS-based model in the following way (41): where , with the forecast horizon.

Regarding the forecasting of clear-sky DNI, the approach described in Section 2.1 is used (42): where , with the forecast horizon, b is a function of the altitude of the site, is the extraterrestrial solar irradiance, m is the relative optical air mass, and is the last-known value of atmospheric turbidity (Nou et al., 2016, 2017).

Results

Performance assessment

Figure 12 shows some 30-min forecasts of DNI produced by the four models included in the comparative study. First, one can note that using a persistence model, either a persistence of DNI or a persistence of the clear-sky index, a 30-min delay appears in the forecasts, which leads to large errors. Both the deterministic model based on , which is calculated on the whole image, and the ANFIS-based model give regularized forecasts, the latter providing the best results due to its ability to take advantage of expert knowledge via fuzzy rules. However, there is clearly room for improvement, especially when DNI variability is high.

Figure 12

30-min forecasts of DNI in Perpignan (southern France). Regarding clear-sky DNI, the approach described in Section 2.1 (this approach combines an existing empirical model with a persistence of atmospheric turbidity) is used. Blue lines are for DNI measurements. Top: October 1, 2015. Middle: October 10, 2015. Bottom: October 15, 2015.

Taking a look at Figure 13, which deals with the Root Mean Square Error (RMSE) value as a function of (over the 20 test days), one can note that the ANFIS-based model outperforms the other models in terms of forecasting accuracy, regardless of the forecast horizon. Both persistence models have similar accuracy: RMSE is approximately equal to , with , and equal to , with . Taking one of these two models as reference, RMSE is reduced of about using the ANFIS-based model, with . It is reduced of about , with . Regarding the deterministic model based on , it outperforms the persistence models for a forecast horizon longer than approximately .

Figure 13

DNI forecasting. RMSE (in W m−2) as a function of the forecast horizon (Δt ⩽ 30 min). Regardless of the horizon, the ANFIS-based model outperforms the other models in terms of forecasting accuracy.

Discussion

Overall, one can highlight the pertinence of the proposed DNI forecasting model (the ANFIS-based model). Using, on one hand, the empirical model proposed by Ineichen and Perez in 2002 (combined with a persistence of atmospheric turbidity) to forecast clear-sky DNI and, on the other hand, sky-imaging data to forecast the clear-sky index has proven to do the trick. As highlighted in this paper, artificial intelligence tools can be used to forecast the clear-sky index () from the fractional cloud cover (). However, the main limitation of such a model relies on considering the whole sky instead of only regions of interest, i.e. regions of the sky with potentially sun-blocking clouds, in the forecasting of the clear-sky index. As a consequence, future work will focus on developing new algorithms for cloud motion estimation and improving the use of sky-imaging data in DNI forecasting.

Estimating cloud motion from a sequence of sky images remains a really challenging task due to the non-linear phenomena of cloud formation and deformation, and the non-rigid motion and structure of clouds (Brad and Letia, 2002). In the literature, several methods for estimating cloud motion have been studied: pel-recursive techniques (Skowronski, 1999), Block Matching Algorithms (BMAs) (Song and Ra, 2000), optical flow algorithms (Chen and Mied, 2013) and fuzzy inference systems (Chacon-Murguia and Ramirez-Alonso, 2015). To date, block matching and optical flow algorithms are the most popular of all motion estimation techniques. In BMAs, successive sky images are partitioned into macro blocks of pixels. Then, for each block in the current image, BMAs seek the best matching block in the previous image, according to some matching criteria (Liaw et al., 2009). The Exhaustive Search (ES) block matching algorithm provides the best performance by matching all possible blocks within the search area. However, ES is hardly usable in real-time applications, since the search process takes significant computation time. In order to reduce the computation time of ES algorithm, fast block matching algorithms for motion estimation have been proposed: three step search (3SS) (Kulkarni et al., 2015), four step search (4SS) (Po and Ma, 1996) and diamond search (DS) (Zhu and Ma, 2000). For their parts, optical flow algorithms are based on a brightness constancy assumption on matching pixels in consecutive images: ideally, optical flow and motion field should be the same. This assumption gives a first equation for estimating two unknowns (horizontal and vertical velocities). This so-called Aperture problem is then regularized with a prior term, generally consisting in a smoothness constraint (Baker et al., 2011).

Regarding cloud motion estimation, a hybrid approach that combines block matching (for dominant motion patterns) and optical flow (for local motion patterns) algorithms will be designed for situations of moderate complexity. For more complex situations, artificial intelligence tools (i.e. machine learning techniques and algorithms) will probably be considered. As another possible improvement, clouds (based on pixel intensity) and situations will be classified. Finally, time-series forecasting techniques will be used to forecast the fractional cloud cover in the region(s) of interest (i.e. the region(s) of the sky with potentially sun-blocking clouds).

Conclusion

Clear-sky DNI is forecasted from DNI measurements, using the empirical model developed by Ineichen and Perez (2002) combined with a persistence of atmospheric turbidity (Nou et al., 2016, 2017). Such a method takes advantage of the fact that changes in this quantity are relatively small throughout the day in comparison to changes in DNI, even when the sky is free of clouds. Clear-sky situations are detected in real time, allowing the last-known value of atmospheric turbidity to be used as model input.

The clear-sky index is forecasted from sky-imaging data, using an Adaptive Network-based Fuzzy Inference System (ANFIS). Indeed, in the framework of the CSPIMP (Concentrating Solar Power plant efficiency IMProvement) research project, PROMES-CNRS has developed a sky imager able to provide High Dynamic Range (HDR) images. HDR images have details both in the circumsolar area and the darkest areas of the sky. A particular care has been taken to maintain a linear sensor response throughout the process, in order to study the sky's radiance from HDR images. Moreover, a hybrid algorithm that takes inspiration from the classification algorithm proposed by Ghonima et al. (2012) when clear-sky anisotropy is known and from the hybrid thresholding algorithm proposed by Li et al. (2011) in the opposite case has been developed to the detection of clouds. Performance is evaluated via a comparative study in which persistence models (either a persistence of DNI or a persistence of the clear-sky index) are included.

Preliminary results highlight that the approach to the forecasting of DNI on which the present paper focuses has the potential to outperform persistence models in terms of accuracy. Both persistence models have similar accuracy. Taking one of these two models as reference, RMSE is reduced of about using the proposed approach (i.e. the ANFIS-based model), with . It is reduced of about , with . The main limitation of such an approach relies on considering the whole sky instead of regions of interest, i.e. regions of the sky with potentially sun-blocking clouds, in the forecasting of the clear-sky index. As a result, future work will focus on developing new algorithms for cloud motion estimation and improving the use of sky-imaging data.

Author contribution statement

Julien Nou, Rémi Chauvin, Julien Eynard, Stéphane Thil, Stéphane Grieu: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Competing interest statement

The authors declare no conflict of interest.

Additional information

There is no additional information available for this paper.