Benjamin Pillot
email: benjaminfp@unisinos.br

Marc Muselli
email: marc.muselli@univ-corse.fr

Philippe Poggi
email: philippe.poggi@univ-corse.fr

Pierrick Haurant
email: haurant@univ-corse.fr

João Batista Dias
email: joaobd@unisinos.br

Development and validation of a new efficient SRTM DEM-based horizon model combined with optimization and error prediction methods

Keywords: horizon model, SRTM DEM, solar radiation maps, error prediction

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

At the present time, the international energy situation is composed of two related main issues: global warming and depletion of fossil resources (Ben Ahmed et al., 2011;[START_REF] Hegerl | Understanding and attributing climate change[END_REF]. Since the advent of the sustainable development concept at the end of the 1980s, many political decisions have been implemented around the world in order to reduce the environmental impact of human development [START_REF] Liébard | Aménagement urbain et développement durable en Europe[END_REF]. Thereby, the use of renewable resources, which emit very few Green House Gases (GHG) and are almost unlimited energies (Ben Ahmed et al., 2011), is one of the first objectives of the Kyoto Protocol established in 1997 (United Nations, 1998). Among many other kinds of renewable energies, direct solar energy is currently used through different conversion technologies mainly to produce heat or electricity (Arvizu et al., 2011), such as photovoltaics [START_REF] Luque | Handbook of photovoltaic science and engineering[END_REF], flat solar thermal collectors [START_REF] Mosallat | Modeling, simulation and control of flat panel solar collectors with thermal storage for heating and cooling applications[END_REF], concentrated solar power [START_REF] Viebahn | The potential role of concentrated solar power (CSP) in Africa and Europe -a dynamic assessment of technology development, cost development and life cycle inventories until 2050[END_REF], solar cooking [START_REF] Wentzel | The development impact of solar cookers: a review of solar cooking impact research in south africa[END_REF] or solar air conditioning (Bermejo et al., 2010).

In order to predict reliability of these systems for supplying energy to rural or urban populations, it is necessary to accurately assess both the significance and distribution of the resource, i.e. to map the solar potential. The consistent solar mapping from ground measurements requires the distance between measuring sensors to be very low [START_REF] Muselli | Utilization of meteosat satellite-derived radiation data for integration of autonomous photovoltaic solar energy systems in remote areas[END_REF][START_REF] Perez | Comparing satellite remote sensing and ground network measurements for the production of site/time specific irradiance data[END_REF], which almost never happens [START_REF] Paulescu | Solar radiation measurements[END_REF]. Therefore, many models have been developed over the past years in order to retrieve solar radiation reaching the ground from satellite observations, like the Heliosat model (Beyer et al., 1996;[START_REF] Rigollier | The method heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF], or the ones developed by [START_REF] Brisson | Development of algorithms for surface solar irradiance retrieval at O&SI SAF low and mid latitudes[END_REF], [START_REF] Zelenka | Effective accuracy of satellite-derived hourly irradiances[END_REF], [START_REF] Perez | A new operational model for satellite-derived irradiances: description and validation[END_REF] or [START_REF] Janjai | Development of a method for generating operational solar radiation maps from satellite data for a tropical environment[END_REF].

Nevertheless, if the accuracy of these models is now well established [START_REF] Paulescu | Solar radiation measurements[END_REF], all of them only estimate ground solar irradiance from the extraterrestrial radiation and a physical characterization of the atmospheric transmittance, and leave out, because of their inherent spatial resolution, interaction of the radiation with the local topography.

Some models, part of Geographic Information Systems (GIS) like ArcGIS or GRASS, already take into account local terrain effects in the computation of solar radiation [START_REF] Ruiz-Arias | A comparative analysis of DEM-based models to estimate the solar radiation in mountainous terrain[END_REF]. However, these models are mainly based on simplified atmospheric parameterizations [START_REF] Ruiz-Arias | Spatial disaggregation of satellite-derived irradiance using a high-resolution digital elevation model[END_REF][START_REF] Ruiz-Arias | A comparative analysis of DEM-based models to estimate the solar radiation in mountainous terrain[END_REF], which still limits their consistency for developing regional solar map databases. As a result, in order to take benefit from the precision of the satellite-derived computation of the atmospheric transmittance, some post-processing procedures integrating topographic effects into satellite-based radiation maps have been proposed [START_REF] Bosch | Solar resources estimation combining digital terrain models and satellite images techniques[END_REF][START_REF] Haurant | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica[END_REF][START_REF] Pillot | The first disaggregated solar atlas of Djibouti: A decisionmaking tool for solar systems integration in the energy scheme[END_REF][START_REF] Ruiz-Arias | Spatial disaggregation of satellite-derived irradiance using a high-resolution digital elevation model[END_REF]. Finally, comparing estimates with ground data in mountainous areas, studies of [START_REF] Ruiz-Arias | Spatial disaggregation of satellite-derived irradiance using a high-resolution digital elevation model[END_REF], [START_REF] Bosch | Solar resources estimation combining digital terrain models and satellite images techniques[END_REF] and [START_REF] Haurant | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica[END_REF] have shown this method improves geographic information and accuracy of the final solar radiation maps.

In order to evaluate terrain effects, this process requires the use of a numerical topographic map of the study region, i.e.

a digital elevation model (DEM). First and main step of the procedure consists in retrieving the local horizon for every pixel of the DEM, in order to estimate diffuse and direct shading effects caused by the terrain onto the radiation [START_REF] Dozier | Rapid calculation of terrain parameters for radiation modeling from digital elevation data[END_REF][START_REF] Quaschning | Irradiance calculation on shaded surfaces[END_REF]. Thus, the higher the DEM resolution, the better the estimate accuracy, but also the longer the global computation, which is a very important parameter to take into consideration when the area of interest gets large, such as a small country [START_REF] Pillot | Planification de l'électrification rurale décentralisée en Afrique subsaharienne à l'aide de sources renouvelables d'énergie : le cas de l'énergie photovoltaïque en République de Djibouti[END_REF][START_REF] Pillot | The first disaggregated solar atlas of Djibouti: A decisionmaking tool for solar systems integration in the energy scheme[END_REF]. As a result, the authors think the corresponding horizon model has to present both accuracy and fast computing features, with in addition the possibility of evaluating model-based as well as DEM-based errors made on horizon estimates.

Consequently, we have developed a new efficient horizon numerical model using the Matlab programming language. In order to be fully implementable with exactly the same accuracy at any world location, it is exclusively based on theoretical assumptions without any empirical approximations. We have also added to this model a characteristic viewing distance computation allowing processing optimization, i.e. reduction of the model's running time by knowingly altering model precision.

We have then compared resulting estimates with in situ data, collected during an original topographic measurement campaign, and with estimates from some other models in order to assess how relevant this new model was. At last, from the ground data comparison, we propose a first approximation method for predicting the DEM-based error achieved on horizon estimates.

This article is divided into 3 main parts: the first focuses on the characteristics of the DEM used in this work; the second depicts the model theoretical basis as well as the related horizon viewing distance and the proposed processing optimization method; the final part presents both the comparison and the DEM-based error prediction studies.

DEM features

Horizon modeling directly depends on which kind of DEM is employed and on which specific Earth's mathematical representation it is defined from. We therefore present here characteristics of the main DEM used at present time, retrieved from the global Shuttle Radar Topography Mission (SRTM), as well as the coordinate reference system in which elevation data are represented.

The Shuttle Radar Topography Mission (SRTM): an interferometric measurement of the Earth's topography

A lot of developed countries own large national cartographic databases which were used before 2000 to produce DEM.

However, these DEM presented some limitations such as, on the one hand, country-specific resolution, scale and geographic projection, and, on the other hand, a low coverage of the global land because of some regions without any good-quality topographic maps [START_REF] Farr | The shuttle radar topography mission[END_REF]. Therefore, the National Aeronautics and Space Administration (NASA) in collaboration with the National Geospatial-Intelligence Agency (NGA) prepared and achieved the SRTM, in order to propose a new global and consistent high-quality representation of the Earth's topographic surface.

The SRTM lasted 10 days in February 2000 and consisted in radar scanning, from the space shuttle Endeavour, of 80 % of the world emerged lands located between 60 • N and 57 • S [START_REF] Farr | The shuttle radar topography mission[END_REF][START_REF] Rabus | The shuttle radar topography mission -a new class of digital elevation models acquired by spaceborne radar[END_REF]. The use of an interferometric synthetic aperture radar (InSAR), operating in the C band (wavelength = 5.6 cm), has permitted the assessment of land topography by measuring phase difference between radar pulses reflected by the surface and received at the same time by 2 distinct antennas [START_REF] Farr | The shuttle radar topography mission[END_REF][START_REF] Jarvis | Practical use of SRTM data in the tropics: comparisons with digital elevation models generated from cartographic data[END_REF]). Earth's DEM surface computed from SRTM data is a regular altitude mesh at 1 resolution (SRTM-1 DEM, about 30 m near the equator), presenting an absolute vertical error lower than 9 m (90 % confidence interval).

