L F Salazar Betancourt
email: luis-fernando.salazarbetancourt@mines-paristech.fr

P Laure

L Silva

M Sager

Thermo-rheological-kinetical Study Of Compression Molding Of Fibre-reinforced Composites

To improve the modeling of fiber reinforced composites, we present in this work numerical methods able to compute both fiber-reinforced composites deformation in squeeze flow and thermal-kinetic evolution. The rheology is given by an homogeneous orthotropic model for fiber composites which describes the anisotropy of the in-plane fiber. The thermics is then extended accounting for the reaction here formulated by the Bailleul's model. Both physics are related since the kinetic evolution as well as the temperature profile modify the rheology of the composites, giving raise to the thermo-rheological-kinetical coupling by means of the viscosity temperature dependence. A study case is presented, where the mold temperature is set to 150 • C with a composite sample at 40 • C. Thermal transfer begins as well as sample compression at constant speed. We present the evolution of the reaction, temperature and viscosity at the core and the surface. Reaction in the core of the material is much quicker than in the surface. Which means that a mapping of viscosity values is presented during the reaction modifying the mechanical response.

INTRODUCTION

Fiber-reinfroced Composites materials represent an interesting approach for lighting automotive structural pieces. The mechanical modeling of such materials has started with the work of Barone et al [1] improved with the study performed by Kotsikos et al [2]. Nevertheless, it was until the work of Le Corre [3] when its anisotropic modeling began. Later, Dumont et al in [4] proposed a model that accounts for the anisotropic behavior exhibited by reinforced materials. By means of structural tensors a full formulation is presented allowing its direct numerical implementation as shown in [START_REF] Salazar Betancourt | Numerical implementation of a rheology model for fiber-reinforced composite and viscous layer approach for friction study[END_REF]. However, current models taking into account the thermo-kinetic nature of such composites remains a complementary key towards the thermo-mechanical modeling of reinforced parts. First studies were accomplished by Lee in [6] and by Maazouza [7] using a model based on the radical polymerisation mechanisms. Currently adapted models to upcoming reinforced materials are needed. Models describing the conversion degree evolution of the composites as well as its direct link to mechanical response.

In this work a thermo-rheological-kinetical computation of the heating-compression of a reinforced sample is presented. In this study we observe not only the squeeze-flow mechanics characteristic of this process, but also the thermal evolution of the sample, including resin reaction. Both, thermal and kinetic evolution are coupled to the viscosity of the material. Thus, affecting the squeeze-flow nature.

MODELING FIBRE-REINFORCED COMPOSITES IN IMMERSED DOMAINS

Compression molding processes involve more complex physical phenomena than solving a fluid mechanics problem. In fact, modeling these processes come along with studying the coupling between the mechanical problem, the thermal

Fiber

Matrix

Mixture

Two phases Rheology

Homogeneised Rheology FIGURE 1. Fibre reinforced composites mainly described by two phases are described here as a homogeneous mixture characterized by its fiber fraction f . All the parameters and variables here described represents the mixture of both phases.

one and for thermosets, the polymerization or degree of curing reaction.

Governing Equations

Let us consider the homogeneous flow of a fiber reinforced composite fluid (Fig. 1) described by the Navier-Stokes equations, with an orthotropic non-Newtonian compressible constitutive law associated to appropriates boundary conditions. We extend the modeling here by adding the thermal evolution of the sample by describing not only its thermal state but also its thermo-kinetics due to resin composition by means of the term q react . We also included in this modeling the evolution of the conversion degree α, which accounts for the measurement of the reaction level. One writes the governing equations of the fully modeling as follows:

ρ ∂ v ∂t + v • ∇ v -∇ • σ = ρ g (1)
∇ • v = 0 (2) ρ C p ∂T ∂t + v • ∇T -∇ • (κ∇T ) -σ : ∇ v = qreact (3)
∂α ∂t + v • ∇α = F(α, T ) (4)
with v the velocity of the mixture, ρ the density of the mixture, σ the stress Cauchy stress tensor and g the force of gravity. The thermal parameters C p states for the heat capacity of the mixture, κ the conductivity, qreac the heat released by the reaction and F(α, T ) the evolution function of the conversion degree.

Immersed Approach

In this work, the level set method for tracking the composite interface Γ is used [START_REF] Coupez | AIP Conference Proceedings[END_REF]. Particularly, the linearized hyperbolic tangent ψ defined in (eq.5). Here, d(x, Γ) stands for the distance of a point to the interface, e is the thickness used to mix the physical properties along the domain and E a thickness where the level set function is truncated. Function ψ is convected to described the interface movement adding a reinitialization velocity v r conserving the metric properties of the function (eq.6).

ψ =                e + E tanh( d(x,Γ)-e E ) if ψ > e ±d(x, Γ) if -e < ψ < e -e + E tanh( -d(x,Γ)+e E ) if ψ < -e (5)
∂ψ ∂t + v + v r • ∇ψ = f (ψ) (6)
When solving (eq.6) volume conservation is not guarantee. For correcting this lost, a predictor-correction scheme is used. Conservation implies Ω H(φ)dΩ -V 0 = 0, where V 0 is the initial volume and H = 1/2 + ψ/2E, the smoothed Heaviside function associated to the level set. First we computed the volume lost each iteration by f n+1 = Ω H(ψ)dΩ-V 0 . Then, the conservative level set φ is obtained by (eq.7). Where δ(α) is a Dirac function associated to level set [START_REF] Zahedi | [END_REF].

φ n+1 = ψ - f n+1 Ω δ(α)dΩ . (7)
In this strategy, the correction is not space dependent. For more enriched methods the reader is referred to [10]. However, we take into account the presence of the mold and punch (tools) in order to not apply the correction causing penetration of the level set into the tools.

Thermo-rheological-kinetical coupling

In this section, we present the expression of the stress tensor σ used in (eq.1), that takes into account the anisotropy nature of the fiber network. Additionally, the expression of the heat rate q react used in (eq.3) is given. These energy released by the reaction represents the transformation of the chemical energy of the resin into heat. Both natures are related through the viscosity expression. During the compression forming, the pseudo-plastic viscosity decreases with the strain rate and the high temperature. However, the evolution of the kinetic trigger the polymer chains and entangled them causing the solidification of the piece, translating into a viscosity increment.

Rheology Model

The rheological model used in this work is an extension of the model found in [4]. The stress tensor is defined by (eq.8) where the anisotropy of the composite is described by the extra contributions associated to the deformation strain rate D described in (eq.9), by the tensor M and the parameters α 0 = 1.03, α 1 = 1.446, α 2 = -1.97. p states for the incompressible pressure.

σ = α 0 η eq D + α 1 M : D M + α 2 2 M • D + D • M -pI (8)
D = 1 2 ∇ v + ∇ t v (9)
D 2 eq = α 0 D : D + α 1 M : D 2 + α 2 (D • M : D) (10)
The model assumes plug flow, causing fibers to spread in an in-plane way. Averaging the fibers' directions modify homogeneously the macroscopic stress, leading to plane isotropy. Here, η eq states for the consistency of the paste under a concentration f of fibers, and states for its power-law exponent m and the equivalent strain rate D eq , here defined as shown in (eq.10). The orthogonal direction n represents the normal to the fibers' in-plane motion (here, the vertical direction) and its effect is represented in the model by means of the tensor M = n ⊗ n. Parameters α, account for fiber concentration f and expressions for these values are found in the original paper [4]. In addition, a numerical comparison of the response of the stress tensor for simple compression, plane strain and shearing cases without taking into account kinetics was performed in [START_REF] Salazar Betancourt | Numerical implementation of a rheology model for fiber-reinforced composite and viscous layer approach for friction study[END_REF].

Thermo-kinetic Model

The kinetic model is giving through the heat released qreact computed by means of the Bailleul's Model. Where the heat released is written by (eq.11):

qreac (T, α) = ρ (T, α) dα dt ∆H tot . ( 11
)
The computation of the evolution of the reaction rate F(α, T ) (eq.12) is usually described by the product of a temperature function K(T ), a function depending on the conversion degree G(α) and the induction period function h(T ) all described in (eq.13). The activation function W(h(t)) stands to zero if the ignition period has not been reached yet and for 1 once all the inhibitors of the resin has been consumed. These expressions read:

F(α, T ) = W(h(t)) K(T ) G(α), (12) K(T ) = k re f exp -A T re f T -1 ; G(α) = a i α i ; h(T, t) = t re f - t 0 exp -B T re f T -1 dt. (13)
In general, all thermal parameter depend not only on temperature T but also on the conversion degree α. From a numerical perspective, we focus on the modeling of the thermo-rhelogical-kinetic coupling of fiber reinforced materials. Expressions for the parameters k re f , A, T re f , a i , t re f , B depend on the resin composition of the composites. Expressions for the evolution of the thermal parameter C p, ρ, κ are not shown in here, due to confidentiality.

Coupling with viscosity

The thermo-kinetic-mechanical coupling is done by factors h(D eq ), f (T ) and g(α). The viscosity η eq is then calculated via η eq = η 0 h(D eq ) f (T ) g(α). Here f (T ) yields for the coupling with the thermal evolution, g(α) the kinetic coupling and h(D eq ) the pseudo plastic coupling. Expression for those contributions are presented in (eq.14, 15, 16) :

h(D eq ) = 1 + 98 f + 980 f 2 2 1 + H 1 + 2H -m+1 2 D eq m-1 (14)
f (T ) = e b( 1 T -1 T o ) (15) g(α) = 1 - α α gel s + g min (16)
The pseudo plastic coupling depends on the fiber fraction f = 14.7% and rheological parameters H = 0.4735, L = 46.87 and the sensibility on the strain rate m = 0.44. The thermal coupling follows an Arrhenius law type depending on activation energy b = 4500K and reference temperature T o = 296K. The kinetical coupling is giving by the evolution of the conversion degree translating into the solidification of the piece. Here, the coupling depends on α gel = 0.05 the reference conversion degree and s = 2.0 a hardening factor depending on the nature of the polymer composites. Here, g min = 0.1 states for a minimum kinetic factor when the conversion degree equals α gel .

COMPRESSION MOLDING OF FIBRE REINFORCED COMPOSITES

As example of fully coupled resolution, we present a test case where a sample is compressed at constant velocity (Fig. 2). The Punch and the Mold have a temperature of 150 • C, whereas the sample has an initial temperature of 40 • C. The sample is compressed at a constant speed of v = 0.004mms -1 while being heat up. The simulation depicts the evolution of the viscosity during the polymerization of the composite. As compression occurs, the hot plates heat up the sample until reaction is triggered. The temperature, shearing and the polymerization interacts modifying the viscosity of the sample. A mesh adaptation technique is implemented, bounded for a minimum mesh size of h min = 0.08mm. Mesh is adapted respect to velocity, temperature and level set fields. An adaptive time step is bounded to avoid a change in temperature per node higher than 5 • C. Also, the time step fulfill the CFL condition for the resolution of (eq. 6).

Evolution of the Conversion Degree

In order to analyze the reaction mechanism, we present here the evolution of the conversion degree α. As evidenced in Fig. 3 the reaction begins at time about 200s for the surfaces and at time equals 380s in the core. The evolution of the α value is well described by our discrete temporal scheme. The reaction in the core of the material occurs faster than in the surface, this due to the higher heat transfer in the center coming from the neighboring. Whereas the reaction on the surface described an extended hyperbolic tangent function, the reaction in the core describe basically a more narrow function. The reaction in the core last really few seconds compare to the reaction in the surface which is smooth. This mechanism is captured by our numerical scheme due to our adaptive time step.

Evolution of the temperature

In Fig. 4, the temperature evolution on the surface and on the center during the heating is depicted. We observe how in the core, temperature remains under the surface temperature until reaction takes place. The heat released is then concentrated in the core causing the temperature of the core to be the highest of the system. increasing even 50 • C more than the temperature of the hot plates. Once reaction is consumed, the thermal balance is achieved towards the temperature imposed by the plates.

Evolution of the viscosity

The thermal and kinetic evolution of the sample make the viscosity evolves. The increase in temperature tends to decrease the viscosity whereas the higher the conversion degree, the viscosity increases since solidification takes place. In Fig5 we compare the evolution of viscosity in the surface and the core. Basically this evolution show the combination of the shearing, thermal, and polymerization process. During about 250s the surface viscosity decreases for the thermal heating and just then raises up. The core of the sample raises up the viscosity much later than the surface. We conclude in this test case that conversion degree governs the viscosity of the sample.

η/η 0 [-]
Evolution of the viscosity through reaction and thermal evolution Core Surface FIGURE 5. Viscosity evolution during compression and heating of the sample. viscosity in the surface and in the core shows the interaction between the thermal softening and the kinetic hardening.

Impact of thermo-kinetics on mechanical response

In this section we present the modification of the mechanical response due to the thermo-kinetic coupling. In Fig. 6 the temperature profile along the sample shows the mapping of values during the reaction. This picture at t = 250s clearly shows the temperature profile during the compression. We see how the heat released by the reaction increases the temperature of the sample to higher values than the imposed temperature. For the time of the picture in Fig. 7, the solidified shell belong to reaction values of 1. Whereas, 0 values states for the non-reacted material. This shape gives us a indicator of how reaction grow in the sample. Begin in the surface, moves towards the center and towards the exterior.

In Fig. 8 the mapping of viscosity is depicted. Lower viscosity values are found outside of the solidification shell and close to the side extreme, where shearing is higher. The coupled of the thermal factor and the kinetic polymerization modify the viscosity. The velocity gradient are found to be naturally higher in the fluid zone. As example, the velocity in a intermediate point of the reaction is presented in Fig. 9. We observe how the velocity profile present strong gradient since part of the piece has solidified. Basically, the solidified shell works as punch and mold squeezing out the fluid towards the exterior. Normally, during compression molding, we would like to have first the filling of a part before the reaction takes place. However, we have seen in here that reaction in the surface takes place during almost the whole process. Of course, mold temperature plays a roll on the beginning of the reaction, but in this simple example, we can see how the reaction mechanism is depth depend, has different reacting times, has a different rate of conversion and modify the pattern of the flow. Details that should be taking into account in the conception phase of industrial piece designing. A more concluding results can be the evolution of the compression force. In Fig. 11 is shown the increment of the force once reaction occurs in the center. At time t = 400s the force has a 'jump'. By looking at the time, it can say that the macroscopic force is indeed modify when reaction in the core takes place. Basically it describes the length of the reaction in the core since it present the same shape as the conversion degree in the core described in Fig. 3.

CONCLUSIONS

This work has presented a coupled thermo-kinetic-mechanical model of fiber-reinforced composites during compression molding. In here, the rheology of the material is described by an one-equivalent-phase orthotropic model that accounts for the anisotropy exhibits for the fiber-network. For the thermal resolution a thermo-kinetic model is used to describe the polymerization mechanism. A thermo-rheo-mechanical coupling is then formulated by means of factors modifying the viscosity of the composites. A compression molding study case is studied showing the evolution of the viscosity, temperature, velocity and stress during reaction. We evidenced that reaction in the core of the material is much quicker than in the surface. That means that a mapping of viscosity values is presented during the reaction modifying the usual the squeeze flow. For our case, the highest temperature of the system was found in the core where reaction takes place as last but quicker than elsewhere. The reaction mechanism moves from the walls towards the center and extended then to the exterior. Thus, having different reacting times causing space-depend solidification zones. We can conclude as well that the length of the reaction in the core can be estimated by the macroscopic force evolution.

FIGURE 2 .

 2 FIGURE 2. Compression Molding of composite sample 40mm × 6mm. Mold and Punch tool at 150 C. Sample with initial temperature of 40 • C. Sample is compressed at constant speed v during reaction.

FIGURE 3 .

 3 FIGURE 3. Evolution of the conversion degree

FIGURE 4 .

 4 FIGURE 4. Temperature profile during temperature surface imposition. Temperature in the core shows the reaction kinetics.