First SRTM datasets were released in 2003 by the NASA, with the SRTM-1 DEM mesh available only for the United States1 . Indeed, the rest of the world is covered by the DEM at 3 resolution (SRTM-3 DEM, about 90 m near the equator), computed by averaging 3 × 3 pixels of the original SRTM-1 DEM product. Furthermore, a significant proportion of elevation data is still missing in the original DEM, corresponding to a total area of 796217 km 2 [START_REF] Jarvis | Hole-filled SRTM for the globe version 4[END_REF] and mainly resulting from the low contrast within the final image (sharp slopes) or from smooth areas (smooth water or sand) which reflect too little energy to the radar [START_REF] Farr | The shuttle radar topography mission[END_REF]. Several methods have been developed to fill in these voids, mainly based on interpolation algorithms and on the use of elevation data from other sources. These schemes have been integrated into a global diffusion tool developed by the Consortium for Spatial Information (CGIAR-CSI), available via the organization website [START_REF] Jarvis | Hole-filled SRTM for the globe version 4[END_REF]. This tool gives access to SRTM-3 DEM 5 • × 5 • tiles we have used in this work.

SRTM DEM coordinate reference systems

A DEM is a numerical mesh of the Earth's topographic surface where geographic coordinates of each pixel depend on the mathematical shape used to represent it. The Earth's shape is not a perfect sphere because of the flattening at the poles and the heterogeneous gravity field at the surface. Therefore, 2 coordinate systems have been used by the SRTM to produce and reference elevation data: the WGS84 ellipsoid of revolution for horizontal position (geographic, or geodetic, longitude and latitude) and the EGM96 geoid for vertical elevation (orthometric height).

Geographic coordinates: the WGS84 ellipsoid

The World Geodetic System 1984 (WGS84) is a geocentric, right-handed and orthogonal Earth coordinate system, where origin O is the globe's center of mass, z axis is the axis passing through the poles and x axis is the intersection between the reference meridian (Greenwich) and the plane normal to the z axis and passing through O (NGA, 2000). The WGS84 ellipsoid of center O is a mathematical approximation of the Earth's shape retrieved by revolution, around its semi-minor axis b, of an ellipse of semi-major axis a and flattening f = a-b a :

a = 6378137.0 m f = 1/298.257223563 (1)
The ellipsoid is fully defined by these 2 parameters, from which we can also express the first eccentricity e [START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF][START_REF] Burkholder | Geometrical geodesy. In: The 3-D global spatial data model -Foundation of the spatial data structure[END_REF]:

e 2 = a 2 -b 2 a 2 = 2 f -f 2 (2)
Let now M be a point of the Cartesian coordinate system (x, y, z) and M * the orthogonal projection of M onto the ellipsoid such as [START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF]. In the geographic coordinate system, M is defined by its latitude ϕ, angle between -→ n and the equatorial plane, its longitude λ, angle between the reference meridian and the meridian passing through M * , and its ellipsoidal height h along the normal. In order to retrieve the coordinates of M in the Cartesian space from the geographic space, it is firstly possible to express the normal to the ellipsoid as in the case of a sphere:

----→ M * M = h • - → n , with - → n the normal to the ellipsoid in M * (
- → n =                  cos λ cos ϕ sin λ cos ϕ sin ϕ                  (3)
Then, we can calculate the Cartesian coordinates x, y, z of any point M onto the ellipsoid from its geographic coordinates ϕ, λ, h [START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF][START_REF] Burkholder | Geometrical geodesy. In: The 3-D global spatial data model -Foundation of the spatial data structure[END_REF]:

                 x = (N + h) cos λ cos ϕ y = (N + h) sin λ cos ϕ z = (N(1 -e 2 ) + h) sin ϕ (4)
Where N is the length of the ellipsoid normal, i.e. the distance between z axis and M * along the ellipsoid normal -→ n [START_REF] Burkholder | Geometrical geodesy. In: The 3-D global spatial data model -Foundation of the spatial data structure[END_REF]:

N = a 1 -e 2 sin 2 ϕ (5)
Figure 1 shows the WGS84 ellipsoid in the Cartesian coordinate system. A(ϕ A , λ A , h A ) and B(ϕ B , λ B , h B ), with corresponding distances N A and N B from z axis, are 2 points of the Cartesian space such as ---→

A * A = h A • -→ n A and ---→ B * B = h B • -→ n B
, and will be used later to explain horizon modeling.

Altitude baseline: the EGM96 geoid

A geoid is a mathematical model defining a specific equipotential surface of the Earth's gravity field around the mean sea level (NGA, 2000), and is used as a baseline to retrieve altitude of the topographic surface. The Earth Gravitational Model to determine surfaces where it remains constant [START_REF] Lemoine | The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96[END_REF]. Undulations of this geoid referring to the WGS84 ellipsoid, i.e. the distance between both the mathematical shapes in any point, can be determined using the model provided by the NGA [START_REF] Nga | NGA/NASA EGM96, N=M=360 Earth Gravitational Model[END_REF]. It is thus possible, for any point of the DEM, to derive the geoid-based altitude or orthometric height H from the ellipsoidal height h considering the geoid undulation ∆ WGS 84-EGM96 [START_REF] Lemoine | The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96[END_REF]NGA, 2000):

φ B B a y x z λ B λ A n A γ AB b N B N A O A h A n B h B Reference meridian Equator A * B * φ A
h = H + ∆ WGS 84-EGM96 (6)
The gap ∆ WGS 84-EGM96 is positive or negative according to the geoid position above or under the ellipsoid, and remains between approximately -107 m and 85 m (NGA, 2000).

Horizon modeling

Obstruction induced by the Earth's topography modifies the local horizon and thus the solar radiation reaching the ground.

It is defined, in horizontal coordinates, by the azimuth γ and the angular elevation α.

Horizon theoretical model

Here we present the horizon model we have developed considering both the points A and B already described in Figure 1 and in Figure 2.

Ellipsoidal height

Our theoretical model is based on both the reference systems previously described (WGS84 and EGM96). Every point of a SRTM DEM is thus defined by its geographic coordinates and its orthometric height [START_REF] Farr | The shuttle radar topography mission[END_REF]. Nevertheless, if the angular elevation α AB depends on the height H B as well as on the height H A , they are not both computed with respect to the same baseline because of the variations of the gravity field near the Earth's surface (geoid). Therefore, in order to accurately evaluate α AB from the DEM, it is necessary to retrieve altitude of the 2 points with respect to a same regular surface such as the WGS84 ellipsoid. As described before, it is possible to derive ellipsoidal heights h A and h B from relation (6):

h A = H A + ∆ WGS 84-EGM96 (A) h B = H B + ∆ WGS 84-EGM96 (B) (7)
∆ WGS 84-EGM96 (A) and ∆ WGS 84-EGM96 (B) are geoid undulations with respect to the ellipsoid in A and B respectively, and are depicted in Figure 3.

B(γ AB , α AB ) (P) h A Geoid Topography h B H B α AB n A n B Ellipsoid ∆ WGS84-EGM96 (A) H A ∆ WGS84-EGM96 (B) B * A A * π/2 -α AB

Angular elevation

Figure 3 is a sectional view of the ellipsoidal surface along the direction AB where belong the geoid, the terrain topography and the horizontal plane P tangent to the ellipsoid in A (orthogonal to -→ n A ). The fastest method to calculate the angular elevation α AB , i.e. the angle between the vector --→ AB and the plane P, is probably to use the scalar product. Indeed, α AB can be computed coupling the scalar product between the normal -→ n A and the vector --→ AB to the angle between them:

-→ n A • --→ AB = -→ n A • --→ AB • cos π 2 -α AB (8)
Obstruction only exists if the angular elevation is positive; so, using the sine, we finally have:

α AB = max           arcsin           -→ n A • --→ AB --→ AB           , 0           (9)
In the Cartesian coordinate system, -→ n A is derived from the relation (3). Then, equation ( 4) computes the coordinates (x A , y A , z A ) and (x B , y B , z B ) of A and B, which can be used to express the vector --→ AB and its norm:

--→ AB =                  x B -x A y B -y A z B -z A                  and --→ AB = (x B -x A ) 2 + (y B -y A ) 2 + (z B -z A ) 2 (10)

Azimuth

In horizontal coordinates, azimuth γ is the rotation angle between the reference north direction and the direction of interest.

Azimuth γ AB of the point B with respect to the observation point A is represented on the ellipsoid of Figure 1 as well as on the horizontal plane of Figure 2. It is possible to retrieve the azimuth γ using the meridian and parallel line elements ds m and ds p depicted in Figure 4, such as [START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF]: [START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF].

tan γ = ds p ds m (
Each line element depends on the considered geographic coordinate system (λ, ϕ), which is not a symmetric parameterization, i.e. the same angular variation (dλ = dϕ) doesn't lead to the same travel onto the surface (ds p ds m ). A symmetric parameterization (λ, L) of the ellipsoid is thus defined as follows:

ds p ds m = dλ dL ( 12
)
Where L is the isometric latitude, undefined at the poles (ϕ = ±π/2), and which is retrieved from the geographic latitude ϕ by applying the first fundamental quadratic form to the ellipsoid of revolution [START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF][START_REF] Eisenhart | A treatise on the differential geometry of curves and surfaces[END_REF]:

L = ln        tan π 4 + ϕ 2 • 1 -e sin ϕ 1 + e sin ϕ e 2        (13)
So, finally, the azimuth γ AB is determined by integrating equation ( 12) along the loxodrome AB, i.e. the constant azimuth line between A and B:

γ AB = arctan λ B -λ A L B -L A (14)

Computing full horizon

We have just defined the angular elevation α AB and the azimuth γ AB of a DEM point B with respect to an observation point A. In order to compute the full horizon, it is therefore necessary to apply the process to all the terrain which can be viewed from this point. Once it is perfectly characterized in the horizontal coordinate system (γ, α), it is then possible to identify, for a complete 360 • azimuth rotation, every maximal angular elevation, and so to retrieve the horizon.

From relations ( 9) and ( 14), a DEM originally referenced in the geographic coordinate system (λ, ϕ, h) can be now represented in the new horizontal coordinate system (γ, α). Furthermore, InSAR technology collects surface altitude: the SRTM DEM is therefore an elevation regular mesh at 1 or 3 resolution, composed of square pixels referenced by the geographic coordinates of their geometric center [START_REF] Farr | The shuttle radar topography mission[END_REF][START_REF] Jarvis | Practical use of SRTM data in the tropics: comparisons with digital elevation models generated from cartographic data[END_REF].

As shown in Figure 5, the first step of the process consists in demarcating a new narrower DEM tile, characterized by longitude width ∆λ and latitude height ∆ϕ, around a given observation point A(λ A , ϕ A , h A ). The size of this new area will determine the horizon computation time and will be defined in the next section. Once the local terrain mesh is perfectly known, the angular elevation α of each pixel is retrieved applying equation ( 9) to the geographic coordinates of its center

M(λ, ϕ, h).
As described in Figure 5, a loxodrome γ AF drawn from A to any boundary point F passes through a finite number of pixels [p 1 , p 2 , . . . , p i ] with angular elevation [α 1 , α 2 . . . , α i ]. Along this line, the point B presents the highest elevation α B and thus belongs to the horizon observed from A. To detect the pixels belonging to the loxodrome, our model calculates for each pixel of center M(γ, α) azimuth of northwest, southwest and southeast corners γ NW , γ S W and γ S E using relation ( 14). Then, the azimuth line passes through the pixel if γ AF is between γ NW and γ S W along a meridian, and between γ S W and γ S E along a parallel. Once elevation of all the pixels of the line is computed, the maximum is easily retrieved and corresponds to the horizon elevation in the γ AF azimuth direction.

λ φ γ AF = cte A(λ A ,φ A ,h A ) Δλ Δφ B(λ B ,φ B ,h B ) M N γ = ct A(λ A ,φ A ) M(γ,α) M(γ NW ,α) M(γ SE ,α) M(γ SW ,α) M(λ,φ,h) F Meridian Parallel DEM tile used to compute horizon M E M O M S
Figure 5: Methodology for computing the horizon around a given point A. Within a specific DEM tile of size (∆λ, ∆ϕ) centered on A, and for each azimuth between 0 • and 360 • , it consists in retrieving the pixel with the highest elevation (here α B ) among all the pixels the given loxodrome passes through (here

γ AF ).
Finally, computing the full horizon consists in repeating this step for a set of azimuth lines [γ AF 1 , γ AF 2 , . . . , γ AF i ]. Obviously, the selected number of lines directly influences horizon accuracy and computing time: in this study, we have considered a horizon azimuth step of 1 • .

Characteristic horizon viewing distance and model optimization

We present here an innovative way to reduce the size of the DEM tile required for accurately computing a single horizon, as well as a related method allowing us to save computation time in the case of large geographic areas, by knowingly altering the precision of multiple estimated horizons.

Characteristic viewing distance

Most horizon numerical models, whether specialized as the Carnaval French software or included in a GIS as the r.horizon model [START_REF] Hofierka | Open source solar radiation tools for environmental and renewable energy applications[END_REF] part of the r.sun add-on of the GRASS GIS software [START_REF] Šúri | A new GIS-based solar radiation model and its application to photovoltaic assessments[END_REF], let the user define the maximal distance used to compute the horizon around the observation point. However, the Earth's roundness inherently limits the minimum altitude of any DEM point: beyond a given distance from the observation point, this point is then necessarily located below the horizontal plane P. We can observe this statement in Figure 3 where, beyond a specific angle between -→ n A and -→ n B , the angular elevation α AB turns negative.

Consequently, we have developed a characteristic viewing distance concept, i.e. the maximum distance beyond which any DEM elevation data is no longer useful to compute a given horizon with the highest possible accuracy. In our model, the DEM tile used to model the horizon has width ∆λ and height ∆ϕ, and is centered on an observation point defined in the geographic coordinate system (λ, ϕ). The aim was therefore to define a characteristic viewing area for the horizon computation, i.e. the largest efficient DEM tile of width ∆λ max and height ∆ϕ max beyond which angular elevation of any DEM pixel is necessarily negative or 0.

To do this, it is firstly possible to express relations allowing calculation of the angle ϕϕ A along a meridian (constant longitude), and the angle λλ A along a parallel (constant latitude), between any point M(λ, ϕ, h) of normal length N and the observation point A(λ A , ϕ A , h A ) of normal length N A , where h > h A , such as α AM = 0. The angular elevation being equal to 0 when the normal -→ n A is perpendicular to the vector --→ AM, it leads to the 2 following equations:

-→ n A • --→ AM λ=λ A = 0 (15) -→ n A • --→ AM ϕ=ϕ A = 0 (16)
From the revolution symmetry of the ellipsoid results an explicit solution to the equation ( 16):

|λ -λ A | ϕ=ϕ A = δλ(h)| ϕ = arccos (N A + h A ) cos 2 ϕ A -(h -h A ) sin 2 ϕ A (N A + h) cos 2 ϕ A (17)
Meanwhile, asymmetry of the meridian ellipse leads to an implicit formula for the angle ϕϕ A :

(N + h) cos(ϕ -ϕ A ) + (N A sin ϕ A -N sin ϕ)e 2 sin ϕ A -(N A + h A ) = 0 (18)
While this equation can be numerically solved (iteration, dichotomy), it is however more interesting to use a local approximation in order to obtain an explicit solution. This is done by regarding the meridian ellipse as a circle where the radius is the radius of curvature in the north-south direction [START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF][START_REF] Burkholder | Geometrical geodesy. In: The 3-D global spatial data model -Foundation of the spatial data structure[END_REF]. The middle term of the previous equation disappears (eccentricity equal to 0) and the lengths of the ellipsoid normal N A and N are replaced by the radius ρ of the local spherical approximation calculated in A:

|ϕ -ϕ A | λ=λ A = δϕ(h)| λ = arccos ρ + h A ρ + h ( 19
)
Where ρ is given by the following relation [START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF][START_REF] Burkholder | Geometrical geodesy. In: The 3-D global spatial data model -Foundation of the spatial data structure[END_REF]:

ρ = a 1 -e 2 1 -e 2 sin 2 (ϕ A ) 3 2 (20)
In that specific case, M can be any point located whether on the meridian or on the parallel passing through A. But we can easily generalize these 2 relations to any boundary point of coordinates (λ A ± δλ(h)| ϕ , ϕ A ± δϕ(h)| λ , h) since we will still have α AM ≤ 0. It is then possible to use this methodology to retrieve the size (∆λ max , ∆ϕ max ) of the largest efficient DEM tile used to build the horizon. Instances of boundary points in the 4 cardinal directions, M N (λ M N , ϕ M N , h M N ) to the north, and are retrieved by iteration. To do this, we firstly initiate heights h M N , h M E , h M S and h M W with the maximum theoretical value of the global DEM (Everest's altitude for instance). Then, we iteratively reduce the size of the optimal DEM tile, by looking for the maximum height in each cardinal direction within the remaining area, as long as one of the following conditions is still met:

                           0 ϕ M N -ϕ A δϕ(h M N )| λ & λ M W λ M N λ M E 0 ϕ A -ϕ M S δϕ(h M S )| λ & λ M W λ M S λ M E 0 λ A -λ M W δλ(h M W )| ϕ & ϕ M S ϕ M W ϕ M N 0 λ M E -λ A δλ(h M E )| ϕ & ϕ M S ϕ M E ϕ M N (21)
Finally, the width and height of the optimal DEM tile required to perfectly characterize the horizon are given by:

∆ϕ max = ϕ M N -ϕ M S ( 22
)
∆λ max = λ M E -λ M W (23)
Obviously, in mapping applications, it is not really efficient to define a characteristic area for every single computed horizon of a given region. But it is possible to define the optimal DEM of a specific region, i.e. the zone beyond which it is no longer necessary to look for elevation data to compute any horizon of the region with the highest possible accuracy. In that case, we no longer consider a single point but a geographic rectangle enclosing the boundaries of the study region; the optimal DEM tile is then retrieved by using the previous method and by regarding h A as the minimum height inside the region. This tile encompasses the region enough for allowing the computation of all horizons with the highest possible accuracy.

Reducing model's running time

Depending on the size of the considered region, the desired accuracy and the purpose of the study, it may be interesting to reduce the model's computation time in order to improve efficiency of the whole process. Moreover, now both the optimal DEM and the maximum characteristic area of a given region are defined, it is possible to knowingly decrease accuracy of a

given horizon by comparing it with the theoretical ideal one. We therefore propose an efficient method allowing, in the case of large geographic areas, computation time to be reduced and resulting horizon precision to be estimated.

In order to avoid an optimization as time-consuming as the area processing itself, it is relevant to use statistical inferences, which consist in assessing the most significant properties of the global error distribution through a given territory by analyzing smaller data sets sampled from it.

Regarding statistical inferences, estimating the parameter θ of a population of observations F consists in determining the estimator Θ, random variable equal to θ for each sample we randomly pull out from F [START_REF] Walpole | Probability & statistics for engineers & scientists[END_REF]. However, in order to accurately assess the standard error and the confidence intervals of the sampling distribution of Θ, the distribution of the population F also requires to be perfectly known [START_REF] Walpole | Probability & statistics for engineers & scientists[END_REF]. Here, it was therefore necessary to look for a non-parametric method, i.e. fully independent of the distribution, remaining accurate in spite of the low sample size.

A well-suited method to that kind of problem is the bootstrap method [START_REF] Efron | An introduction to the bootstrap[END_REF], which relies on the plug-in statistical principle. This principle consists in determining the parameter θ of any distribution F, such as θ = t(F), from the estimate θ of an empirical distribution F, which is a sample composed of n observations randomly drawn from F, such as θ = t F . In addition, it is also necessary to apprehend the corresponding error, and here comes the main idea of Let ε be the error between an horizon estimate retrieved for a given viewing area (∆ϕ, ∆λ) and the ideal theoretical horizon corresponding to the optimal area (∆ϕ max , ∆λ max ). Let then E k be the population of all the errors achieved on the calculation of each k-horizon, i.e. computed using a specific k-area of size (∆ϕ k , ∆λ k ), within a given region. In order to determine the main components of the resulting error, i.e. both the mean µ E k and the variance σ 2 E k of the population E k , a plug-in estimate of n independent errors Êk = {ε 1 , ε 2 , . . . , ε n } k is gathered by computing the errors of n k-horizons corresponding to n random DEM pixels.

First of all, in order to define the horizon viewing area of a given region of interest, it is relevant to consider a low number of observations (n ≤ 50) and to repeat the previous step for different k-areas. Then, by using a box plot representing the error against the computation time, it is possible to select the best compromise between them. As an example, Figure 6 shows the box plot we have used to determine the most relevant area size for computing all horizons of the Republic of Djibouti [START_REF] Pillot | Planification de l'électrification rurale décentralisée en Afrique subsaharienne à l'aide de sources renouvelables d'énergie : le cas de l'énergie photovoltaïque en République de Djibouti[END_REF]. for different k-areas [START_REF] Pillot | Planification de l'électrification rurale décentralisée en Afrique subsaharienne à l'aide de sources renouvelables d'énergie : le cas de l'énergie photovoltaïque en République de Djibouti[END_REF]. In that example, the error determined with respect to the ideal horizon is the root mean square error.

Once the desired k-area is known, it is then possible to evaluate the likely error performed on any corresponding horizon computed within the given region. In other words, it means estimating the mean µ E and the standard deviation σ E of the distribution E of the modeling error ε, which follows a probability law ε ∼ µ E , σ 2 E . To perform this, a new empirical Finally, a simple method for estimating confidence intervals of the population parameters (µ E and σ E ) is the percentile method. However, asymmetry and bias of the standard deviation distribution [START_REF] Walpole | Probability & statistics for engineers & scientists[END_REF]) make an approach such as the bias-corrected and accelerated (BC a ) procedure, developed by [START_REF] Efron | Better bootstrap confidence intervals[END_REF], better suited. The classical percentile approach gives, for a significance level ν, the (1 -2ν) % confidence interval of the given estimation θ as [ θ * (ν) , θ * (1-ν) ], where θ * (ν) is the 100ν th percentile of the distribution of the B bootstrap estimates θ * (1) , θ * (2) , . . . , θ * (B) . In order to correct both the bias and the asymmetry of the given estimator, the BC a method enhances this first approximation by using other interval boundaries, as [ θ * (ν 1 ) , θ * (ν 2 ) ] where ν 1 and ν 2 are adjusted with respect to ν and 1 -ν using the standard normal cumulative distribution function [START_REF] Efron | Better bootstrap confidence intervals[END_REF][START_REF] Efron | An introduction to the bootstrap[END_REF].

random sample Ê = {ε 1 , ε 2 , . . . , ε n } composed of
Accordingly, the resulting distribution E of the horizon modeling error ε will follow a probability law ε ∼ µ E , σ 2 E such as:

         µ E ∈ [ ε * (ν 1 ) , ε * (ν 2 ) ] σ E ∈ [s * (ν 1 ) ε , s * (ν 2 ) ε ] (24)

Model's efficiency and DEM-based error prediction

We have presented the theoretical basis of our horizon model and proposed a method to choose the best compromise between the computation time and the corresponding model's accuracy. Nevertheless, this precision is only related to the ideal DEM-based horizon and not to the real horizon, so the error induced by the not fully accurate DEM representation of the Earth's topography has not been considered yet. As a result, we assess here the global accuracy of the final horizon estimates and compare it to other existing models in order to evaluate our model's performances. Finally, we also propose a method for predicting the horizon's DEM-based error.

Model's efficiency: accuracy vs computation time

The efficiency of any numerical model can be regarded as the ability to describe accurately a phenomenon in the shortest possible time. Therefore, we have firstly assessed the quality of our model by statistically comparing estimates with ground measurements. Then, we have confronted these results and the corresponding computation time to those retrieved using other typical models.

Model's accuracy: comparison with ground measurements and results from other existing models

In order to evaluate the accuracy of our model, we have compared horizon estimates with 10 in situ horizons collected during a measurement campaign achieved in Corsica in 2010; both the geographic distribution and the topographic situation of these sites are presented in Figure 7. The availability of an adapted material was one of the main reasons for choosing Corsica, the other one being the existence of many different types of terrain within a small area, which allowed retrieving a significant sample of measurements without difficulties. In order to measure the different horizons, we have used a topographic mapping device, the Leica Builder 100 theodolite (azimuth accuracy = 9 ; elevation accuracy = 6 ). The operating principle of this accurate device, mainly used by topographers for measuring terrain slope and aspect, is depicted in Figure 8. Furthermore, in order to endorse our model, it was also relevant to compare it to other ones. Thus, we have meanwhile analyzed the results of both the French software Carnaval and the r.horizon add-on of the GRASS GIS open source. Currently released by Sober Software (http://www.sober-software.com), Carnaval is based on the SRTM-3 DEM and only computes single horizons intended to be used with other solar simulation software. On the contrary, r.horizon is closer to our own model since it is dedicated to the integration of terrain effects into solar potential maps computed by the solar radiation model r.sun [START_REF] Ruiz-Arias | A comparative analysis of DEM-based models to estimate the solar radiation in mountainous terrain[END_REF][START_REF] Šúri | A new GIS-based solar radiation model and its application to photovoltaic assessments[END_REF]. Besides, this model allows 2 different computing configurations using a DEM implemented by the user: the first one retrieves shading maps by computing angular elevation for a specific azimuth; the second configuration, and the one we have considered in this work, estimates the full horizon around a given point.

The experiment achieved in 2010 consisted in measuring horizons on 10 different sites across Corsica, representing a full set of 615 points. Moreover, in order to compare the exact same points from both estimated and measured horizons, we have linearly interpolated modeled horizons, at 1 • azimuth step, over measured azimuths. Then, in order to assess the accuracy of each model, we have considered the following well-known indicators: the root mean square error (RMSE), the mean bias error (MBE) and the correlation coefficient (CC). All final results of the comparison are described in Table 1, where each location is defined by the geographic coordinates and the DEM-based altitude of the given observation point.

Globally, these results show higher accuracy of our model, while Carnaval presents the lowest one. Thus, our model and r.horizon present the most homogeneous correlation for all samples with 2 close correlation coefficients (respectively 0.9762 and 0.9761). Meanwhile, our model also exposes a better RMSE (1.555 • against 1.712 • ) and a lower negative bias (-0.416 • against -0.587 • ). Finally, Carnaval presents the less fitted results regardless of the considered indicator: a correlation coefficient of 0.9564, a RMSE of 2.116 • and a bias of -0.705 • . Furthermore, we can note that the increase by 10 of the horizon azimuth resolution of our model only reduces RMSE by 0.1 %, from 1.555 • to 1.547 • , which validates the initial choice of a 1 • azimuth step.

Model's running time

By measuring horizons and comparing them with the estimates of other models, we have demonstrated the accurateness and the relevance of the one developed in this study. Most of all, its theoretical basis will allow accurate estimation of any horizon in the world, as long as the SRTM-based DEM exists. Nonetheless, as previously depicted, this model has been primarily designed for considering shading effects within satellite-based solar radiation maps [START_REF] Bosch | Solar resources estimation combining digital terrain models and satellite images techniques[END_REF][START_REF] Haurant | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica[END_REF][START_REF] Pillot | The first disaggregated solar atlas of Djibouti: A decisionmaking tool for solar systems integration in the energy scheme[END_REF][START_REF] Ruiz-Arias | Spatial disaggregation of satellite-derived irradiance using a high-resolution digital elevation model[END_REF]. Therefore, as with the previous optimization method, used to reduce the area size required to compute all horizons of a given region, everything has been made to minimize computation time.

Firstly, as our model is written with the Matlab programming language, it takes benefit from the software specific matrix calculations, which saves time for handling of large arrays. Thus, on the same machine and for the same viewing distance

Model

Running time (s)

This study ≈ 0.5 r.horizon 1 ≈ 2.0 Carnaval2 2.0

1 Given by the GRASS GIS command window.

2 Estimated by Sober Software.

(50 km), Table 2 shows our model is, per thread and per horizon, 4 times faster or more than r.horizon and Carnaval 2 . Besides that, our model is also directly incorporated into the process implemented for considering topography within radiation maps [START_REF] Pillot | Planification de l'électrification rurale décentralisée en Afrique subsaharienne à l'aide de sources renouvelables d'énergie : le cas de l'énergie photovoltaïque en République de Djibouti[END_REF][START_REF] Pillot | The first disaggregated solar atlas of Djibouti: A decisionmaking tool for solar systems integration in the energy scheme[END_REF], which avoids extra time consumption during data transmission with external models. Finally, our model's source code used to estimate all horizons of a given region is parallelized, which reduces running time as much as the number of available threads within the computer's processor.

DEM-based error prediction

A more detailed reading of Table 1 reveals that, regardless of the given model, error and bias also dramatically vary from one horizon to another. Indeed, some sites present very good estimates (Bastia, Solenzara or Calvi) while some others show both high RMSE and negative bias (Tavignano, Venaco or Castellare). If the original DEM error obviously affects the final quality of the horizon estimates, it is also interesting to know whether this error depends on some specific topographic parameters, which would have influenced the original SRTM precision. In that case, the resulting horizon error would also depend on some specific features; consequently, we propose a method for predicting this error, or at least to detect if the error of any given horizon will probably be significant.

Leaving aside possible measurement errors (angular height, azimuth or geographic coordinates of the observation point), the main idea is to connect the horizon error to the DEM error. To perform this, it is firstly necessary to know the kind of error existing within the SRTM DEM. It is possible to enumerate 2 main errors: vertical and georeferencing errors, respectively achieved on the altitude and the geographic coordinates of a DEM point [START_REF] Farr | The shuttle radar topography mission[END_REF][START_REF] Rodríguez | A global assessment of the SRTM performance[END_REF]. In this study, we have only considered the vertical error, easier to implement and the one which mainly induces the angular elevation error.

In order to predict the error generated by the DEM quality, it is firstly relevant to express the RMSE of a given horizon consisting of n points:

RMSE = n i=1 α i -α i 2 n ( 25
)
With α i the ith estimated angular elevation and α i the ith measured one. Let us then focus on the ith point of this horizon, of angular height α i ; it is known this point has a vertical error on its DEM altitude h i , with respect to the true altitude h i , which is propagated on the angular error α i -α i . So, if it is possible to predict the vertical error h ih i , it is also possible to assess the resulting angular error, by estimating α i from h i using relation [START_REF] Hofierka | Open source solar radiation tools for environmental and renewable energy applications[END_REF]. Essentially, we have thus been led by the main idea of looking for an empirical expression, allowing us to merely define the vertical error as a function of SRTM DEM-derived parameters.

Validation of the SRTM has shown the vertical error can be regarded as a random variable following an almost normal probability law, centered on 0 and having a standard deviation specific to each continent [START_REF] Farr | The shuttle radar topography mission[END_REF][START_REF] Rodríguez | A global assessment of the SRTM performance[END_REF]. Meanwhile, the Miliaresis and Paraschou ( 2005) and [START_REF] Gorokhovich | Accuracy assessment of the processed SRTM-based elevation data by CGIAR using field data from USA and thailand and its relation to the terrain characteristics[END_REF] studies have highlighted the error correlation with topographic features of a given terrain, typically both the slope and aspect of a DEM pixel. In this study, we have only considered the slope S , since the aspect is less significant but also more difficult to model as it is highly dependent on the study region. The slope is a function of the altitude gradient to the north (y) and to the east (x) [START_REF] Zhou | Analysis of errors of derived slope and aspect related to DEM data properties[END_REF]:

S = arctan          ∂h ∂x 2 + ∂h ∂y 2          (26)
Several methods exist for calculating S , among which the one developed by Fleming and Hoffer presents a good precision [START_REF] Jones | A comparison of algorithms used to compute hill slope as a property of the DEM[END_REF]. It is based on the formulation of the partial derivatives of the height h [START_REF] Zhou | Analysis of errors of derived slope and aspect related to DEM data properties[END_REF]:

∂h ∂x = h E -h W 2R x DEM ; ∂h ∂y = h N -h S 2R y DEM (27)
Where h E , h W , h N , h S are respectively heights of eastern, western, northern and southern pixels contiguous to the pixel of interest, R x DEM the DEM metric resolution along x axis, and R y DEM the resolution along y axis. Table 3 summarize, for each site, the main statistical parameters (mean, standard deviation and range) of the horizon slope, that is the slope of all DEM pixels composing the horizon.

Both [START_REF] Wiley | Vertical accuracy of the SRTM DTED level 1 of Crete[END_REF] and [START_REF] Gorokhovich | Accuracy assessment of the processed SRTM-based elevation data by CGIAR using field data from USA and thailand and its relation to the terrain characteristics[END_REF] studies have emphasized the increase of the mean vertical error with the slope S (overestimation of the DEM altitude compared to the true altitude). [START_REF] Wiley | Vertical accuracy of the SRTM DTED level 1 of Crete[END_REF] especially have shown that, even though the error distribution was translated to the right, this distribution was still close to the one for low slopes (S < 2 • ) with about the same flattening and asymmetry. Considering these assumptions, it is then possible to define the vertical error as: Where r is a normally distributed random variable centered on 0 such as r ∼ N(0, σ), resulting of the SRTM validation [START_REF] Rodríguez | A global assessment of the SRTM performance[END_REF], and f (S ) any function of the pixel's slope S . From the Miliaresis and Paraschou ( 2005) study, we can then consider, in first approximation, the error linearly increases with the slope. Finally, the expected value of the vertical error hh is given by:

h -h = r + f (S ) (28)
E h -h = E [r] + βS = βS ( 29
)
Where β is a strictly positive constant (distribution translated to the right). Getting h from relation (29), it is then possible to predict the value of the corresponding measured angular elevation α, using equation [START_REF] Hofierka | Open source solar radiation tools for environmental and renewable energy applications[END_REF]. In order to evaluate the constant β, we have minimized the quadratic error between the measured RMSE of the full sample in Table 1 and the RMSE estimated with this method. Carnaval and r.horizon presenting different errors, the constant β is therefore different from one model to another, and is respectively equal to 0.357, 0.398 and 0.494 for our model, r.horizon and Carnaval. In addition, the correlation between the prediction and the measure for each model is given by the correlation coefficient. It is equal, for a 1 % significance level, to 0.9331, 0.9611 and 0.9800 for our model, r.horizon and Carnaval respectively. Finally, the consistency of the prediction method is depicted in Figure 9, with on the one hand the estimates vs measurements scatter plot, and on the other hand the comparison between 2 horizons presenting low (Corte) and high error (Tavignano).

Some aspects of this analysis must be emphasized. First of all, we have not taken into consideration the vertical error of the horizon observation point, although it also interferes in the global angular error. Indeed, if the many points composing the horizon allow smoothing the approximate formula (29), the possible bias introduced by the observation point would, on the contrary, impact all the points, and so could strongly distort the final estimate. This point leads to another aspect of this study: the greater the number of points composing the horizon, the more accurate the quadratic error prediction will be. Furthermore, vegetation is also a subject of concern for the error estimation since depending on its type and density, the wavelength of the C band might penetrate into the canopy [START_REF] Farr | The shuttle radar topography mission[END_REF]. Accordingly, in vegetated areas, SRTM DEM elevation is located between the ground and the top of the canopy, at a distance which varies significantly with the SRTM DEM variability [START_REF] Carabajal | SRTM C-band and ICESat laser altimetry elevation comparisons as a function of tree cover and relief[END_REF]. As the shading effects are influenced by the canopy, this might lead to some extra underestimation, mainly when the horizon is close. Finally, linear approximation between vertical error and slope as well as the decision to leave aside DEM aspect and georeferencing may also explain the shift between measured and predicted values. Even so, this method shows quite accurate estimates, and can further be generalized to all the SRTM-3 DEM, by mainly regarding the predicted RMSE as a sensitivity coefficient allowing the horizon's level of quality to be assessed.

Error propagation to solar radiation

As depicted in previous studies, the integration of shading effects improves both the accuracy and spatial information provided by satellite-derived solar radiation maps [START_REF] Haurant | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica[END_REF][START_REF] Pillot | Planification de l'électrification rurale décentralisée en Afrique subsaharienne à l'aide de sources renouvelables d'énergie : le cas de l'énergie photovoltaïque en République de Djibouti[END_REF][START_REF] Ruiz-Arias | Spatial disaggregation of satellite-derived irradiance using a high-resolution digital elevation model[END_REF]. However, as described above, the DEM-based horizon computation results itself in some error, which necessarily propagates to the solar irradiation. In order to assess how this error affects the final irradiation estimation, we have used the disaggregation methodology we have developed in another study [START_REF] Pillot | Planification de l'électrification rurale décentralisée en Afrique subsaharienne à l'aide de sources renouvelables d'énergie : le cas de l'énergie photovoltaïque en République de Djibouti[END_REF], coupled with the ESRA clear-sky model developed by [START_REF] Rigollier | On the clear sky model of the ESRA -european solar radiation atlas -with respect to the heliosat method[END_REF] and corrected in altitude by [START_REF] Geiger | A web service for controlling the quality of measurements of global solar irradiation[END_REF]. This methodology includes shading effects into the global solar radiation by applying horizon-derived factors to its diffuse isotropic, diffuse circumsolar and beam components.

Besides, depending on both the sun path and the sky cloudiness, the final influence of the horizon on the solar radiation can be significantly different. In order to apprehend this double dependence, we have estimated the clear-sky daily irradiation throughout the year, and analyzed the influence of the solar radiation's diffuse component by acting on the key parameter of the ESRA model: the Linke turbidity factor T L [START_REF] Geiger | A web service for controlling the quality of measurements of global solar irradiation[END_REF][START_REF] Rigollier | On the clear sky model of the ESRA -european solar radiation atlas -with respect to the heliosat method[END_REF]. This coefficient usually varies from 1 to 10, T L = 2 corresponding to a very clear atmosphere, and T L = 8 corresponding to a polluted atmosphere [START_REF] Eltbaakh | Issues concerning atmospheric turbidity indices[END_REF]. In order to evaluate the solar irradiation loss potentially caused by the horizon, we have determined the daily clear-sky irradiation reaching the Tavignano site, with and without shading effects computed from the measured horizon. Figure 10 depicts the final relative loss according to the day of the year and the Linke turbidity factor. As expected [START_REF] Quaschning | Irradiance calculation on shaded surfaces[END_REF], shading effects affect more the beam component than the diffuse radiation: the irradiation blocked by the horizon decreases when T L increases. When the sun is high however, between days 100 and 250, the loss remains low and similar (< 10%) regardless of the Linke coefficient. In that case, the most significant loss (> 50%) appears at the end of the autumn and the beginning of the winter, as a consequence of the low sun path.

Solar irradiation loss

Horizon error propagation

In order to assess how the RMSE of 1.555 • (Table 1) ends up affecting the solar estimation, as well as the evolution of this influence, we have statistically analyzed the daily irradiation with shading effects reaching the 10 sites of Table 1, using either the horizon measurements or the estimates from our model. To do this, we have computed the relative root mean square error (RRMSE) between both the results, according to the day of the year and the Linke turbidity factor, depicted in Figure 11.

Again, as expected, the beam radiation is more impacted by the horizon shadowing than the diffuse radiation, as the final error decreases with the rise of the Linke coefficient. Also, the RRMSE is significantly variable over the year, compared to the smooth tendency of the solar radiation loss; it can probably be explained by the original variability of the sign of the bias between the estimate and the measured horizon elevation depending on the azimuth. Essentially, in this study, depending on the day and the diffuse fraction, the horizon mismatching results in a final error ranging from approximately 0.5% (high sun path) to 5.5% (low sun path) of the daily solar irradiation. 1, with respect to the day of the year and the Linke turbidity factor.

Conclusion

Many models now retrieve solar radiation from satellite data, but are still limited to an atmospheric characterization without any consideration of the interaction between radiation and Earth's surface. Therefore, in order to improve the accuracy of these satellite-based estimates, some post-processing methods have been proposed for integrating terrain-based effects into final solar irradiance or irradiation maps using DEMs. The main process of all these methods consists in estimating the obstructed horizon around every DEM pixel, in order to evaluate direct and diffuse shadowing. Nevertheless, depending on the size of the study area, computation can also be really time consuming. That is why we have developed and validated a new efficient horizon model, combined with an optimization method and a DEM-based error prediction study.

We have developed this model using the global DEM retrieved from the SRTM, as it is currently the most consistent and accurate available DEM. We have furthermore only based this model on the mathematical statements used to represent the DEM, avoiding any region-specific empirical assumption so as to allow the model to be used at any world location with the same expected accuracy. Then, we have assessed its efficiency by comparing its precision and computation time with those of other models. Precision has been retrieved using a statistical comparison between estimates and in situ horizons measured with a theodolite. Finally, we have demonstrated the implementation of our model was relevant, since it is both more accurate (RMSE of 1.555 • against 1.712 • or more for other presented models) and less time consuming (at least 4 times faster than other models).

In addition, we have also proposed an optimization method for reducing horizon computation time in the case of large areas, as well as an empirical approach for predicting the DEM-based error on final estimates. The first one takes into consideration the Earth's roundness in order to limit the maximum zone required to fully describe any horizon of a given region. Hence, it allows knowingly decreasing of modeling accuracy with running time, and thus get the best compromise between them depending on the purpose of the regarded study. The second evaluates the impact of the SRTM DEM original error on the final horizon precision by correlating it to the DEM slope, and has been validated using the in situ horizon sample collected in Corsica.

Finally, we shall note that, in this study, we have not considered the distortion effects of the atmospheric refraction onto the apparent position of objects on Earth, and so on the resulting angular elevation, known as terrestrial refraction [START_REF] De Graaff Hunter | Formulae for atmospheric refraction and their application to terrestrial refraction and geodesy[END_REF]. This approach seems relevant since, in mapping applications, the additional computation time required for modeling this phenomenon, and its related parameters such as object distance or temperature and pressure conditions, would not necessarily be worth the final accuracy gain. Nevertheless, future research could still improve the accuracy of our current model by focusing on this part.

∆

  WGS 84-EGM96 geoid undulation with respect to the ellipsoid (m) ∆λ longitude width of the DEM tile used in horizon computation (rad) ∆ϕ latitude height of the DEM tile used in horizon computation (rad

Figure 1 :

 1 Figure 1: WGS84 ellipsoid of revolution used as the reference geographic coordinate system for the SRTM DEM.

Figure 2

 2 depicts the horizon observed from A(ϕ A , λ A , h A ) and including B(ϕ B , λ B , h B ), defined in the horizontal plane P with respect to A by its azimuth γ AB and its elevation α AB . DEM-based horizon modeling consists therefore in computing elevation α and azimuth γ of all the points of the DEM around the observation point, and then retrieving the maximum α for every azimuth between 0 and 360 • .

Figure 2 :

 2 Figure 2: Obstructed horizon, induced by the terrain topography, represented in the horizontal coordinates (γ, α) of the plane P centered on A.

Figure 3 :

 3 Figure 3: Sectional view of the ellipsoid surface in the AB direction.

Figure 4 :

 4 Figure4: Azimuth γ of a line drawn on the ellipsoid and calculated from meridian (ds m ) and parallel (ds p ) line elements[START_REF] Bosser | Géométrie de l'ellipsoïde[END_REF].

  the bootstrap concept. It consists in re-sampling F, i.e. in drawing with replacement B n-sized random samples F * (b) in F, with b = 1, 2, . . . , B, called bootstrap samples. Then, the distribution of the B boostrap sample estimates θ * (1) , θ * (2) , . . . , θ * (B) around θ, such as θ * (b) = t( F * (b) ), is used to calculate the standard error and the confidence interval of the estimate θ (Efron[START_REF] Efron | An introduction to the bootstrap[END_REF][START_REF] Singh | Bootstrap method[END_REF].

Figure 6 :

 6 Figure 6: Example of the box plot used to optimize horizon calculation in the Republic of Djibouti, where model's error is plotted against mean running time

  a larger set of observations (for instance n = 200) is implemented. From this empirical population, B = n 2 bootstrap samples Ê * (b) = ε * 1 , ε * 2 , . . . , ε * n , where b = 1, 2, . . . , B, are drawn with replacement (Singh and Xie, 2010). Thereafter, both the mean ε * (b) and the standard deviation s * (b) ε of each sample are calculated. Their distribution allows assessing the accuracy of the mean ε and the standard deviation s ε of the empirical population, which infer the parameters µ E and σ E of the global population.

Figure 7 :

 7 Figure 7: Topographic SRTM-DEM map of Corsica island (meters) with location of the measurement sites.

Table 1 :Figure 8 :

 18 Figure8: Theodolite operating principle. This device allows both the azimuth γ and the angular height α of a given point to be accurately measured.

igure 9 :

 9 Predicting SRTM-3 DEM-based horizon error. Above is the scatter plot of predicted against measured RMSE values for the 10 horizons depicted in Table 1. Below is an example of horizon measurements and estimates over Corte and Tavignano, respectively for a low and a high error.

Figure 10 :

 10 Figure 10: Relative irradiation loss (%) caused by the horizon at Tavignano site, with respect to the day of the year and the Linke turbidity factor.

Figure 11 :

 11 Figure 11: RRMSE (%) between solar irradiation with shading effects computed from horizon measurements and estimates for all sites ofTable 1, with

Table 2 :

 2 Approximate horizon per-thread computation time for a 50 km view distance around observation point.

Table 3 :

 3 Main characteristics of the horizon slope.

	Site	Slope ( • )
		Mean Std Range
	Ajaccio	10.1 7.4 0.0 -33.5
	Calvi	11.1 9.4 0.6 -40.5
	Lozari	8.7 8.7 0.0 -46.9
	Bastia	6.4 7.8 0.4 -35.3
	Solenzara	7.9 8.5 1.2 -39.3
	Tavignano	21.7 5.0 9.3 -31.1
	Castellare	16.0 6.2 4.0 -35.8
	Venaco	14.4 6.9 0.5 -35.8
	Corte	16.0 8.7 1.4 -39.8
	Soveria	17.3 7.3 2.0 -38.0

At the time this study was driven, worldwide SRTM-1 DEM was not available yet. It has been meanwhile released and is now available from the NASA/USGS website (https://lta.cr.usgs.gov/SRTM1Arc).

So as to get an idea of the computation time over areas, at 2 s/horizon, it would spend about

days, 2000 days or 250 years to retrieve all horizons from Corsica, France or Africa respectively.

Acknowledgments

We are grateful to Dr Michael Mercier-Finidori and the Corsica University Institute of Technology for providing the theodolite we used in the experimental part of this work.

Arvizu, D., Balaya, P., Cabeza, L. F., Hollands, T., Jäger-Waldau, A., Kondo, M., Konseibo, C., Meleshko, V., Stein, W., Tamaura, Y., Xu, H., Zilles, R., 2011. Direct solar energy. In: IPCC special report on renewable energy sources and climate change mitigation. Cambridge University Press, Cambridge, U.K. and New York, NY, USA, Ch. 3.

Ben Ahmed, H., Multon, B., Thiaux, Y., 2011. Consommation d'énergie, ressources énergétiques et place de l'électricité.

Techniques de l'ingénieur, Paris, France. Bermejo, P., Pino, F. J., Rosa, F., Aug. 2010. Solar absorption cooling plant in Seville. Solar Energy 84 (8), 1503-1512. Beyer, H. G., Costanzo, C., Heinemann, D., 1996. Modifications of the heliosat procedure for irradiance estimates from satellite images. Solar Energy 56 (3), 207 -212.


Jeferson A L De Souza
email: jeferson.souza@usp.br

Diolino J Santos Fo
email: diolinos@usp.br

Reinaldo Squillante
email: reinaldo.squillante@usp.br

Fabricio Junqueira

Paulo E Miyagi
email: pemiyagi@usp.br

Mitigation Control of Faults in Critical Production Systems

Keywords: Critical Systems, Mitigation Control System, Safety Instrumented System, Fuzzy Logic

The inherent complexity of critical production systems, coupled with policies to preserve people´s safety and health, environmental management, and the facilities themselves, and stricter laws regarding the occurrence of accidents, are the motivation to the design of Safety Control Systems that leads the mitigation functionality. According to experts, the concept of Safety Instrumented Systems (SIS) is a solution to these types of issues. They strongly recommend layers of risk reduction based on hierarchical control systems in order to manage risks, preventing or mitigating faults, a n d to lead the process to a safe state. Additionally some of the safety standards such as IEC 61508, IEC 61511, among others, guide different activities related Safety Life Cycle design of SIS. The IEC 61508 suggests layers of critical fault prevention and critical fault mitigation. In the context of mitigation control system, the standard provides a recommendation of activities to mitigate critical faults, by proposing control levels of mitigation. This paper proposes a method to implement the mitigation layer based on the risk analysis of the plant and the consequences of faults of its critical components. The control architecture, based on distributed and hierarchical control systems in a collaborative way, will make use of the techniques of risk analysis raised and mitigation actions, based on the knowledge of an expert, implemented by fuzzy logic.

Introduction

In this first decade of the century XXI many studies have indicated that automation processes are undergoing transformations that have been strongly influenced by the advance of technology and computing resources, becoming increasingly complex due to their dynamic and needed to address issues such as global market competitive production and technology used, among other factors [START_REF] Chen | Design and high-level synthesis of hybrid controller[END_REF], [START_REF] Santos Filho | Aspectos do Projeto de Sistemas Produtivos[END_REF], [START_REF] Wu | Service-oriented communication architecture for automated manufacturing system integration[END_REF]. Given this new scenario, industrial processes and their control are becoming more complex. Additionally, organizations have focused on policies to achieve and to demonstrate people's safety and health, environmental management system, and controlling risks. In this context, any industrial system, as modern and innovative as can be, could be considered to pose a serious risk to people's health, the environment and equipment [START_REF] Sallak | A fuzzy probabilistic approach for determining safety integrity level[END_REF].

Although many studies have been presented for diagnosis and treatment of faults, a review of fault-tolerant reconfigurable control system can be found in [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF], b u t a ccidents still occur. These issues are fully justified because there is no zero risk in process industries since: (i) physical devices do not have zero risk of failure [START_REF] Summers | Common cause and common sense, designing failure out of your safety instrumented systems (SIS)[END_REF], (ii) human operators do not have zero risk of error and (iii) there is no computational software project developed that can predict all the possibilities [START_REF] Miyagi | Controle Programável -Fundamentos do controle de sistemas a eventos discretos[END_REF].

According to experts, the concepts of safety instrumented systems (SIS) is a solution to these types of issues and strongly recommend layers of risk reduction based on control systems organized hierarchically in order to manage risks by either preventing or mitigating faults, and to bring the process to a safe state. In this sense, some safety standards such as IEC 61508 [START_REF] Iec | Functional safety of electrical / electronic / programmable electronic safety-related systems (IEC 61508)[END_REF], IEC 61511 [START_REF] Iec | Functional safety -safety instrumented systems for the process industry sector -part 1 (IEC 61511)[END_REF] among others, guide activities related with a SIS Safety Life Cycle (SLC), such as design, installation, operation, maintenance, tests and others [START_REF] Lundteigen | Architectural constraints in IEC 61508: Do they have the intended effect? Reliability Engineering and System Safety[END_REF], [START_REF] Bell | Introduction to IEC 61508[END_REF].

According to IEC 61508, the term "fault" is defined as an abnormal condition that can cause a reduction or loss of the ability of a functional unit, and is defined two layers of SIS: the prevention layer and the mitigation layer. Recently, [START_REF] Squillante | Mathematical method for modeling and validating of safety instrumented system designed according to IEC 61508 and IEC 61511[END_REF] proposes the implementation of a SIS prevention layer.

This work is initially proposed a systematic for modeling and validating layer of mitigation control within SIS. This approach considers the cause of the fault, its severity and its consequence for the system, through the application of risk analysis techniques such as Failure Modes and Event Analysis -FMEA, Fault Tree Analysis -FTA [START_REF] Modarres | Reliability Engineering and Risk Analysis: a practical guide[END_REF], and the What-If technique [START_REF] Souza | O treinamento industrial e a gerência de riscos -uma proposta de instrução programada[END_REF], based on a database of occurrence of faults or on knowledge of an expert or operator. The effects and the consequences of the occurrence of a critical fault, listed on the risk analyses study, are monitored and treated by the SIS sensors and actuators, respectively, independently of the BPCS devices, as predicts the IEC 61508 [START_REF] Squillante | Safety in Supervisory Control for Critical Systems[END_REF], [START_REF] Cavalheiro | Design of supervisory control system for ventricular assist device[END_REF]. The effect of every critical fault, or safety instrumented function (SIF), results in mitigation actions, determined by the What-If technique yet implemented.

Fuzzy logic is utilized for the generation of the control algorithm. It has the advantage of not using differential equations or complex mathematical models for determining the dynamic behavior of the system [START_REF] Popa | A PI-Fuzzy controller designated for industrial motor control applications[END_REF], and can therefore use the proposed mitigation actions, which in turn were the result of applying the techniques of risk analysis, for the determination of fuzzy control algorithms. Another advantage of fuzzy logic is the analysis of the parameter time derivative [START_REF] Legaspe | Open source fuzzy controller for programmable controllers[END_REF], thus contributing to an anticipatory action of the proposed mitigation layer. The generation of control codes for programmable logic controllers (PLC) can be made based on IEC 61131-7 [START_REF] Iec | Programmable Controllers IEC 61131-7: Fuzzy Control programming[END_REF], which deals with the conversion of the generated fuzzy logic algorithm to conventional PLC languages, based on IEC 61131-3. [START_REF] Iec | Programmable controllers IEC 61131-part 3: Programming languages[END_REF].

Relationship to Collective Awareness Systems

Recent political awareness with focus on sustainability and recycling, the use of resources and raw materials from renewable resources, along with the practices of waste management and emission control of pollutants, coupled with more rigid and punitive laws to production systems that do not meet the new regulations, results in new control systems in manufacturing plants.

In this context, a mitigation control layer is needed because, in addition to the market requirements of competitiveness, such policies result in a choosing collective awareness of choice by consumers for companies committed to focus on the environment and sustainability. Besides the increase in the complexity of the control systems for these new requirements, the proposed mitigation layer, based on IEC 61508, is precisely the preservation of men and the environment, just required for productive systems that wishing to adapt to these new practices of collective awareness recently observed.

3

Proposal of Layer of Mitigating Control System

Description of the Proposed Method

The proposed method is summarized by the following five steps, described on sections 3.1.1 to 3.1.5 below:

Determination of the Critical Elements

To determine the critical elements of the process under study we utilize the risk analysis techniques FTA, FMEA and What-If. The FMEA, to associate a severity level to the occurrence of fault of a component indicates which components must be monitored in the mitigation layer. Faulted components that pose risks to operators, the environment and equipment, besides violating the legislation are classified to maximum severity. Furthermore, components which fault under no danger considerable not part of our analysis.

Because the FMEA to be centered on the component, combination of faults and a possible domino effect over other components may be analyzed by the FTA in conjunction to What-If technique. It is possible, to determine the how and the why of the fault, therefore rendering a more comprehensive study.

Detection of Effects Caused by the Occurrence of Faults of the Critical Elements

Each effect arising from the fault of a critical component must be monitored by a specific sensor for its fault mode. According to IEC 61508, such sensors must be independent of the BPCS. To avoid spurious faults and reading errors, it is recommended to use redundant architectures [START_REF] Bell | Introduction to IEC 61508[END_REF], such as the criteria voting 2oo3 (two of three).

Mitigation Actions of the Effects of Faults of Critical Elements

For each effect of a critical fault, detected by the SIS mitigation sensors, a mitigation action must be implemented by SIS mitigation actuators, controlled by the SIS mitigation control layer, aiming to preserve people, environment and equipment.

To determine de mitigation actions will make use of What-If technique, based on human knowledge and records of occurrence of faults, its effects and the actions proposed to mitigate its effects.

Some mitigation actions can be matched to faults occur in different components, but not necessarily input signals from SIS mitigation sensors are the same. That is, for different input signals, from different sensors, mitigation actions may be the same. In this step, besides determination of the mitigation sensors can arrive at the conclusion that the sensors would be the same prevention. In this case, it is recommended doubling of signals from sensors for prevention PLC we use in our mitigation.

After this study and compilation of mitigation actions, will determine which actuators required for each mitigation action.

Construction of Models for Implementation ofControl Algorithms

In this step will be used results of the What-If technique already implemented in section 3.1.3 to determine the level or percentage of the measured variable values for activation layer of prevention and / or mitigation using the absolute value of the measured variable and its temporal variation or derivative of the measured parameter. The results of this study will form the basis of fuzzy algorithms for mitigation control layer.

Control Codes Generation based on IEC 61131-7

For each mitigation action determined by the fuzzy control algorithm, the next step is to convert the generated control algorithm for a language of IEC 61131-3 to implement in the Safety PLC for mitigation.

The IEC 61131-7 deals with the implementation of fuzzy algorithms in FCL (Fuzzy Control Language), based on IEC 61131-3 [24] ST (Structured Text) for the implementation in conventional PLCs.

Example of Application

To illustrate the method proposed, an application example for critical faults to be mitigated by SIS Mitigation layer in a natural gas compression station is presented. Natural gas is a mixture of highly flammable hydrocarbons. To be extracted from the environment must be pressurized in compressor stations to its carriage due to consumer centers.

Process Description

The natural gas station has one or more natural gas supply lines, called suction, from a gas pipeline which transports this natural gas. At the station entrance, natural gas goes through filters equipment before being compressed by the turbo compressor machine. A portion of this gas is directed to the utility unit. The utility unit accounts for controlling the gas temperature and pressure for use in the compression station, such as fuel gas for the turbo-compressor machine, gas heaters and gas power generators. After the natural gas is compressed by turbo compressor machine, it is sent back to the gas pipeline through discharge lines, called headers. The PI&D of a turbo compressor uni tis shown in Fig. 1.

Application of the Proposed Method

We apply the proposed method, based on the SIS prevention layer proposed by [START_REF] Bell | Introduction to IEC 61508[END_REF] for the case of the compressor gas compression station. A more elaborate study should be done, considering all critical components indicated by the application of FMEA and FTA techniques. This work presents an example for a system component, in order to exemplify the application of the proposed method.

Determination of the Critical Elements

Applying the FMEA technique can be seen that the compressors are critical to effectively our system, because they operate under high temperature, pressure and speed, in addition to use as fuel the compression fluid itself, which is natural gas, just explosive. A fault in this equipment certainly put under unacceptable risk operators, the environment and the equipment itself, besides violating government standards for safety. Hence its severity is maximum, and must be entered in our mitigation layer. Table 1 illustrates a FMEA for compressor, and Fig. 2 the FTA for the top event "High Temperature Lubricants Shaft". Both FMEA and FTA found that an effect of occurrence of fault in the compressor is to increase the temperature on the cooling fluid turbine shaft, being able to have also an increase in temperature of the working fluid in the discharge line. We will perform our study on the mitigation system as a function of monitoring the temperature parameter for this component. Other effects can be measured as changes in discharge pressure coming from a lower performance of the compressor operating under fault.

Detection of Effects Caused by the Occurrence of Faults of the Critical Elements

For the effects of faults listed in the previous step, we have temperature sensors coolant axis of compressors, independently of the BPCS. Such sensors will be designated TAT 211 -Temperature Axis Turbine -for each unit present in the natural gas station. So we have the TAT 211 A, TAT 211 B, TAT 211 C and TAT 211 D as input signals our mitigation PLC. Again, a redundant architecture of these sensors as well as the implementation of algorithms for detecting spurious faults [START_REF] Bell | Introduction to IEC 61508[END_REF] should be implemented.

Mitigation Actions of the Effects of Faults of Critical Elements

To mitigate the effects caused by the occurrence of a fault in the compressor, beside the action of shutdown from the prevention layer, suggested action to mitigate the effects is the forced cooling of the compressor, if preventive layer is not sufficient or if the temporal variation of temperature proves too high. Will be used both carbon dioxide cylinders large that are already installed in natural gas station, and have the purpose of fire combat if an outbreak of fire. The release of carbon dioxide is currently done manually, through the action of fire brigade teams, specially trained for this purpose. The proposal would be the installation of pipelines leaving the cylinders to compressors with proportional valves connected to the outputs of mitigation´s Safety PLC. As the intensity of mitigation action, the valve would release the carbon dioxide in the same proportion.

Construction of Models for Implementation of Control Algorithms

From mitigation proposals have the construction of the control algorithms implemented by fuzzy logic, from the What-If technique already implemented, based on the expertise of a specialist. To illustrate the algorithm, the expert reports that 150% of the temperature set point would be unacceptable to the turbo compressor. So we adopted a range of 110% to 130% for the prevention layer. Above 120% mitigation layer already comes into operation in a proportional action. Note that the temporal variation in temperature is part of the algorithm´s control input. Fig. 3 illustrates the proposed model for temperature: The rules of the fuzzy algorithm, according to What-If technique are as follows in Table 2. The output signal, or the proportional action of mitigation, here designated by proportional valve opening, can be seen by the generated surface on Fig. 6.

We can see from the graphs of anticipatory mitigation action due to the temporary increase of the measured variable. This results in better efficiency of the system, thus contributing to a further reduction of the inherent process risk.

Control codes generation based on IEC 61131-7

From the algorithms based on fuzzy logic implementation has the control codes to Safety PLC for mitigation, considering the anticipatory model, as shown in Fig. 7.

Conclusions

A method for the implementation of mitigation layer in critical industrial systems was proposed, based on the IEC 61508 and IEC 61511 standards, which recommend layers of risk reduction based on cooperative and hierarchical control prevention and mitigation of critical faults. Based on the results of applying the risk analysis techniques can be evaluated, due to the effects of their faults, what the critical elements present in the process. Based on the knowledge of an expert and making use of the What-If technique already deployed, implement corresponding mitigation actions using fuzzy logic, becoming such an algorithm in industrial PLCs languages based on IEC 61131-7. This layer proposal, coupled with the prevention layer, contributes to reduce the inherent risk in the process and adding to the temporal analysis of the variable associated with the effect of a critical component fault results in anticipatory mitigation action, resulting in a higher process risk reduction.

A refinement of this method can be accomplished by inserting a larger set of terms for de derivative membership function, such as PS (Positive Short), PM (Positive Medium) and PH (Positive High) and adopting the same procedure for negative derivative. Intermediate values of the actuator, eg 30% may be associated with these new values, which will surely determine new fuzzy rules in the algorithm. Other mitigation actions can be proposed, and this model must be implemented for the other critical elements of plant. Such elements may have other parameters that indicate the fault component and also other mitigating actions.

Fig. 1 .

 1 Fig. 1. P&ID of a compressor unit of the gas compression station.

Fig. 2 .

 2 Fig. 2. Suggested FTA for the top event "High Temperature Lubricants Shaft".

Fig. 3 .

 3 Fig. 3. Fuzzy Membership functions for temperature

Fig. 4 .

 4 Fig. 4. & Fig. 5. Membership functions to temperature derivative and percentage of valve opening.

Fig. 6 .

 6 Fig. 6. Surface generated by the fuzzy algorithm by the fuzzy rules defined for the mitigation model.

Fig. 7 .

 7 Fig. 7. Control code generated, implemented in FCL (Fuzzy Control Language) according to IEC 61131-7.

Table 1 .

 1 Proposed FMEA for temperature increase of the lubricating oil of shaft bearing compressor

Table 2 .

 2 Fuzzy rules for mitigation layer.

	IF	TEMPERATURE	OPERATOR	DTEMPERATURE	THEN	VALVE
	1	MITIGATE				HIGH -100% OPEN
	2	BASIC CONTROL				ZERO -CLOSED
	3	PREVENT	AND	P		MEDIUM -50%
	4	MITIGATE1	AND	N		MEDIUM -50%

Acknowledgments. The authors would like to thank the Brazilian governmental agencies CNPq, FAPESP, and CAPES for their financial support to this work.