FIGURE 6 .

 6 FIGURE 6. Compression molding of reinforced composites. Temperature field along the sample in Kelvins at t = 250s. We observe the hot zones where the reaction just took place.

FIGURE 7 .

 7 FIGURE 7. Conversion degree α. We observe the front of the reaction at t = 250s. The reaction begins in the surfaces and moves towards the center and towards the external boundaries.

FIGURE 8 .

 8 FIGURE 8. Viscosity field along the sample in Pascal. We observed zones with higher viscosity due to solidification at t = 250s.

FIGURE 9 .

 9 FIGURE 9. Velocity magnitude during reaction in mm s -1 at t = 250s. Solidification modifies the viscosity of the sample close to the surfaces. The common squeeze flow is not longer seem. during reaction, the material with lower viscosity tends to be squeezed.

FIGURE 10 .

 10 FIGURE 10. Stress magnitude at t = 250s in Pascal under reaction during compression.

FIGURE 11 .

 11 FIGURE 11. Compression Force evolution. When compressing the sample the squeeze force increased by the increasing of the surface in contact. The force slope is affected as well for the hardening due to polymerization.

ACKNOWLEDGMENTS

The authors would like to acknowledge to Romain Cardinaud for providing the kinetic-model used in this work and its advises during the numerical implementation.


G Roussel

X Ceamanos

X Briottet

A sun/shadow approach for the classification of hyperspectral data

Shadows constitute a major issue for various remote sensing applications, from unmixing and anomaly detection to image segmentation and classification. It is particularly problematic in urban areas, characterized by a high heterogeneity in terms of both relief, with the presence of numerous tall buildings, and materials, which implies that a single shadow can hide a wide variety of different objects. Shadow induces a strong weakening of the signal which can lead, in the context of a classification process, to the creation of a parasit shadow class for unsupervised cases, or to an increase of class variability for supervised ones. This paper introduces a new classification algorithm where sunlit and shadowed pixels are processed separately. Using a shadow mask, sunlit pixels are isolated and classified using a SVM or a Kmeans algorithm. Then, a centroid is processed for each classe and used by a spectral angle mapper (SAM) method to classify the shadowed pixels. This approach showed promissing results, especially in an unsupervised context.

Introduction

When shadows are omnipresent in an image, a de-shadowing preprocessing step is highly recommended in order to improve the classification results. Shadow processing is usually separated into two phases : shadow detection and shadow compensation. Adeline et al. [START_REF] Adeline | Shadow detection in very high spatial resolution aerial images: A comparative study[END_REF] proposed a very exhaustive taxonomy on shadow detection methods, which are classified in four categories including property-based, model-based, physics-based and machine learning methods. As for the shadow compensation phase, Li and al. [START_REF] Li | An adaptive nonlocal regularized shadow removal method for aerial remote sensing images[END_REF] have identified two kinds of methods, those that focus on the intensity domain [START_REF] Yamazaki | Characteristics of shadow and removal of its effects for remote sensing imagery[END_REF][START_REF] Lorenzi | A complete processing chain for shadow detection and reconstruction in vhr images[END_REF][START_REF] Arbel | Texture-preserving shadow removal in color images containing curved surfaces[END_REF] and those working in the gradient domain [START_REF] Finlayson | Entropy minimization for shadow removel[END_REF]. However, these methods either rely on a prior manual selection of paired regions, are imprecise when the shadowed areas are composed of a wide variety of materials or tend to alter significantly the non-shadowed areas. Moreover, they have been widely tested on natural multispectral images and rarely on aerial hyperspectral data, where the distance between the sensor and the target can be very important and imply distorsion and absorption effects from the atmosphere and reflection effects from the surrounding 3D structure. Yet, a third categorie of shadow compensation method, more suited to aerial images, can be considered. Indeed, a few atmospheric correction methods, which aim to process a reflectance image through an inverse radiative transfer equation, take shadows into account. The ICARE [START_REF] Lacherade | Icare: A physically-based model to correct atmospheric and geometric effects from high spatial and spectral remote sensing images over 3d urban areas[END_REF] algorithm can correct both hard shadows, casted by tall objects, and vague shadows, induced by slope orientations, using a precise 3D model of the observed scene. However this algorithm is very expensive in terms of computation time and requires a digital elevation model of the scene which is not always available. Then, Chen et al. [START_REF] Chen | Efficient empirical reflectance retrieval in urban environments[END_REF] introduced a faster and simpler method able to compute reflectance in shadow areas using only a radiance image and the atmospheric and viewing conditions, by approximating some of the radiative terms involved in the inversion.

Usually, shadow compensation and classification are two distinct processes. A shadowfree image is produced by the first step, and then processed by regular a classification algorithm, which can leads to misclassification when the shadowed areas aren't perfectly compensated. In this paper, we propose a new hybrid method where sunlit and shadowed pixels are classified sequentially in order to overcome the flaws of the shadow compensation step. First, an hyperspectral adaptation of the Chen et al. [START_REF] Chen | Efficient empirical reflectance retrieval in urban environments[END_REF] empirical method is introduced. Then our hybrid classification approache is detailed, followed by a discussion on the classification results.

Empirical reflectance processing method

This preprocessing step aim to compute for each pixel a reflectance value, unitless and wavelength-dependant, representing the effectiveness of a material in reflecting radiant energy. The method we chose in this study is an empirical algorithm introduced by Chen et al. [START_REF] Chen | Efficient empirical reflectance retrieval in urban environments[END_REF]. The main asset of this algorithme is its ability to take shadows into account, which implies that it can be considered as a proper shadow compensation method. Assuming a lambertian hypothesis and for given illumination conditions, the reflectance can be expressed as :

ρ = π × (R tot -R env -R atm ) (I dir + I dif + I coup + I ref l ) × τ ↑ dir (1)
where R tot is the radiance value received by the sensor. R env is the portion of R tot due to the neighborhood of the target pixel and R atm the one scattered by the atmosphere to the tagert withour interacting with the ground. I dir , I dif , I coup and I ref l are portions of the overall irradiance striking the target pixel. I dir stands for the photons coming directly from the top of atmosphere, I dif for those which are scattered at least once by the atmosphere, I coup for those which make at least one round trip between the ground and the atmosphere and I ref l for those which are reflected by the environment before hitting the ground. Finally, τ ↑ dir is the direct upward transmission. All these terms are pixel dependant except for I dir , I dif , R atm and τ ↑ dir which are constant for the whole scene. In the original algorithm, these constant terms were computed by the radiative transfert processing software 6S [START_REF] Vermote | Second simulation of the satellite signal in the solar spectrum, 6s: an overview[END_REF] but we chose to replace it by MODTRAN [START_REF] Berk | Modtran : A moderate resolution model for lowtran 7[END_REF], more suited to hyperspectral data due to its better spectral resolution.

As a first step of the algorithm, a shadow mask is processed using an automatic shadow detection method introduced by Nagao et al. [START_REF] Nagao | Region extraction and shape analysis in aerial photographs[END_REF], which was pointed fot its efficiency by Adeline et al. [START_REF] Adeline | Shadow detection in very high spatial resolution aerial images: A comparative study[END_REF]. The mask is obtained by applying an automatic histogram thresholding to the following combination of bands :

R comb = 2R λr + R λg + R λ b + 2R λ nir 6 ( 2
)
where R λr , R λg , R λ b and R λ nir are respectively the radiance values associated to the red, green, blue and near infrared channels of an image. In order to discriminate water from shadows, a final step has been added where the detected shadow pixels are thresholded a second time, only on the green channel which showed the best ability for this discrimination. Shadowed and sunlit pixels are then processed separately. In sunlit areas, I dir and I dif constitute most of the irradiance received by a pixel. We consider then that the coupling and reflective terms can be neglected. In shadow areas, there is no direct irradiance and the coupling term is also put aside. Therefore, the total irradiance is reduced to a diffuse and a reflective term. However in this case, because of the proximity of buildings, the diffuse irradiance is weighted by a skyview factor α sky . Using a DEM of the downtown of Toulouse, France, we observed that the shadowed areas of this kind of environment are mainly composed of narrow streets and half-opened areas, and that skyview factors belong to [0.45, 0.65]. A mean value of 0.55 is then chosen. As for the reflected irradiance, it is processed as the sum of direct and diffuse irradiances weighted by an attenuation factor κ and a hypothetical mean reflectance of the scene ρ :

I ref l = κ × (I dir + I dif ) × ρ (3)
Chen et al. [START_REF] Chen | Efficient empirical reflectance retrieval in urban environments[END_REF] recommend to set κ at 0.2 for cluttered urban environment and 0 for suburban areas. As our data focus on the downtown of Toulouse, only the first value will be used. For ρ, a combination of several urban materials (concrete, asphalt, tile and vegetation) gathered from spectral libraries is used. At this point, only one radiative term remains, the environment radiance R env . Contrary to S, it is not a direct outpout of MODTRAN, si it needs to be computed iteratively as follows :

R env,t = τ ↑ dif f • I tot • ρt 1 -ρt • S (4)
ρ t = π(R tot -R env,t -R atm ) I tot × τ ↑ dir (5)
where ρt+1 = (u,v)∈V (x,y)

ρ t (u, v) • du • dv (6)
is the mean reflectance spectrum associated to the neighborhood of the target for the iteration t + 1 and ρ0 = ρ. τ ↑ dif f , the diffuse upward transmission and S, the spherical albedo, are processed by MODTRAN. Usually, two iterations of the process suffice to reach the convergence. The graphic on figure 1 shows four spectra corresponding to two vegetation pixels (one in the shadow and the other in the sun light) and two atmospheric correction methods, the studied empirical reflectance processing method and COCHISE [START_REF] Miesch | Direct and inverse radiative transfer solutions for visible and near-infrared hyperspectral imagery[END_REF], an approach which doesn't take shadows into account. For the empirical method, the reflectance spectra of the shadowed spectrum is slightly weaker than that of the sunlit spectrum, but the shape of a vegetation spectrum is still easily recognizable. This is not the case for the COCHISE approach, where the shadowed spectrum is strongly flattened. The benefits of the empiral method can also be noted by comparing a luminance image and a reflectance image on the figure 2. However, we can see that the shadow areas are not perfectly corrected, which is why they have to be taken into account in the classification process.

Sun/shadow classification method

In order to complement the empirical atmospheric correction method, we propose a solution implying a sequential process (cf. figure 3) where sunlit and shadow pixels are classified separately. Using the shadow mask described in section 2, the sunlit pixels are isolated and classified through a Kmeans or a SVM algorithm, with a number of classes C fixed by the Then, the cendroids of each class are preserved and used as input for the classification of shadow pixels. In this case, a spectral angle mapper (SAM) classification has been chosen for its ability to focus on the shape of spectra, which remain globally similar for a pair of pixels associated to the same material, even if one is located in a sunlit area and the other in a shadow area. For each pixel p and each centroid µ i , i ∈ [1, C], a SAM value is processed as follow :

SAM (p, µ) = cos -1       B b=1 p(b)µ(b) B b=1 p(b) 2 B b=1 µ(b) 2       ( 7
)
where B is the number of bands of the image. The class chosen for pixel p is the one minimizing the spectral angle mapper between p and its representative centroid. It should be noted that for this method to be efficient, the same classes must be present in both sunlit and shadow areas.

Results

The classification algorithm is applied on a section of a dataset acquired during the Umbra campaign over the downtown of Toulouse (cf. figure 4a). Two Hyspex hyperspectral sensors have been used in order to cover a wavelengths interval starting from 0.4µm to 2.5µm, with a spatial resolution of 1.7m. The noisy bands located at the extremities of the spectrum or in the middle of water vapor absorption bands (centered on 0.94, 1.1, 1.38 and 1.87µm) have been removed, resulting in a set of 255 bands. Four classification approach have been applied to the image : a support vector machine algorithm (SVM), a K-means (KM), and our new method with a SVM (SS-SVM) and a K-means (SS-KM) for the sunlit pixels classification. This latter is applied on the whole set of 255 bands but to process shadowed pixels, we used only the hundred first bands since after 0.9µm, shadowed areas are too noisy. Five classes have been considered : Vegetation (green), water (blue), tile (orange), asphalt (black) and gravel (white). The results on figure 4 show that the SVM and KM methods don't work well in shadow areas, where most pixels are classified as asphalt. The two sun/shadow classifications on the other hand manage to discriminate water and vegetation pixels even in the shadow. The only issue remains in gravel and asphalt classes which still present strong confusions. The improvement induced by the sun/shadow approach is also observable in the increase of overall accuracies from simple SVM and K-means to sun/shadow SVM and K-means (cf. table 1).

Conclusion

A new sun/shadow approach have been introduced for the classification of hyperspectral data. This method involve several preprocessing steps, including an atmospheric correction algorithm and the creation of a shadow mask, before classifying sunlit and shadowed pixels In either cases, it appears that taking shadows into account in the classification process can strongly improve its accuracy. Therefore, this approach is a good complement to shadow compensation techniques, especially in a hyperspectral context where these latter are still rare and perfectible.

Figure 1 :Figure 2 :

 12 Figure 1: Comparison between a shadow and a sunlit reflectance spectrum obtained with two correction algorithms : COCHISE and the empirical method

Figure 3 :

 3 Figure 3: Sun/shadow classification

  Sun/shadow SVM (e) Kmeans (f) Sun/shadow Kmeans

Figure 4 :

 4 Figure 4: Classification results obtained on four points of interest in the image

Table 1 :

 1 Overall accuracies and Kappa for the four classification methods user.

		SVM SSVM	KM	SSKM
	OA	89.2% 95.1% 80.9% 94.7%
	Kappa	0.86	0.93	0.75	0.93

Acknowledgment

The authors would like to acknowledge the ANR HYEP for the financial support of this work, the Onera institute for the provision of the hyperspectral data and the IGN institute for the coregistration of the two HYSPEX sensors.


R Clerc

B Bouthinon

M Mohankumar

P Rannou

J Vaillant

T Maindron

B Racine

Y-F Chen

L Hirsch

J M Verilhac

A Pereira

A Revaux

Current status and challenges of the modeling of organic photodiodes and solar cells Invited paper

Progress in the modeling of charge transport in solution processed solar cells and photodiodes is reviewed.

Through several examples involving modeling and original experiments, the role of intentional doping, structural defects, and oxygen contamination are discussed.

I. INTRODUCTION

Organic electronics is an emerging technology using conjugated molecules or polymers (i.e. organic semiconductors) [START_REF] Facchetti | Π-conjugated polymers for organic electronics and photovoltaic cell applications[END_REF] to produce electronic devices. Batteries, photovoltaic, (bio) sensors, displays, lighting are among the most relevant fields of applications of this technology [START_REF] Lupo | OE-A roadmap for organic and printed electronics[END_REF][START_REF] Deen | Flexible electronicsopportunities and challenges[END_REF][START_REF] Jacob | Large area sensing surfaces: Flexible organic printed interfacing circuits and sensors[END_REF]. One particular challenge of this domain consists in taking advantage of solution-processed organic materials to produce low-cost and large area devices by either roll-to-roll, sheet-tosheet or printing technologies [START_REF] Deen | Flexible electronicsopportunities and challenges[END_REF]. In this context, solution processed photodiodes are of particular interest, either for outdoor [START_REF] Kang | Bulk-Heterojunction Organic Solar Cells: Five Core Technologies for Their Commercialization[END_REF][START_REF] Zhao | Efficient organic solar cells processed from hydrocarbon solvents[END_REF] and indoor [START_REF] Lechêne | Organic solar cells and fully printed supercapacitors optimized for indoor light energy harvesting[END_REF] photovoltaic or imaging applications [START_REF] Verilhac | Recent developments of solution-processed organic photodetectors[END_REF][START_REF] Tedde | Imaging with organic and hybrid photodetectors[END_REF][START_REF] Pace | Printed photodetectors[END_REF]. In the field of solar energy, External Power Conversion Energy (EPCE) exceeding 10 % have been reported using Bulk HeteroJunction (BHJ) polymer based devices [START_REF] Zhao | Efficient organic solar cells processed from hydrocarbon solvents[END_REF] (approaching its fundamental limits [START_REF] Janssen | Factors limiting device efficiency in organic photovoltaics[END_REF]) and more than 20 % using solution processed hybrid inorganic/organic perovskite semiconductors [START_REF] Yi | Perovskite Photovoltaics With Outstanding Performance Produced by Chemical Conversion of Bilayer Mesostructured Lead Halide/TiO2 Films[END_REF]. Organic imagers are also attractive, either to improve silicon CMOS imagers [START_REF] Mori | Thin organic photoconductive film image sensors with extremely high saturation of 8500 electrons/µm 2[END_REF], or to realize large-area, flexible, low-weight full organic imagers. Dark current level comparable with silicon photodiodes have already been reported [START_REF] Pierre | High Detectivity All-Printed Organic Photodiodes[END_REF] [START_REF] Saracco | Work function tuning for high-performance solution-processed organic photodetectors with inverted structure[END_REF]. In particular, large-area organic imagers, either deposited on the top of a thin-film transistor (TFT) or organic transistor backplane are attractive alternatives to amorphous silicon imagers in X-ray detectors for medical applications [START_REF] Tedde | Imaging with organic and hybrid photodetectors[END_REF] [START_REF] Blakesley | Modeling the imaging performance of prototype organic x-ray imagers[END_REF]. In all cases, despite the instability of most organic material in presence of UV, water and/or oxygen, the reliability of polymer-based devices has greatly improved and is no longer considered as a major obstacle, provided that suitable materials, processes, and encapsulations are used [START_REF] Peters | High efficiency polymer solar cells with long operating lifetimes[END_REF] [START_REF] Cheng | Stability of organic solar cells: challenges and strategies[END_REF].

With the emergence of this technology, device modeling is more and more needed, either to support device processing or to predict system performances. The aim of this paper is to review the existing physical models, to discuss the main ingredients needed for device simulations, and to give example of application of modeling, in close comparison with experiments.

II. CHARGE TRANSPORT IN DISORDERED ORGANIC

SEMICONDUCTORS.

It is generally admitted that -conjugated small molecules molecules and polymers behave as disordered semiconductors [START_REF] Bässler | Charge transport in organic semiconductors[END_REF], even if the exact origins of such disorder are still under debate, and are material/process dependent. Structural heterogeneities in amorphous or polycrystalline semiconductors, induced by processes and interfaces, are of course a major sources of disorder, and many studies have reported significant mobility improvements by optimization of the material morphology [START_REF] Tsutsui | Unraveling Unprecedented Charge Carrier Mobility through Structure Property Relationship of Four Isomers of C12-BTBT-C12[END_REF]. However, it seems that disorder may have also intrinsic origins, induced either by chemical impurities, by changes in the surrounding polarization during charge transport [START_REF] Bässler | Charge transport in organic semiconductors[END_REF], or as recently proposed, induced by large thermal molecular motions (transient localization), a direct consequence of the weak strength of van der Waals intermolecular interactions [START_REF] Fratini | The transient localization scenario for charge transport in crystalline organic materials[END_REF]. This later theory is supposed to explain the low temperature mobility improvement reported in some high mobility (> 1 cm 2 .V -1 .s -1 ) organic single crystals, and sometimes referred as "band like" transport [START_REF] Sakanoue | Band-like temperature dependence of mobility in a solution-processed organic semiconductor[END_REF].

In the literature, two main approaches have been used to model transport in disordered semiconductor. The first one consists in solving the hopping master equation, either by Monte Carlo [START_REF] Bässler | Charge transport in disordered organic photoconductors a Monte Carlo simulation study[END_REF] or by direct solving [START_REF] Jungemann | DC, AC and noise simulation of organic semiconductor devices based on the master equation[END_REF]. Quantum chemistry can be used to calculate hopping site energy distributions [START_REF] Kordt | Modeling of organic light emitting diodes: From molecular to device properties[END_REF]. In many cases however, the disordered energy positions of hopping sites is assumed to follow a Gaussian law, and hopping probabilities are modeled by Miller-Abrahams expressions, leading to the so-called Gaussian Disorder Model (GDM) [START_REF] Bässler | Charge transport in organic semiconductors[END_REF]. This model predicts an increasing mobility with temperature, electric field and carrier concentrations, in agreement with several experiments [START_REF] Pasveer | Unified description of charge-carrier mobilities in disordered semiconducting polymers[END_REF]. There are exceptions, such as P3HT [START_REF] Mozer | Negative electric field dependence of charge carrier drift mobility in conjugated semiconducting polymers[END_REF] or the high mobility "band like" single crystal for instance [START_REF] Sakanoue | Band-like temperature dependence of mobility in a solution-processed organic semiconductor[END_REF]. Numerically expensive, this approach is typically used to get insight on transport [START_REF] Pasveer | Unified description of charge-carrier mobilities in disordered semiconducting polymers[END_REF], recombinations [START_REF] Groves | Bimolecular recombination in polymer electronic devices[END_REF] noise properties [START_REF] Jungemann | DC, AC and noise simulation of organic semiconductor devices based on the master equation[END_REF] or blend morphology [START_REF] Meng | An improved dynamic Monte Carlo model coupled with Poisson equation to simulate the performance of organic photovoltaic devices[END_REF] of organic materials (rather than to simulate complete device).

An illustration of the capability of GDM models to reproduce a large set of experimental data is shown in Fig. 12345. In these experiments, several asymmetrical -NPB p-only devices have been realized by evaporation, featuring different thicknesses in the 20-110 nm range. -NPB is a hole transport material (HTM), classically used in OLED [START_REF] Brütting | Device physics of organic light-emitting diodes based on molecular materials[END_REF]. Interestingly, the conductivity of -NPB has been tuned by molecular doping of MoO 3 , as in [START_REF] Kröger | P-type doping of organic wide band gap materials by transition metal oxides[END_REF]. In agreement with GDM theory, the conductivity shows a clear exponential dependency with doping concentration (Fig. 2), despite significant process induced variability. Moreover, using low-field conductivity extracted on each device, I-V experiments can be satisfactory reproduced with the same set of GDM parameters, provided that the mobility electric-field dependence is included [START_REF] Brütting | Device physics of organic light-emitting diodes based on molecular materials[END_REF].

III. DRIFT DIFFUSION (DD) MODELS FOR ORGANIC SEMICONDUCTOR DEVICES.

The second commonly used approach to model organic device is the well-known Drift Diffusion model, numerically efficient and versatile. Disorder is included via Gaussian Bands (or exponential bandtails), shallow and deep traps and GDM inspired mobility models, as it was done in the past for amorphous silicon solar cells [START_REF] Hack | Theoretical modeling of amorphous silicon-based alloy p-i-n solar cells[END_REF]. It can account for electron and holes transport in heterostructures [START_REF] Lundstrom | Numerical analysis of heterostructure semiconductor devices[END_REF] and interface charges and dipoles [START_REF] Brütting | Device physics of organic light-emitting diodes based on molecular materials[END_REF], which both play an important role in organic devices. Light propagation is performed by solving Maxwell equations [START_REF] Sievers | Modeling optical effects and thickness dependent current in polymer BHJ solar cells[END_REF]. DD requires detailed calibration with experiments, and especially band structure parameters (using UPS), mobility (using Time of Flight TOF), optical index and recombination (deep and shallow traps, bandtails). Even if the applicability of DD to model organic semiconductors is questionable [START_REF] Stallinga | Electronic transport in organic materials: comparison of band theory with percolation/hopping theory[END_REF], the efficiency of this approach to tackle many sophisticated phenomena such as transient low-light intensity in organic photodiodes [START_REF] Popescu | Modeling and simulation of organic photodetectors for low light intensity applications[END_REF], impact of blend morphology in BHJ solar cell [START_REF] Ray | Random vs regularized OPV: Limits of performance gain of organic bulk heterojunction solar cells by morphology engineering[END_REF] [START_REF] Bouthinon | Impact of blend morphology on interface state recombination in bulk heterojunction organic solar cells[END_REF], role of electron/hole transport layer (sometimes causing S shape [START_REF] Tress | Influence of Hole-Transport Layers and Donor Materials on Open-Circuit Voltage and Shape of I-V Curves of Organic Solar Cells[END_REF]), analysis of impedance spectroscopy experiments [START_REF] Bisquert | Theory of impedance and capacitance spectroscopy of solar cells with dielectric relaxation, drift-diffusion transport, and recombination[END_REF] (and many others) is remarkable.

In that respect, let us give a recent example illustrating the strengths and weaknesses of the drift diffusion model (Fig. 678). The p-only polycrystalline perovskite MAPbI 3 has been processed and characterized at low temperature. Accounting for volume and interface shallow traps, I-V curves at 10 K have been nicely fitted by a simple DD model (Fig. 7), extracting data from electroluminescence, I-V, TOF and C-V experiments (not shown here). Interestingly, with no additional fitting parameters, the model also correctly predicts the higher temperature I-V curves, however with a lower accuracy (Fig. 8). A closer look to experiments has indeed revealed that the experimental temperature activated conductivity was in better agreement with hopping models, suggesting the occurrence of transport from traps to traps, not included in the DD model.

IV. AN EXAMPLE OF DEVICE SIMULATION: FRONT AND BACK SIDE QUANTUM EFFICIENCY DIFFERENCES IN SEMI-TRANSPARENT ORGANIC PHOTODIODES CAUSED BY OXYGEN CONTAMINATION

To conclude, let us give another example of the capability of device simulation to support device processing. Semitransparent organic photodiodes may feature significantly different external quantum efficiency (EQE) when measured from front and back sides (see Fig. 910111213). It was initially suspected that this difference was caused by the different optical properties of top and bottom layers. However simulations has shown that optics can hardly explain of fraction of this discrepancy, and not its thickness and voltage dependency. On the other hand, the introduction of deep acceptor traps, inducing a negative charge within the active layer in presence of light generated carriers was able to explain this effect (Fig. 11&12). In fact, the introduction of a single level acceptor deep traps, inducing both recombination and space charge effects, was able to reproduce a large set of experimental data (Fig. 12). As oxygen is known to potentially induce negative charges and trapping [START_REF] Seemann | Reversible and irreversible degradation of organic solar cell performance by oxygen[END_REF], new experiments have been performed in absence of oxygen. In this later case, in agreement with simulations (Fig. 13), the difference of EQE has almost disappeared. The comparison between experiments and DD simulations has thus been of great help, indicating oxygen contamination as the potential origin of this effect.

V. CONCLUSIONS

Despite the complexity of disordered organic materials, Drift Diffusion based device simulation remains a powerful tool to analyze organic solar cells and photodiodes, and in particular the role of intentional doping, structural or induced defects states and interfaces. It requires however experimented users, capable on one hand to improve models thanks to the progress of more sophisticated approaches (such as variable range hopping model and quantum chemistry) and on the other hand, to perform detailed calibrations with experiments.

Fig. 1 :

 1 Fig. 1 : Band Diagram of the p-only asymmetrical -NPB device Fig. 2 : Low-field conductivity extracted from I-V experiments on each device versus nominal molecular weight (error bar indicates variability)

Fig. 6 : 3 Fig. 7 :

 637 Fig. 6 : Schematic band diagram of the p-only polycrystalline perovskite MAPbI 3

Fig. 8 :

 8 Fig. 8 : I-V curves for several temperature of the device of Fig. 6. Simulations are in good agreement with experiments, but do not perfectly capture the thermally activated low field regime.

Fig. 9 :

 9 Fig. 9 : Composition of the semitransparent organic photodiode, lighted by either the back side (substrate) or the front side.

Fig. 10 :Fig. 11 :Fig. 12 :Fig. 13 :

 10111213 Fig. 10 : Experimental ratio between the back side external quantum efficiency EQE and the front side @ -2V versus wavelength

ACKNOWLEDGMENT

This work benefited from the support of the project TAPIR ANR15CE24002405 of the French National Research Agency.


Isabelle Dahman

Nicolas Jeannin Demr

Philippe Arbogast

Bouchra Benammar

Attenuation Forecasts Model Exploiting Short Range Probabilistic Weather Forecasts

Keywords:

This paper proposes a method to control the data rate of an Earth Space link using VCM at Ka band based on the use of meteorological forecasts. Specifically, a probabilistic model of propagation impairments forecast has been developed from the analysis of rain rate measurements and outputs of the Météo-France short-range probabilistic forecast model. The approach is validated against propagation data acquired on a GEO to ground link and applied to the case of a LEO data down link using VCM, where it is shown to outperform significantly the data rate achievable using a purely statistical margin, while keeping the same availability target.

I. INTRODUCTION

High resolution satellite imagery programs require ever increasing data rates for data download links which are no longer compatible with the limited bandwidth available at X band (8-10 GHz), the frequency band usually used to retransmit data on Earth. In order to increase the capacity of such links, the use of higher frequency bands, for which larger bandwidths are available, becomes a necessity for future Earth observation satellites. It is therefore envisaged to operate communications between satellites and ground stations at Ka band (around 26.5 GHz) and in later future at Q band (around 40 GHz). A counterpart to the large bandwidth available at those frequency bands is a higher sensitivity of radiowaves to rain, with rain attenuations reaching potentially several tens of dB some tenth of percent of the time. In addition, considering the fluctuations of the elevation angle during the visibility period, the magnitude of the propagation impairments is strongly varying. To avoid the severe constrains set by the use of a single modulation and coding for the link in terms of power margins required disregarding the elevation angle, different techniques can be used. Firstly, the power margin to fulfill a given availability objective can be adapted by conditioning it to the elevation angle. This technique is called variable coding and modulation (VCM [START_REF] Cossu | Effects of link availability on the achievable performance with variable coding modulation earth observation satellites[END_REF]) and requires an adaptation of the physical layer depending on the elevation of the data download link. No feedback on the state of the link is required and the on-board modulator is programmed based on orbital data. This mechanism still require potentially high margin to cope with rather infrequent propagation conditions. To increase the efficiency of the link, the adaptation in quasi real time of the waveform to the state of the channel through adaptive coding and modulation is extremely efficient [START_REF] Holm | Adaptive Code Modulation Performance and Channel Estimation Tools for Flat Fading Channels[END_REF], but requires a quasi-instantaneous feedback of the state of the channel to be sent to the modulator through a return channel. The existence of this return channel raises various issues in terms of payload complexity, security of the telecommand (TC), and required infrastructure. It will need significant development before entering in operation. An intermediate solution would be to adapt the waveform for the next contacts with the ground station network each time the satellite is in visibility of the network control center through standard TC link. This requires a probabilistic forecast of the state of the propagation channel for the following contact periods with the receiving ground station network. This forecast can be sent at each contact period by the network control center to program the on-board modulator for the next visibility period for which data have to be downloaded on ground stations.

One of the difficulty related to the practical implementation of such a strategy that is referred to as dynamic VCM in [START_REF] Toptsidis | Link adaptation for Ka band low earth orbit earh observation systems : a realistic performance assessment[END_REF] (D-VCM) is to develop a probabilistic forecast model of the state of the channel as this information is not directly related to the information contained for instance in numerical weather forecasts. In addition, the accuracy of the weather forecast has to be taken into account. One of the issues is also to get an evaluation of the order of magnitude of the achievable improvements in terms of performance with regards to conventional VCM approach.

Although the presence of cloud, water vapor and oxygen also alters the radiowave during its propagation into atmosphere; this study is focused as a starting point on the forecast of rain attenuation that causes the higher impairments.

In this work, a model to forecast rain attenuation on an Earth space channel for a time horizon of a few hours or days is presented. It makes use of a probabilistic weather forecast data from a model developed by Meteo-France named PEARP, and is parameterized using propagation data collected in the south of France by ONERA and using the support of rainfall estimates given by rainfall radar measurements.

In a first stage of this paper, details on the proposed model are given. Then it is validated against fixed Earth space propagation data with a GEO satellite as no propagation data are available for the LEO to ground channel. In a later stage, the performances of D-VCM technique are assessed on a typical Ka band data download link of LEO satellite.

The main originality of this work lies on the use of concurrent and independent data sources for the propagation and weather forecast data and it can thus be representative of the performance of a system actually using those models.

II. ATTENUATION FORECAST MODEL

A. Approach

The proposed model aims at computing the probability to exceed a given rain attenuation threshold, knowing rain meteorological forecasts. This information is then used to deduce the power margin required to match a given availability level.

The learning process of the model requires rain forecast, rain rate measurements and rain attenuation measurements.

Meteorological Forecasts

The accuracy of the attenuation forecasted by the model mainly depends on the quality of the meteorological forecasts. Deterministic models are the most common type of meteorological models: analysis of initial conditions result in a single forecast such as the forecast of a rain fall rate. There is no assessment of the uncertainty of the prediction. However, knowing the trustworthiness of data used to predict the impairments of the propagation channel could greatly help to limit false alarms and non-detection cases.

Preference is given to probabilistic approach: the probabilistic models PEARP compute 35 realizations of the forecast in order to characterize both of uncertainty on initial conditions and error in the modelization. The closest the forecasts of the ensemble are, the more confidence one might have in the forecast.

PEARP [START_REF] Descamps | The Météo-France short-range ensemble prediction system[END_REF] benefiting of a global coverage with a horizontal resolution of 10 km over France (60 km on the opposite side of the globe) and forecast up to 36 h with a step of 3 hours, it appears to be suitable for studies presenting variables requirements such as different spatial locations.

In the whole study, a good compromise between the granularity of the forecasts and the possibility to offer a statistic link between rain rate and rain attenuation has been find selecting an area of 100 km around the receiving station (latitude 43.5 °N and longitude 1.5 °E) and averaging PEARP data on it.

Rain Measurements

Rain rate measurements are also provided by MeteoFrance in the form of radar quantitative precipitation estimates, called PANTHERE [START_REF] Parent-Du Chatelet | The PANTHERE project and the evolution of the French operational radar network and products: Rain-estimation, Doppler winds, and dual-polarization[END_REF]. PANTHERE product has a spatial resolution of 1 km² over France and a temporal one of 5 minutes. An illustration is given on figure 1.

Rain rate distribution conditioned to rain forecasts is

Rain Attenuation Measurements

Ka band attenuation measurements have been obtained by ONERA analyzing power fluctuation of the Ka band beacon of Astra 3B geostationary satellite. The link between Astra 3B and the receiving site has an elevation angle of 35°. Although attenuation data are available with a temporal resolution of 1 second for the full year 2014, data have been averaged on 5 minutes to save computational resources.

Errors in weather forecasts are firstly characterized analyzing the correspondence between rainfall radar measurements PANTHERE and the rain forecasts PEARP. In a second stage, a relationship between PANTHERE data and the Ka band rain attenuation measurements is established.

The expression of the probability of exceeding a given rain attenuation threshold A* knowing the rain meteorological forecast F is then given by (1).

𝑃(𝐴 > 𝐴

* |𝐹) = ∑ 𝑃(𝐴 > 𝐴 * |𝑐 = 𝑐 𝑗 )𝑃(𝑐 = 𝑐 𝑗 |𝐹) 𝑗 (1)
With A, the attenuation due to rain, c the rainfall measurements Panthere and F rain forecasts PEARP. The probabilistic forecast model of the state on the channel is therefore based on the computation of two conditional probability distributions:

-The rain attenuation distribution conditioned to rainfall radar measurements.

-The rain fall measurements distribution conditioned to rain forecasts.

B. Rain rate distribution conditioned to rain forecasts

To allow a comparison between PANTHERE data and PEARP outputs, PANTHERE data have been averaged on space time domains covering 100 km² and 1 hour. The cumulative distribution functions of rain fall measurements conditioned to PEARP rain forecast is presented in figure 2.

The PANTHERE rainfall probability distributions conditioned to PEARP demonstrates unsurprisingly that the higher the forecasted rain rate, the greater the probability to have a high rain rate.

C. Rain attenuation distribution conditioned to rain measurements

The estimation of rain attenuation based on PEARP forecasts is obtained coupling the rain fall measurements distribution presented above with the rain attenuation distribution conditioned to PANTHERE measurements.

The ccdf of Ka band rain attenuation conditioned to PANTHERE derived rain amount presented in figure 3.

Rain attenuation distribution conditioned to rain forecasts D. Rain attenuation distribution conditioned to rain forecasts

As described by equation ( 1), the above conditional probability distributions leads to the computation of the probability to exceed a given rain attenuation threshold, knowing rain meteorological forecasts. The result is illustrated on figure 4. For each realization of the forecast computed by the PEARP model, an estimation of the attenuation not exceeded p% of the time can be deduced from the previous graph. The mean of the attenuation forecasts thus obtained gives us a good estimation of the power margin that shall be provided to fulfill availability equal to 1-p % of the visibility duration. This information will allow optimizing the transmission adapting the waveform to the state of the channel.

A method to increase the efficiency of the link selecting the most suitable modulation and coding based on the attenuation forecasted is presented in the next section.

III. SYSTEM SCENARIO

Communications with LEO satellites are constrained by the visibility periods of the receiving stations. The deterioration of propagation conditions during this time window, especially for systems using Ka bands or above, may result in transmission outage and consequently in loss of data. The requirements in terms of minimum availability can be satisfied by applying to the link a couple modulation/coding that close the link budget for the given propagation margin that corresponds to the requested availability.

The ModCods considered in this study are presented in Table 1.

The idea is to select the ModCod leading to the higher efficiency 𝜂 while keeping an error free transmission. To achieve this, the carrier to noise ratio of the link must exceed a preset threshold C/N required .

The expression of the Carrier to Noise Ratio of the link not exceeded for p% of the visibility time and an elevation angle θ is given by equation [START_REF] Holm | Adaptive Code Modulation Performance and Channel Estimation Tools for Flat Fading Channels[END_REF].

The most effective strategy is to adapt in real time the transmission scheme to propagation conditions. A detailed description and an evaluation of the limits of this solution, called ACM (for Adaptive Coding and Modulation solution), are available in [START_REF] Holm | Adaptive Code Modulation Performance and Channel Estimation Tools for Flat Fading Channels[END_REF] and [START_REF] Jeannin | Sizing and optimization of high throughput radiofrequency data down link of earth observation satellites[END_REF].

As illustrated in the figure 5, the choice of the ModCod is made at the level of the receiving station on the actual instantaneous attenuation A converted into effective carrier to noise ratio. The satellite is then informed of the most efficient ModCod usable under current propagation conditions. The communication between the satellite and the Earth stations is optimized and much higher mean capacities are offered by the link. From DVB-S2X [START_REF]Second generation framing structure, channel coding and modulation systems for broascasting, interactive services, news gathering and other broadband satellite application Part II : s2-extensions (DVB-S2X)[END_REF].

ModCods Modulation

With EIRP the equivalent isotropic radiated power, expressed in dB. B is the Bandwidth of the link. G/T represents the figure of merit in clear sky conditions of the receiving station. kB is the Boltzmann constant. FSL represents the Free Space Loss which are function of θ. A is the power margin required to ensure an availability p% of the time for an elevation 𝜃.

However, a return channel, from the receiving station to the satellite, is required to transmit the state of the channel to the on board modulator. This ACM is not yet operational and requires further developments. Therefore, in absence of a return link, a pre-programming of the ModCod for the next orbit is necessary.

Free space losses significantly grow with decreasing elevation angle. The simplest strategy to guarantee the required availability is to assume a worst case link budget corresponding to the lowest elevation suitable for transmission. Assuming that a minimum elevation angle of 5° is required to operate the transmission, the power margin A(p) is equal to the attenuation not exceeded p% of the time at 5° of elevation. The ModCod selected is thus the most efficient that fulfill the required ( A first degree of adaptation consists in associating to each elevation angle the most efficient ModCod that close the link budget.

A schematic representation is given in figure 6. At the level of the control station, a plan of ModCods that statistically fulfill, for each elevation angle θ, the required availability p(%) is established and communicated to the satellite. During the visibility period of the receiving station, the satellite applies to the transmission the pre-programmed ModCods. Compared to CCM, this technique called Variable Coding and Modulation (VCM [1][6]) allows operating the link more efficiently but still do not allow any flexibility with regards to weather conditions and limits the throughput of the link.

A good trade-off between performances and complexity could be to program in advance the transmission scheme, according to weather forecast. A method to predict the attenuation on the link from probabilistic weather forecasts has been developed in the previous section. This Dynamic VCM offers the ability to adapt the link to propagation impairments without the need of a return link. The figure 7 summarizes the steps involved in the implementation of D-VCM strategy. Each day, the control station is informed of the current PEARP forecasts which allow the calculation of the attenuation forecasts A(p, PEARP) for the next communication with the receiving station. This information can be exploited at the level of the control center to build a plan of ModCods for the next visibility period of the receiving station. This solution allows a first degree of adaptation to weather conditions optimizing the capacity of the link C(p,PEARP) to weather forecasts.

IV. RESULTS

A comparison of VCM, ACM and D-VCM strategy is proposed in the next part, first for a fixed Earth space link, second simulating a data download link with a LEO satellite.

The set up and validation of the model requires Ka band attenuation measurements. No data are available for LEO satellites. It prevents the direct validation of the attenuation forecast model in this context and fixed Earth space propagation data are first considered. In a second step, a comparison of D-VCM over VCM is proposed to demonstrate the interest of such method in case of downlink LEO transmission. To assess the performances of the D-VCM strategy, one year of transmission between the geostationary satellite Astra 3B and a receiving station located in Toulouse, France, at latitude 43.5 °N and longitude 1.5 °E has been simulated, using measured attenuation data for the channel. The EIRP of the satellite is assumed to be 50 dBW, the G/T of the receiving station 30 dB/K. A bandwidth of 540 MHz is also considered for the link.

The PEARP forecasts of 2014 have been used as input of the D-VCM decision algorithm. The ACM one was based on the Ka band measurements of the same year assuming no switching losses or margins. The throughputs obtained have been compared to the ones provided by CCM link.

The table [START_REF] Holm | Adaptive Code Modulation Performance and Channel Estimation Tools for Flat Fading Channels[END_REF] illustrates the mean capacities associated to each of these methods for an availability of 99.9 %.

As expected, the best performance is obtained for an ACM strategy. Due to the limits of the weather forecasts, the D-VCM strategy proposed is not able to offer the same level of mean spectral efficiency. It is nevertheless clear that, with an increase of more than 48 %, this technique highly improve the throughput ensured by constant coding and modulation without any prior on the instantaneous state of the propagation channel.

Despite these achievements, the D-VCM is not of great interest in case of a geostationary satellite since either a return is available making it possible to applied ACM, or the absence of communication between the control center and the satellite impose a constant coding and modulation. However, to study the application of D-VCM in the context of a LEO satellite allowed testing the algorithm from end to end.

B. Application to LEO satellite:

It has been shown that opting for D-VCM rather than CCM strategy in context of GEO satellite allows a significant gain in terms of capacity. However, the main application of such method is certainly the case of LEO satellites for which the communication with the control station is occurring once every few hours and for which propagation prediction between those contacts may be needed. A comparison of D-VCM and VCM system is now proposed.

The attenuation probability distribution conditioned to PEARP (figure 4) has been computed from Astra 3B beacon measurements and is thus only representative of the attenuation expected for an elevation angle of 35 °. However, since propagation attenuation increases with decreasing elevation angles, the attenuation probability distribution must be available for the full range of elevation angle.

From the simplification given in equation ( 3) the attenuation probability distribution conditioned to PEARP and to the elevation angle θ of the link is generated.

𝐴(𝜃) = 𝐴(𝜃 = 90°) sin(𝜃) (3)
The elevation time series using for the simulation has been obtained considering a LEO satellite with an altitude of 800 km. Figure 8 illustrates the power margins computed for both systems and the resulting capacities. The forecast of a rain event followed by a lull has been isolated (8-a). A forecasted rain amount equal to 2 mm leads to a power margin ( 8-c) closes to the one obtained for VCM system (8-b). By contrast, in clear sky condition, the power margin of the D-VCM system is substantially lower that the VCM one, resulting in a better capacity (8-e). The peaks visible for both systems are the result of the evolution of the elevation angle of the satellite. On average over the year, D-VCM allows to reach a mean capacity of 821.9 Mbps against 638.8 Mbps for VCM (Tab 3), representing an increase of more than 30 %.

V. DISCUSSION

It has been shown that the implementation of an attenuation prediction system based on probabilistic weather forecasts could constitute an appropriate response to the problematic of the Ka band transmissions operated in the context of Low Earth Orbit satellites. An increase closes to 30% in the mean capacity of a LEO system has in fact been obtained using D-VCM instead of VCM strategy. With some adjustments, these already encouraging results could be improved. First, the model, currently limited to the rain attenuation forecast, could as well include the forecast of cloud, water vapor and oxygen attenuation, that are much more accurately predicted by weather forecasts model and for which the link between the attenuation on the link and the forecasted quantities is easier.

Similarly, a better use of the probabilistic weather forecasts PEARP could be foreseen. As an example, the standard deviation of the members of the ensemble gives us a good estimation of the realizations confidence that can be placed in the forecasted attenuation. This information would allow adapting the choice of ModCods, not only to the weather forecast, but also to its likelihood.

Finally, it would be interesting to study the impact of the spatial resolution of weather forecasts on the performances of the model proposed. Finest resolution allows locating in a more accurate way the rain cells but in a same time and potential error in the weather forecast become thus more penalizing. A tradeoff between resolution and error of forecasts must be found.

The major drawback of the model is probably its lack of flexibility to meet different receiving stations localizations and orbits parameters. In fact, the development of the model is based on beacon measurements, on rain fall measurements and on weather forecasts which characterize the situation around a given receiving station. Furthermore, the rain fall measurements are radar images which often present a limited reliability, in particular due to measurement artifact. A further study will address the use of supervised learning techniques to avoid the use of radar measurements.

Once the model is developed, the only precondition for the implementation of the model is the availability, on a daily basis, of weather forecasts.

Perspectives of such study are numerous and far from being limited to the specific case of LEO satellite. As developed in [START_REF] Biscarini | Optimizing Data Volume Return for Ka-Band Deep Space Links Exploiting Short-Term Radiometeorological Model Forecast[END_REF] and [START_REF] Davarian | Deep space Ka-band link management and Mars reconnaissance orbiter: Long-term weather statistics versus forecasting[END_REF], Deep space Ka Band link also suffers from the propagation impairments and their anticipation, based on weather forecasts, would allow more frequent updates and an increase in the rain rate available on the link during a pass.

Another potential application is the Ka band broadcasting reconfigurable satellites antenna [START_REF] Paraboni | Assessment of Performance of a Ka-Band Boradcasting recondigurable satellite antenna with adaptive mitigation for atmospheric attenuation[END_REF]. Known in advance thanks to weather forecasts, the area with risk of rain could benefit from a greater concentration of power obtained with an appropriate redirection of the power flux.

Figure 1 :

 1 Figure 1: Example of PANTHERE rainfall measurement. The red square represents the study area

Figure 2 :

 2 Figure 2: Rain fall measurement cumulative distribution conditioned to rain forecasts from PEARP

Figure 4 :

 4 Figure 4: Ka band attenuation measurements inverse cumulative distribution conditioned to PEARP rain forecasts. The elevation angle of the link is 35 °.

Figure 5

 5 Figure 5: ACM system principles

Figure 6 :

 6 Figure 6: VCM system principles Figure 7: D-VCM system principle

Tab 2 :

 2 Mean capacities obtained for a target availability of 99.9% and considering a GEO satellite. Comparison of CCM D-VCM and ACM strategies.

Figure 8 :

 8 Figure 8: VCM (a) and D-VCM (b) power margin and resulting capacity time series (c) and (d) computed for a LEO satellite and a target availability of 99.9%. The red lines show the mean capacities obtained for each system.