Olivia Bory Devisme

Frédéric Alicalapa

Denis Genon-Catalot

Pierre-Olivier Lucas De Peslouan

Jean-Pierre Chabriat

Design and implementation of an experimental DC microgrid for energy self-sufficient isolated sites

à la diffusion de documents scientifiques 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.

Olivia BORY DEVISME* , Frédéric ALICALAPA , Denis GENON-CATALOT², Pierre-Olivier LUCAS DE PESLOUAN , Jean-Pierre CHABRIAT *Contact : olivia.bory-devisme@univ-reunion.fr

measurement bench not backed up

  1 ENERGY-lab, Université de la Réunion, 15 Av. René Cassin CS 92003, 97400 Saint-Denis Cedex 9, La Réunion, France 2 LCIS, Université de Grenoble, 50 rue Barthélémy de Laffemas CS 10054, 26902 Valence cedex 9, France electricity produced by solar panels, the energy losses due to transport and the multiple conversions of sources from direct current (DC) to alternating current (AC) and vice versa, considerably reduce the electrical yield. If we go back to Edison's proposal, a DC electrical architecture has many advantages in the case of local energy production since energy losses in conversions and transport are considerably reduced. In order to quantify precisely with measures potential energy gains of a DC-DC architecture in comparison to a classical DC-AC-DC architecture, currently deployed for isolated sites, the objectives of this project are to deploy a DC

	SOLAR PANEL 48 VDC In the context of climate change, the French law on energy transition for green growth (LTECV) has set a target of electrical autonomy for non-interconnected areas such as Reunion Island by 2030. Despite the significant LifePO4 48 VDC Battery increasing use of STORAGE Keywords	LOADS LEDS -FAN -LAPTOP
			EMS SLAVES PoE Drivers
	Power	Power measurements

by generators and then develop a predictive EMS (Energy Management System) including data from environmental indoor and outdoor sensors and PoE (Power Over Ethernet) drivers for each equipments.

:

  DC microgrid, solar energy, storage, isolated sites, PoE (Power Over Ethernet), EMS (Energy Management System). Pilots equipments trough PoE drivers. -Exchanges with a database. -Sets control indicators in user interface.-Makes decisions between storage and consumers.

	MPPT 48 VDC
	SENSORS GHI (W/m²), T°C, H%, Lux, Presence	EMS MASTER MyRIO, CRIO
	Legend :	Communications
		First measurements carried out made it possible to highlight a percentage of 62% of the energy consumption linked to conversions for lighting, within a

-Receives data from sensors. -DC-AC-DC electrical architecture. The perspectives of the project consist in including in the EMS other types of electrical loads and propose a predictive control model based on the available solar resource. This model will also have to integrate hybrid storage (battery, hot water, hydrogen reserve) in order to optimize the battery drain.


domain: spi math

text:
Viktor Kosin

Amélie Fau

Clément Jailin

Benjamin Smaniotto

Thomas Wick

François Hild

A projection-based approach to extend digital volume correlation for 4D spacetime measurements

Keywords: digital volume correlation (DVC), in-situ tests, spacetime analyses, tomography academie-sciences

experiments are generally based on scans reconstructed from a large number of projections acquired under constant deformation of samples. Standard digital volume correlation (DVC) methods are based on a limited number of scans due to acquisition duration. They thus prevent analyses of time-dependent phenomena. In this paper, a modal procedure is proposed that allows time-dependent occurrences to be analyzed. It estimates spacetime displacement fields during the whole loading history.

The spatial modes are based on standard DVC, which is subsequently enriched using projection-based digital volume correlation (P-DVC) to measure the temporal amplitudes. The method is applied to two cases, namely, a virtual experiment mimicking wedge splitting and an actual shear test on a pantographic metamaterial inducing large motions. With the proposed method, the temporal amplitude in the real test was measured for each projection leading to a temporal resolution of one tenth of a second and the analysis of 16,400 time steps. For the proposed algorithm, the sensitivity to the acquisition angle of the sample was investigated and measurement uncertainties were assessed.

Introduction

X-ray computed tomography is an imaging technique that enables to reconstruct volumes of material absorption from series of projections [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF][START_REF] Maire | Quantitative X-ray tomography[END_REF]. The developments of in-situ testing machines, which are placed inside X-ray scanners, give access to images of loaded samples [START_REF] Buffière | In Situ Experiments with X ray Tomography: an Attractive Tool for Experimental Mechanics[END_REF]. The quantitative analysis of sample deformations may be assessed with Digital Volume Correlation (DVC) [START_REF] Bay | Methods and applications of digital volume correlation[END_REF][START_REF] Sutton | Computer Vision-Based, Noncontacting Deformation Measurements in Mechanics: A Generational Transformation[END_REF][START_REF] Buljac | Digital Volume Correlation: Review of Progress and Challenges[END_REF] with the measurement of displacement fields between discrete states. One significant limitation of such a procedure is the acquisition time of each scan. The reconstruction of a digital volume requires the sample to be steady during acquisitions of typically one thousand projections [START_REF] Buffière | In Situ Experiments with X ray Tomography: an Attractive Tool for Experimental Mechanics[END_REF]. Depending on the imaging device, this acquisition may take few minutes up to few hours. This observation forbids in-situ experiments to image time-dependent phenomena such as creep or relaxation that occur when viscoelastic/viscoplastic materials are investigated. Furthermore, relaxation in in-situ experiments may result in motion artifacts. Currently, this is generally avoided by starting the acquisition of projections only after the measured force has stabilized [START_REF] Buffière | In Situ Experiments with X ray Tomography: an Attractive Tool for Experimental Mechanics[END_REF][START_REF] Vargas | Calibration of cohesive parameters for a castable refractory using 4D tomographic data and realistic crack path from in-situ wedge splitting test[END_REF].

To circumvent the aforementioned restrictions, two routes are followed. First, fast scan acquisitions are carried out in synchrotron facilities with extremely bright X-ray beams [START_REF] Ludwig | In-situ Three-Dimensional Microstructural Investigation of Solidification of an Al-Cu Alloy by Ultrafast X-ray Microtomography[END_REF][START_REF] Lhuissier | Continuous characterization by X-ray microtomography of damage during high-temperature deformation of magnesium alloy[END_REF] using high-speed cameras [START_REF] Maire | 20 Hz X-ray tomography during an in situ tensile test[END_REF][START_REF] García-Moreno | Using X-ray tomoscopy to explore the dynamics of foaming metal[END_REF] at the expense of degraded reconstruction quality (due to blur induced by vibrations [START_REF] Maire | 20 Hz X-ray tomography during an in situ tensile test[END_REF]). Another route consists in performing projection-based measurements leveraging the projection temporal sampling instead of their reconstruction. Compared to volume-based registration, the number of projections used to measure displacement fields in projection-based registrations can be reduced by typically three orders of magnitude [START_REF] Leclerc | Projection Savings in CT-based Digital Volume Correlation[END_REF][START_REF] Taillandier-Thomas | Soft route to 4D tomography[END_REF][START_REF] Jailin | In situ CT-scan Mechanical Tests: Fast 4D Mechanical Identification[END_REF][START_REF] Khalili | Discrete Digital Projections Correlation: a reconstruction-free method to quantify local kinematics in granular media by X-ray tomography[END_REF], which allows sub-minute in-situ experiments to be conducted in a lab scanner [START_REF] Jailin | Sub-minute in situ Fracture Test in a Lab CT-scanner[END_REF]. For such approaches, the spatial resolutions of the displacement fields are reduced [START_REF] Leclerc | Projection Savings in CT-based Digital Volume Correlation[END_REF][START_REF] Khalili | Discrete Digital Projections Correlation: a reconstruction-free method to quantify local kinematics in granular media by X-ray tomography[END_REF], a temporal regularization is introduced [START_REF] Jailin | In situ CT-scan Mechanical Tests: Fast 4D Mechanical Identification[END_REF], or mechanical regularization is added to decrease the measurement uncertainties [START_REF] Taillandier-Thomas | Soft route to 4D tomography[END_REF][START_REF] Jailin | In situ CT-scan Mechanical Tests: Fast 4D Mechanical Identification[END_REF].

In projection-based approaches, no standard DVC measurements were considered. This paper presents an enhanced DVC procedure by significantly increasing the temporal sampling of measured displacement fields. Thus, the spatial discretization of the displacement field is totally controlled by standard DVC procedures [START_REF] Buljac | Digital Volume Correlation: Review of Progress and Challenges[END_REF]. By acquiring on-the-fly projections during the loading phases between full scans, it is shown that it is possible to measure temporal amplitudes of the displacement field by projection-based registrations [START_REF] Leclerc | Projection Savings in CT-based Digital Volume Correlation[END_REF]. This method is essentially restricted to the number of frames a tomograph can acquire, which resulted in a temporal resolution of one tenth of a second in the present study. Following a continuous Galerkin finite element discretization in space [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF], a discontinuous Galerkin discretization in time is developed such that the new method is compatible with standard practices of in-situ experiments and their analyses. By computing the displacement field at certain loading steps and not with a projection-based measurement strategy, the presented approach allows for high spatial resolutions and high temporal resolution while keeping the uncertainties low. Conceptually, a 4D spacetime framework is designed in which temporal and spatial components are formulated simultaneously. In modern numerical simulations, spacetime concepts and the corresponding discretizations are currently intensively investigated [START_REF]Space-time methods: Application to Partial Differential Equations[END_REF][START_REF] Wick | Space-time methods: formulations, discretization, solution, goal-oriented error control and adaptivity[END_REF][START_REF] Larsson | Space-time Methods for Time-dependent Partial Differential Equations[END_REF]. In projection-based DVC (P-DVC), 4D concepts were proposed as well by using proper generalized decomposition concepts [START_REF] Jailin | Fast 4D tensile test monitored via X-CT: Single projection based Digital Volume Correlation dedicated to slender samples[END_REF][START_REF] Jailin | Fast four-dimensional tensile test monitored via X-ray computed tomography: Elastoplastic identification from radiographs[END_REF]. The following study introduces an alternative route constructing spatial modes via DVC, and then temporal modes via P-DVC. To this end, one of the goals is to formulate such 4D concepts within a spacetime setting.

The outline is as follows. In Section 2, the principle of the proposed procedure is presented by introducing the corresponding variational formulations of DVC and Projection-based DVC (P-DVC) methods, and formulating the sequential solution algorithm. The method is tested on a synthetic experiment and applied to a shear test on a pantographic metamaterial in Section 3.

Projection-enhanced DVC

Spacetime framework

The general aim of kinematic measurements is to determine how some material domain Ω ⊂ R 3 (or region of interest) deforms during time interval [0, t end ], which forms the spacetime cylinder [START_REF]Space-time methods: Application to Partial Differential Equations[END_REF][START_REF] Wick | Space-time methods: formulations, discretization, solution, goal-oriented error control and adaptivity[END_REF][START_REF] Larsson | Space-time Methods for Time-dependent Partial Differential Equations[END_REF]. One is seeking a displacement field u : Ω × (0, t end ) → R 3 that describes the change of the reference domain Ω to the deformed configuration Ω(t ), such that

Ω(t ) = x = x + u( x, t ) : x ∈ Ω, t ∈ [0, t end ] , (1)
where x denotes the position of any material point in the reference configuration Ω, and x the position at time t . The displacement at spatial location x ∈ Ω(t ) is defined by

u(x, t ) := x -x = u( x, t ), (2)
which completes the notations needed to describe the motion in Lagrangian and Eulerian frameworks.

A classical approach to measure u requires Ω and Ω(t ) to be represented by a series of 3D images [START_REF] Bay | Methods and applications of digital volume correlation[END_REF][START_REF] Buljac | Digital Volume Correlation: Review of Progress and Challenges[END_REF][START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF]]

I : D × [0, t end ] → R ≥0 , (3)
where D ⊂ R 3 is the field of view, with Ω(t ) ⊂ D for all t ∈ [0, t end ]. In practice, an image I t (x) := I (x, t ) cannot be acquired instantaneously at time t (Figure 1). To reconstruct such an image by X-Ray computed tomography, hundreds if not thousands, of radiographs (i.e. projections) are needed over a full revolution (or turn) of the sample with respect to the detector [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. Conveniently, standard in-situ tomography experiments are designed to have phases when the loading is kept steady and the sample is rotated for one turn to acquire enough projections for a full 3D reconstruction of volumes in deformed configurations [START_REF] Buffière | In Situ Experiments with X ray Tomography: an Attractive Tool for Experimental Mechanics[END_REF]. These 3D images are no longer an instantaneous representation of the configuration Ω(t ). Instead, they are reconstructed over the time needed for a turn T (Figure 1). To emphasize the different time scale, an image reconstructed from a set of radiographs acquired over a turn T is denoted by I T . For displacement fields describing the 3D image I T of the deformed configuration, capital letters will be used i.e.

U( x, T ).

To measure the spacetime displacement field u, separable variables are assumed. Hence, u is decomposed as

u( x, t ) = α(t ) W( x), (4)
with α : (0, t end ) → R, and W : Ω → R 3 . Let T f and T g be turns in [0, t end ] with T f < T g , U( x, T f ) and U( x, T g ) the corresponding displacement fields, respectively, which are measured using standard DVC methods. Then, the space component of the displacement field is defined from DVC analyses as follows

W( x) = U( x, T g ) -U( x, T f ). (5)
Having W( x) assessed, only the temporal amplitude α(t ) is left unknown.

Remark 1. The freedom to choose T f and T g allows the displacement field to be computed at any point in time, before, after or even during reconstruction intervals. For example, when a sample is tested until failure, performing a final full scan is difficult, whereas it is simple to perform scans before and only acquire projections until failure to measure pre-failure displacement fields.

A full acquisition of the initial state can also be challenging since the sample may need to be preloaded to ensure minimal motion for the first scan. Here one can exploit the fact that measurements are possible when t < t T f .

Digital volume correlation

As for most of the algorithms used for digital volume correlation, the conservation of gray levels of I is assumed. This hypothesis states that the gray level of each material point x is considered constant in time. Hence, for each image I T g ( x) := I ( x, T g ) there is a displacement field U( x, T g ) such that

I T f ( x) = I T g ( x + U( x, T g )) ∀ x ∈ Ω, T g ∈ [0, t end ], (6)
where I T f is the image of the reference state, i.e. U( x, T f ) ≡ 0. Then, by considering another image I T g ( x) and applying decomposition (4) with condition (5), it follows that U( x, T g ) = W( x).

Consequently, assumption (6) reads as follows. There is a displacement field W such that

I T f ( x) = I T g ( x + W( x)) =: I T g , W ( x) ∀ x ∈ Ω. (7)
One seeks for a solution of the following minimization problem [START_REF] Fedele | Global 2D digital image correlation for motion estimation in a finite ele-ment framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF].

Find W ∈ V( Ω) := L 2 ( Ω) such that J A ( W) = inf V∈V( Ω) J A ( V) with J A ( V) := 1 2 ∥I T f -I T g , V ∥ 2 L 2 ( Ω) . (8)
Here, W is a solution to Equation ( 8) only if it fulfills the following first-order optimality condition

A( W)( V) := Ω (I T g , W -I T f )( ∇ ∇ ∇I T g , W • V ) d x = 0 ∀ V ∈ V( Ω). (9)
This is a nonlinear problem, and a Gauss-Newton scheme (see e.g. Ref. [START_REF] Deuflhard | Newton Methods for Nonlinear Problems[END_REF]) will be utilized for the minimization problem ( 9), the Hessian is needed

A ′ ( W)(δ W, V) = Ω ( ∇ ∇ ∇I T g , W • δ W)( ∇ ∇ ∇I T g , W • V) + (I T g , W -I T f ) δ W • ( ∇ ∇ ∇ 2 I T g , W ) • V d x. (10)
This solution strategy is standard and explained in more details in Ref. [START_REF] Fedele | A regularized, pyramidal multi-grid approach to global 3Dvolume digital image correlation based on X-ray micro-tomography[END_REF].

Projection-based digital volume correlation

In the same spirit as the presented DVC framework, the underlying hypothesis of P-DVC is again the conservation of gray levels. Contrary to Equation ( 6), the problem is not stated in the Lagrangian framework, but in the Eulerian system [START_REF] Leclerc | Projection Savings in CT-based Digital Volume Correlation[END_REF]. With I T f (x) := I (x, T f ), there is a displacement field u(x, t ) such that

I T f ,u(t ) (x) := I (x -u(x, t ), T f ) = I (x, t ) ∀x ∈ Ω(t ), t ∈ [0, t end ]. (11)
Instead of assuming that there is another 3D image, only one projection is considered to be available for any considered time t . Let Π θ be the linear projection operator of the underlying X-ray setup (Figure 2), which maps any 3D gray level image I t (x) := I (x, t ) defined on the volume domain D to a 2D gray level image Π θ [I t ] defined on the detector domain Ξ at angle θ(t

) ∈ [0, 2π].
Then, assuming Equation ( 11) holds, it follows that

Π θ(t ) [I T f ,u(t ) ](r) = Π θ(t ) [I t ](r) ∀r ∈ Ξ, t ∈ [0, t end ], (12)
where r denotes a pixel coordinate defined on the detector plane Ξ. The projections acquired during an experiment are defined as Let ω(t ) := Π θ(t ) [Ω(t )] ⊂ Ξ be the domain of the projection of Ω(t ) at time t . Then, the gray level conservation assumption for the P-DVC problem is expressed as follows. There is a timedependent amplitude α(t ) such that

p t (r, θ(t )) := Π θ(t ) [I t ](r). (13)
Π θ(t ) [I T f ,α,W ](r) = p t (r, θ(t )) ∀r ∈ ω(t ), t ∈ [0, t end ], (14)
with

I T f ,α,W (x, t ) := I T f (x -α(t )W(x)). ( 15
)
Since W is computed using DVC and W(x) = W( x), only α(t ) needs to be computed to assess u(x, t ). Consequently, the corresponding minimization problem consists in finding α ∈ Q := L 2 ((0, t end )) such that

J B (α) = inf β∈L 2 ((0,t end )) J B (β) with J B (β) := 1 2 ∥Π θ(t ) [I T f ,β,w ] -p t ∥ 2 L 2 ((0,t end );L 2 (ω(t ))) , (16)
where the norm is defined as

∥ • ∥ 2 L 2 ((0,t end );L 2 (ω(t ))) := (0,t end ) ω(t ) • 2 dr dt . ( 17
)
This time, the corresponding optimality condition is given by

B (α)(β) = (0,t end ) ω(t ) β(Π θ(t ) [I T f ,α,W ] -p t )Π θ(t ) [∇ ∇ ∇I T f ,α,w • W] dr dt ! = 0, ∀β ∈ L 2 ((0, t end )). ( 18
)
To use again a Gauss-Newton scheme to solve this nonlinear problem, the Hessian is computed

B ′ (α)(δα, β) = (0,t end ) ω(t ) δαβΠ θ(t ) [∇ ∇ ∇I T f ,α,W •W] 2 +δαβ(Π θ(t ) [I T f ,α,W ]-p t )Π θ(t ) [W•(∇ ∇ ∇ 2 I T f ,α,W )•W]dr dt . (19)
Remark 2. In theory, the projected domain ω(t ) and the projection p t represent a projection of an image I t of the deformed volume Ω(t ). But I t is not known. Further, Ω(t ) is unknown as well, and thus ω(t ). This feature has to be taken into account in the devised algorithm.

Solution algorithm

Summarizing the previous sections, two minimization problems are solved successively. For both minimization problems, a Gauss-Newton scheme is considered. The procedure is described in Algorithm 1. Solve

A ′ ( W l )(δ W, V) = -A( W l , V) ∀ V ∈ V (20) 5:
W l +1 ← δ W + W l 6: end while

7: W(x) ← W( x) 8: while ∥u l -u l -1 ∥ > T OL u do ▷ P-DVC minimization problem 9:
Approximate Ω(t ) by

Ω l (t ) ← x + u l ( x, t ) : x ∈ Ω; u l ( x, t ) = α l (t ) W( x) . ( 21
)
10:

ω l (t ) ← Π θ(t ) [Ω l (t )] 11: Solve B ′ (α l )(δα, β) = -B (α l )(β) ∀β ∈ Q (22) 12:
α l +1 ← α l + δα 13:

u l +1 ( x, t ) ← α l +1 (t ) W(x)
14: end while Remark 3. In practice, the Hessians A ′ and B ′ given in Equations ( 10) and ( 19) are not exactly computed as they were defined, but only approximated by their first term. First, the second term is not always positive and therefore the Hessian itself could get indefinite. This could lead to a search direction that does not follow a descent direction. Second, it includes the Hessian of the current image, which itself is quite expensive to compute. To summarize, the two approximate Hessians read

A ′ ( W)(δ W, V) ≈ Ω ( ∇ ∇ ∇I T g , W • δ W)( ∇ ∇ ∇I T g , W • V) d x (23)
and

B ′ (α)(δα, β) ≈ (0,t end ) ω(t ) δαβΠ θ(t ) [∇ ∇ ∇I T f ,α,W • W] 2 dr dt . (24)
A detailed discussion about this point and other approximation approaches of the Hessian matrix is found in Ref. [START_REF] Neggers | On image gradients in digital image correlation[END_REF].

Discretization

The variational formulations [START_REF] Larsson | Space-time Methods for Time-dependent Partial Differential Equations[END_REF] and [START_REF] Jailin | Fast four-dimensional tensile test monitored via X-ray computed tomography: Elastoplastic identification from radiographs[END_REF] in the previous sections are not computationally accessible. Therefore Galerkin-type methods are introduced to solve the minimization problems in a numerical manner. In space, a finite-element discretization is defined on a tetrahedral discretization T h of Ω with elements K and piecewise linear functions

P 1 V h ( Ω) := V h ∈ C ( Ω) ∩ L 2 ( Ω) : V h | K ∈ P 1 for all K ∈ T h . ( 25
)
Let { Φ i } be the basis of V h ( Ω), then each element

V h ∈ V h ( Ω) is of the form V h ( x) = i a i Φ i ( x). (26)
Then, the Galerkin ansatz (i.e. approximating partial differential equations in weak form using a finite-dimensional function space [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) is used to replace Equation ( 20), and solve instead the following linear system

M l i j a j = m l i , (27)
with

M l i j = Ω ( ∇ ∇ ∇I T g , W l h • Φ j )( ∇ ∇ ∇I T g , W l h • Φ i ) d x and m l i = Ω (I T f -I T g W l h )( ∇ ∇ ∇I T g , W l h • Φ i ) d x. (28)
In time, one wants to analyze each time step or more precisely each projection p t k separately.

Hence, the interval (0, t end ) is subdivided into sub-intervals (t k -d t /2, t k + d t /2) for k = 1, ..., n t , and piecewise constant basis functions are selected. This decomposition results in the following (discontinuous) finite element space

Q h ((0, t end )) := β h ∈ L 2 ((0, t end )) : β h (t )| (tk-dt/2,tk+dt/2) = const., k = 1, ..., n t . ( 29
)
As a basis of Q h , {ψ k } is chosen as

ψ k (t ) =        1 t ∈ (t k -d t /2, t k + d t /2) , 0 otherwise. ( 30
)
This setting allows each element β h ∈ Q h to be written as

β h (t ) = k b k ψ k (t ). (31)
Last, instead of solving Equation ( 22), the following linear system is considered

N l i j b j = n l i ( 32
)
with

N l i j = (0,t end ) ω l (t ) ψ j ψ i Π θ(t ) [∇ ∇ ∇I T f ,α l h ,W h • W h ] 2 dr dt n l i = (0,t end ) ω l (t ) ψ i (p t -Π θ(t ) [I T f ,α l h ,W h ])Π θ(t ) [∇ ∇ ∇I T f ,α l h ,W h • w h ] dr dt . ( 33
)
Remark 4. Assembling N i j and n i defined in Equation ( 33) is not straightforward. In a standard finite element code, one would simply choose some quadrature schemes. Without any preliminaries, this would require that one can select arbitrary quadrature points t q in each interval

(t k -d t /2, t k + d t /2
). Unfortunately one cannot evaluate p t at any point in time, but only at the grid points t k , as p t is given by projections acquired during an experiment. In the present case, this issue is bypassed by having constant shape functions and using the midpoint rule for integration. Matrix [N] then becomes diagonal, such that each degree of freedom is solved independently and Equation (33) reduces to

b i = n l i N l i i . ( 34
)
One can freely choose the quadrature points when integrating over ω(t ) or Ω by interpolating over the grid.

Application

In this section the proposed method is tested. First, only the projection-based Gauss-Newton scheme (Equations ( 21) and ( 22)) is applied to a synthetic experiment, where the displacement field W is predefined. Then, the complete algorithm is run on a real experiment. The ASTRA toolbox [START_REF] Van Aarle | The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography[END_REF][START_REF] Van Aarle | Fast and flexible X-ray tomography using the ASTRA toolbox[END_REF] was used for the cone beam projection operator Π and the DVC calculations were performed with the FE-based DVC library Correli 3.0 [START_REF] Leclerc | Correli 3.0[END_REF].

Synthetic experiment

For a first assessment of the proposed algorithm, a study of a predefined 3D image I T f of size 120 × 90 × 200 vx (Figure 3) was performed. given by a uniform value

I T f | Ω ≡ 5 • 10 -3 . ( 35
)
The shape and motion of this synthetic case were inspired by wedge splitting tests [START_REF] Vargas | Calibration of cohesive parameters for a castable refractory using 4D tomographic data and realistic crack path from in-situ wedge splitting test[END_REF], but for the sake of simplicity the sample did not rotate. The displacement field was chosen as

u( x, t ) = α(t ) W( x), (36)
where

W( x) = 10 L z (z -30) 2 L x (x -60.5) e x , (37)
and

α(t ) = 3τ 2 -2τ 3 with τ = t t end , ( 38
)
with all spatial numbers expressed in voxels. The analytical displacement field W was chosen such that it was perfectly represented by the selected spatial finite element discretization (Figure 3 To validate the accuracy and speed of the computation of the temporal amplitude α with the proposed method, the absolute difference of the solution α and the computed amplitude α h is presented in Figure 5(a). Since each amplitude was computed by minimizing the L 2 -norm of the residuals ρ t , the convergence of the root mean square (RMS) of the residuals is displayed in

|α(t ) -α h (t )| τ = 0 τ = 1/4 τ = 2/4 τ = 3/4 τ = 1 x -x

∥Π θ(t

) [I f ,α h ,w ] -p t ∥ L 2 (ω(t )) τ = 0 τ = 1/4 τ = 2/4 τ = 3/4 τ = 1 x -x (b)

Shear test on pantographic metamaterial

The proposed algorithm was also studied on an actual experiment. An in-situ shear test on a pantographic metamaterial, which was made of Polyamid PA2200 powder [START_REF] Dell'isola | Advances in Pantographic Structures: Design, Manufacturing, Models, Experiments and Image Analyses[END_REF], was performed.

The geometry of the sample is shown in Figure 6(b) and the pantographic structure was of size 105 mm × 35 mm × 45 mm. The sample, which was mounted horizontally in the in-situ tensiontorsion-compression (TTC) testing machine (Figure 6(a)), was scanned in the North Star Imaging X50+ tomograph of LMPS at a binned definition (i.e. 1944 × 1536 pixels, see Table 1). The shear displacement was induced by clamping one side and pushing the other one down, such that a maximum stroke of 30 mm was reached at the end of the test (Figure 6(c)). The clamps were made of (3D printed) ABS and are visible in the reconstructed volumes as well as in the projections.

Therefore, one has to take the motion of the clamps into account to have good results for the projection-based measurements. After flat-field corrections, the projections were further binned by a factor of four such that their final definition was 486 × 344 pixels, and the reconstructed volumes had a resolution of 372 µm per voxel. The reconstructions were performed with the simultaneous iterative reconstruction algorithm (SIRT) for cone beams [START_REF] Van Aarle | Fast and flexible X-ray tomography using the ASTRA toolbox[END_REF]. In the reference configuration, the pantographic structure was in contact with the support. This contact zone was not ideally meshed such that there the gray level assumption could not be fulfilled anymore (i.e. I T f , u h and Π θ(t ) [I T f , u h ] were only approximated). This phenomenon led to oscillations in α h .

To exclude these zones, Π θ(t ) [Ω] was not chosen as the region of interest, but a slightly smaller subdomain of Π θ(t ) [Ω] (Figure 9(a)). The in-situ experiment protocol was devised as follows:

• The sample was turned continuously with a rotation speed of 3 • /s during two revolutions, and 2,400 projections were simultaneously acquired at a rate of 10 fps (i.e. n = 1, 200, and d t = 0.1 s). Scan 0 corresponds to the second turn.

• Then, the sample was loaded with a constant vertical stroke rate of 56 µm/s until a vertical displacement of 10 mm was reached while rotating at 3 • /s for two more turns acquiring 1,200 projections per turn. Scan 1 corresponds to turn 6.

• The procedure was repeated until a maximum vertical displacement of 30 mm was reached considering 2 turns for each loading phase, and 2 turns for each stage with steady stroke.

The finite element mesh had 7,234 elements with an average volume of 18 3 vx (Figure 7). For the analysis of the experiment, algorithm 1 was used with T OL W = 10 -3 vx and T OL u = 10 -3 vx. As the solution strategy was decoupled (Remark 4), the minimization was performed for each projection separately. This separation allows the solution of the current projection to be used for the following one, such that convergence was reached after 7 iterations on average. The computation of each amplitude α took about one minute. In Figure 8 (a), the measured maximum vertical displacement is presented. One clearly sees the planned loading history with steady and loading phases. Due to continuous rotation of the sample, the projections were acquired at an increment of 0.3 • , which resulted in visible fluctuations because not every angle of the sample was well suited for the registration procedure. For example, projections at some angles have more artifacts than others and the norm of the projected gradient is much lower such that they became more sensitive to acquisition noise.

Interestingly, the fluctuations did not increase over time and remained very low. The standard deviations of the amplitudes α during the turns for which the stroke was steady were equal to 7.4 × 10 -3 , 8.7 × 10 -3 , 7.7 × 10 -3 and 7.1 × 10 -3 ), the corresponding standard deviations of the maximum displacement in the loading direction (i.e. u z ) were 0.20 vx, 0.23 vx, 0.21 vx and 0.19 vx, respectively.

Comparing the RMS of the projection residuals

ρ t = Π θ(t ) [I t ,α, W ] -p t (39)
with the eigenvalue of the sensitivity B ′ (Equation ( 19)) at convergence (Figure 8(b)), it is observed that the angles for which the RMS residuals were high were associated with low values of the sensitivity and vice versa. This result indicates that the choice of the acquisition angle has an impact on the quality of the solution. In Figure 9(a), one can see that the clamps are in front and behind the sample resulting in projections low contrast. Yet, the residuals were low with respect to the dynamic range of the original radiographs even for angles that were not optimal (Figure 9(b-e)). The remaining residuals, visible in the zoomed view of scan 0 and scan 3 (Figure 10), show that the highest values are at the end of the beams of the pantograph. In these areas the macromesh, which was used for the analysis (Figure 7), was not able to match the kinematic details of the beams. Similar residuals were observed in the DVC analyses (for the same reasons). By acquiring projections not only during a scan but also during the loading phase, it was possible to increase the temporal discretization of 4D measured displacement fields by more than three orders of magnitude (i.e. from 4 scans to 16,400 time steps) with projection-enhanced DVC.

Conclusion

Since the choice of the acquisition angle had an impact on the quality of the solution, methods for choosing optimal angles are highly interesting and should be investigated. Yet, having already accurate results and the possibility of measuring a displacement of a volume at a speed of one tenth of a second indicated that the proposed methodology may help to measure real timedependent phenomena such as, for example, viscous effects or crack propagation.

Conflicts of interest

The authors declare no competing financial interest.

Figure 1 .

 1 Figure 1. DVC and P-DVC temporal frames. In present case, n radiographs are acquired per turn. The temporal delay between two radiographic acquisitions is d t .

Figure 2 .

 2 Figure 2. Schematic view of X-Ray Computed Tomography with a point source. A projection operator Π θ maps the volume I t defined on the domain D onto a 2D gray level image p t defined on the detector plane Ξ.

Algorithm 1 P 2 : 3 :

 123 -DVC enhanced DVC algorithm 1: Choose an initial guess W 0 Choose an initial guess α 0 while ∥δ W∥ > T OL W do ▷ DVC minimization problem 4:

Figure 3 .

 3 Figure 3. Sections of the phantom used in the synthetic experiment. The mesh is composed of 198 T4 elements whose mean volume is equal to 18 3 vx.

Figure 4 .

 4 Figure 4. Artificial projections of the 3D image I τ at dimensionless times (a) τ = 0, (b) τ = 1/4, (c) τ = 1/2, (d) τ = 3/4, (e) τ = 1.

Figure 5 (

 5 Figure 5(b). Even though no specific initialization was chosen for all projections (i.e. α 0 (t k ) = 0), the algorithm converged extremely fast, i.e. exponential convergence rates were reached for all time steps (Figure5(a)), and with very low residuals (Figure5(b)) essentially due to interpolation errors. The Newton solver converged after a maximum of 10 iterations to a tolerance of 10 -9 in terms of amplitude corrections. A lower tolerance was not reached, because interpolation errors started to dominate the residual. For the last projection, this result indicates that the algorithm was not only fast but also robust for larger displacements (up to 10 vx in the present case).

Figure 5 .

 5 Figure 5. Convergence of the amplitude α (a) and the projection (b) as functions of Newton iterations. The (cyan) dashed lines show that exponential convergence is reached.

Figure 6 .

 6 Figure 6. Sample inside the testing machine (a). Reconstructed volumes for scans 0 (b) and 3 (c).

Figure 7 .

 7 Figure 7. Finite element mesh used for the DVC analyses of the pantograph. The mean element volume is 18 3 vx.

Figure 8 .

 8 Figure 8. (a) Measured maximum vertical displacement u z during the whole test. (b) Comparison of the angular eigenvalue (EV) of the Hessian matrix [N] and the RMS residual with respect to the angular position averaged over the first 6 revolutions (scans 0 to 1 analysis).

Figure 9 .

 9 Figure 9. (a) Flat-field corrected projections with the region of interest ω (blue) at 0 • (top) and 90 • (bottom) angles. (b)-(e) Final residuals ρ for scans 0 to 3 for the 0 • (top) and 90 • (bottom) angles. The projections and final residuals are scaled by the maximum gray level in the region of interest ω.

Figure 10 .

 10 Figure 10. Zoomed view of the residuals ρ for scans 0 (a) and 3 (b). At the top are the 0 • views and at the bottom the 90 • views. The gray level color bar is scaled to the full dynamic range of the residuals inside the region of interest.

Figure 11 .

 11 Figure 11. (a) Measured maximum vertical displacement u z for scans to 1 compared to the prescribed stroke. (b) Difference between both signals for each turn.

A

  new spacetime developed to measure 4D displacement fields by combining DVC and P-DVC. Using a spatiotemporal separation of variables, the computation of the spatial component (via DVC) was based on reconstructed volumes, and that of the temporal component (via P-DVC) on the projections acquired on-the-fly during the whole experiment.

Viktor Kosin, Amélie Fau, Clément Jailin, Benjamin Smaniotto, Thomas Wick and François Hild

Acknowledgments

The financial support of the French-German University through the French-German Doctoral college "Sophisticated Numerical and Testing Approaches" (CDFA-DFDK 19-04) is acknowledged. This project was also partially supported by the French "Agence Nationale de la Recherche" through the "Investissements d'avenir" program (ANR-10-EQPX-37 MATMECA Grant). Discussions are acknowledged within the framework of the International Research Training Group on Computational Mechanics Techniques in High Dimensions GRK 2657 funded by the German Research Foundation (DFG) under Grant Number 433082294.

Funding. Financial support from the French-German University through the French-German Doctoral college "Sophisticated Numerical and Testing Approaches" (CDFA-DFDK 19-04) is acknowledged. This project was also partially supported by the French "Agence Nationale de la Recherche" through the "Investissements d'avenir" program (ANR-10-EQPX-37 MATMECA Grant). This article is a draft

Dedication

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Appendix: Hardware parameters

The hardware parameters of the in situ setup are gathered in Table


Nícolas De

Araújo Moreira

Vasconcelos Rubem
email: rubem070@gmail.com

Pacelli

Yuri Carvalho

Barbosa Silva

Tarcísio Ferreira Maciel

Ingrid Simões
email: ingrid@gtel.ufc.br

João César

Moura Mota
email: mota@gtel.ufc.br

Cerine Hamida
email: chamida@mit.edu

Rodrigo Zambrana Prado
email: rodrigo.zambrana@weatherforce.org

Emilie Caillault
email: ecaillault@gmail.com

Modeste Kacou
email: modeste.kacou@ird.fr

Marielle Gosset
email: marielle.gosset@ird.fr

Convolutional Long-Short-Term Memory Networks (ConvLSTM) for Weather Prediction using Radar and Satellite Images ⋆

Keywords: Artificial intelligence, machine learning, ConvLSTM, image processing, weather forecasting

Artificial Intelligence techniques, mainly machine and deep learning ones, are becoming the most common approach for data prediction. In this context, using these techniques instead of approaches based on classical statistics has shown interesting and important contributions to weather prediction. The present paper discusses the prediction of rainfall and clouds direction based on a sequence of 10 and 14 frames of radar images with a loss inferior to 0.06. Two different Convolutional Long-Short Term Memory Networks configurations were tested and this work presents the estimated frames resulting from these algorithms and presents comparisons between them, real data, and with the performance of other works. The results show that these algorithms can be suitable for short-term weather forecasting.

INTRODUCTION

Climate change is one of the most important concerns of the 21st century and it is a well-known fact that weather modeling and prediction are a mathematical and computational challenge. Meteorological phenomena are complex dynamic systems and involve different fields of physics and mathematics.

According to [START_REF] Shi | Convolutional lstm network: A machine learning approach for precipitation nowcasting[END_REF] and [START_REF] Koutník | Convcast: An embedded convolutional lstm based architecture for precipitation nowcasting using satellite data[END_REF], nowcasting, i.e. short-term prediction of rainfall intensity in a local region, involves non-uniform modeling of meteorological phenomena over time. More precisely, rainfall is a continuous but intermittent random process usually recorded as a cumulative quantity over constant time intervals as is declared by [START_REF] Ramesh | Temporal modelling of short-term rainfall using cox process[END_REF]; [START_REF] Poornima | Prediction of rainfall using intensified lstm based recurrent neural network with weighted linear units[END_REF] while precipitation fields are complex and irregular according to [START_REF] Rasmussen | Stochastic simulation of high-resolution daily precipitation using Gaussian Markov Random Fields[END_REF]. All these features increase the difficulty of a precise nowcasting.

⋆ The authors would like to thank Institut de Recherche pour le Développement for the financial support.

Rainfall has an influence on ecosystems, agriculture, energy production, and even tourism. It is an important resource for human life: long drought periods, monsoons, and floods have a significant impact on the economy and may lead to real disasters. Hence, meteorological prediction can help public managers and monitoring organizations to make decisions to prevent such major impacts as is shown in [START_REF] Poornima | Prediction of rainfall using intensified lstm based recurrent neural network with weighted linear units[END_REF]. The improvement in predicting the behavior of hurricanes, for example, can impact directly early warning systems in regions affected by this kind of event and the anticipation of the behavior of hurricanes can help local governments and emergency teams in planning mitigation schemes for reducing impacts, avoiding bigger tragedies, and loss of lives.

Until the end of the 20th century, the modeling and prediction of meteorological events were based on classical probability and statistics methods, such as Coxer process and Neyman-Scott process as is shown in [START_REF] Smith | Statistical inference for point process models of rainfall[END_REF]; [START_REF] Ramesh | Temporal modelling of short-term rainfall using cox process[END_REF] (usually used to model frontogenesis processes), Gaussian Random Field (GRF, used for representing precipitation in time and space) as is shown in [START_REF] Liu | Advances in gaussian random field generation: a review[END_REF]; [START_REF] Rasmussen | Stochastic simulation of high-resolution daily precipitation using Gaussian Markov Random Fields[END_REF], fractals, as is discussed by [START_REF] Sivakumar | Fractal analysis of rainfall observed in two different climatic regions[END_REF]; [START_REF] Breslin | Fractal dimensions for rainfall time series[END_REF]; [START_REF] Lovejoy | Fractal properties of rain and a fractal model[END_REF]; [START_REF] Pathirana | Fractal modeling of rainfall: downscaling in time and space for hydrological applications[END_REF], etc. For instance, the authors in [START_REF] Lovejoy | Fractal properties of rain and a fractal model[END_REF] state that the projection areas on the Earth of clouds and rain have shapes whose boundaries are fractal curves while the spatiotemporal structure of rain has frequently hyperbolically distributed features (random variables), which also characterizes a fractal. The authors in Gyasu [START_REF] Gyasu Agyei | Interpolation of daily rainfall networks using simulated radar fields for realistic hydrological modelling of spatial rain field ensembles[END_REF] used GRF for interpolating daily rainfall using simulated radar fields for realistic hydrological modeling of spatial rain fields.

However, the use of artificial intelligence techniques is becoming more and more popular in environmental monitoring and meteorological survey, mainly in weather forecasting and predicting the behavior of meteorological phenomena such as hurricanes and typhoons as discussed in [START_REF] Ruttgers | Prediction of a typhoon track using a generative adversarial network and satellite images[END_REF]. These methods were shown to be simple, feasible, and present a good level of accuracy. In this context, radar and satellite images play an important role in addition to historical data.

The objective of the present work is to verify the performance of Convolutional Long-Short-Term Memory (Con-vLSTM) networks for digital image processing for predicting the short-term path of clouds and rainfall obtained from a sequence of frames of meteorological radar and/or satellite images.

The remainder of this paper is organized as follows. Section II presents a literature review and the state of the art related to weather prediction and other applications of artificial intelligence, machine and deep learning, statistics, probability, and stochastic processes on meteorology, as well as an overview of key concepts related to this topic. Section III presents in detail how Convolutional Networks, LSTM, and ConvLSTM work and the methodology used. Section IV presents and discusses the results. Finally, in Section V, the paper presents the conclusions and some perspectives for future works.

LITERATURE REVIEW AND STATE OF THE ART

The authors in [START_REF] Shi | Convolutional lstm network: A machine learning approach for precipitation nowcasting[END_REF] Image series prediction using convolutional LSTM applied to precipitation nowcasting is also discussed in [START_REF] Wu | Sequential images prediction using convolutional LSTM with application in precipitation nowcasting[END_REF].

Some other recent results on deep learning-based approaches for precipitation nowcasting are summarized in [START_REF] Gao | Deep learning and the weather forecasting problem: Precipitation nowcasting[END_REF], where the authors initially describe mathematically the forecasting problem as a spatiotemporal problem, present a literature review on the state of the art and introduce a systematic benchmark for evaluating the performance of the discussed works. Other approaches related to video frame prediction are discussed in [START_REF] Hong | Psique: next sequence prediction of satellite images using a convolutional sequence-to-sequence network[END_REF].

THEORETICAL INTRODUCTION

In this section, the key concepts of Convolutional, Recurrent, Long-Short Term Memory, and ConvLSTM networks will be introduced.

Convolutional Neural Networks (CNNs)

Convolutional Neural Networks (CNNs), or ConvNets, are feed-forward deep learning algorithms designed to work with bi-dimensional data, particularly images, that present many advantages such as less demanding pre-processing and image reduction without loss of significant features.

It is used mainly for image processing due to its ability to learn a large number of filters in parallel. It is usually applied to predictions and image classification.

CNNs are multilayer architecture composed by connected nodes. Each layer is composed of a two-dimensional plan containing many different neurons. The input layer receives the initial data and each neuron extracts the local features using a convolution kernel. The training process does a continuous modification of this convolution kernel weights. The convolution operation extracts high-level features such as edges. The network is composed by one or more pairs of convolutional layers and a max pooling layer. The max pooling enhances the effect of translation invariance and tolerance to deviations, reduces the spatial size (dimensionality reduction) of the convolved feature, reducing the demanded computational power, and suppresses the noise. The max pooling returns the maximum value from the portion of the image covered by the kernel.

The addition of a fully-connected layer is a simple way to learn non-linear combinations of the high-level features. The input image is flattened into a column vector and then fed to a feed-forward neural network and backpropagated in every iteration of training. Over a series of epochs, the model can distinguish dominating and low-level features and classify them using softmax classification.

3.2 Long-Short-Term Memory Networks (LSTM)

Recurrent Neural Networks are networks with a feedback loop that allows information to persist. So, from an architectural point of view, it is a kind of densely connected neural network with the output layer on the recurrent neural network connected to itself and with the output of the hidden layer passing through a delay block. We can create a chain of these blocks with multiple copies of these structures. Denoting by U and V the matrices connecting inputs and recurrent outputs, and with x t , h t-1 denoting an input and the output of the previous cell, we have:

h t = σ(Ux t + Vh t-1 ).
(1) Long Short-Term Memory (LSTM) networks are a type of recurrent neural network used in deep learning and widely applied tn time-series analysis. The main feature which distinguishes it from other types of neural networks is the introduction of a forget gate that controls which states are remembered or forgotten.

A gate is a structure that adds, removes, or transmits information. They are composed of a layer of a sigmoid neural network (here denoted by σ) and a Hadamard product (denoted by •). The sigmoid layer outputs a number within the interval [0, 1], where 0 denotes not allow the transmission of any information, and 1 denotes that all information is re-transmitted. Each cell is composed of five parts:

• Cell state (c t ): represents the internal memory of the cell that stores short and long-term memories. • Hidden state (c h ): decides to retain short and/or longterm information on cell state to make the prediction; • Input gate (i t ): conditionally decides which input values will be updated on memory state, or, in other words, decides which information to forget and which will be set on the input of cell state; • Forget gate (f t ): decides how much information of current input and from previous cell state is forwarded to current cell state; • Output gate (o t ): conditionally decides which will be the output according to input and memory block.

So, the first step of LSTM is to decide which information will be ignored and then which new information will be stored.

The previous cell c t-1 is updated into a new cell c t and then the old state is multiplied by a forgetting factor f t .

A new candidate value is generated, so it is necessary to decide what will be transmitted to output based on the cell state by applying a sigmoid layer. The tanh layer puts the value on the interval [-1, 1]. Denoting U g , V g the weights for input and previous cell output, respectively, and b g an input bias, we have:

g = tanh(b g + x t U g + h t-1 V g ).
(2) The input gate is basically a hidden layer of sigmoid activation nodes with weighted inputs x t , h t-1 . So, the expression for the input gate is given by:

i = σ(b i + x t U i + h t-1 V i ). ( 3
)
Denoting by s t the internal state of a LSTM cell, the forget gate is a set of sigmoid activation nodes multiplied by s t-1

to determine what must be remembered (output of forget gate close to 1) and what must be forgotten (output of forget gate close to 0). So, LSTM can learn according to the context and the forget gate expression is given by:

f = σ(b f + x t U f + h t-1 V f ).
(4) Then, the output of the forget gate acts as a weight for the internal states. So, the output of this state s t is expressed by:

s t = s t-1 • f + g • i.
(5) The output gate is the final stage of the LSTM cell. The output gate has two components: a sigmoid layer and a hyperbolic tangent (tanh) layer (which creates a vector with new candidates to be added to a ct state). The output is expressed by:

o = σ(b o + x t U o + h t-1 V o ).
(6) So, the final expression for cell output is: Other works deal with the problem of prediction using other approaches. Predictive Recurrent Neural Networks (PredRNN) are introduced in [START_REF] Wang | Predrnn: recurrent neural networks for predictive learning using spatiotemporal lstms[END_REF]. Pre-dRNN is based on the idea that prediction strategies should store spatial and temporal data on a unique memory pool. This idea demands a new network architecture called Spatiotemporal LSTM (ST-LSTM) that extracts and memorizes spatiotemporal representations simultaneously. The authors in [START_REF] Kwon | Predicting future frames using retrospective cycle gan[END_REF] use Generative Adversarial Networks for predicting video frames.

h t = tanh(s t ) • o (

Convolutional LSTM (ConvLSTM)

Convolutional LSTM is an adequate approach to spatiotemporal prediction: while LSTM is used for a sequence of collected data along time that is modeled as a temporal series, Convolutional Neural Networks are widely used for image processing. Thus, for a sequence of images, a Convolutional LSTM (ConvLSTM) combines the features and advantages of LSTM and CNN, which makes ConvLSTM the best approach for image prediction.

Cell outputs, hidden states, and gates of the ConvLSTM are represented by 3D tensors whose last two dimensions are spatial dimensions (rows and columns). ConvLSTM calculates the future state of a cell by the inputs and past states of its local neighbors using convolution in the stateto-state and input-to-state transitions. The convolution is preceded by a padding operation to ensure the dimension compatibility between states and inputs as is stated in [START_REF] Shi | Convolutional lstm network: A machine learning approach for precipitation nowcasting[END_REF].

The first step in frame prediction is the visualization and analysis of files, extracting rainfall image data. An input example for estimators is shown in Figure 2

METHODOLOGY

Two approaches were used and compared, with different settings applied to a dataset composed of 14 colored and 10 grayscale images of rainy days.

Application of ConvLSTM Models, First Approach

In the present work, two different scenarios were tested, during February 2019 and February 2020 for the Fortaleza region with slightly different Convolutional LSTM architecture configurations for frame prediction given an image sequence.

On the first approach, the selected coordinates were 7°00'00.0"S 40°00'00.0"W, 2°00'00.0"S 37°00'00.0"W, and the samples were obtained between February 1st, and 2nd 2019.

Three 2D convolutional LSTM networks using tanh activation function and a convolution window of size 7 × 7 were used. The number of output filters in the convolution for the first network was 128 and for the remaining was 64. A last 3D convolutional LSTM network using only one output filter and with a convolution window with dimension 1×1×1 and using ReLU activation function was used, as precipitation level assumes only values equal or superior to zero, so, the negative activations were mapped to zero.

Adam algorithm was used as optimization function. Adam optimization is a stochastic gradient descent technique used for large data and parameter problems based on adaptive estimation of first and second-order moments. It is computationally efficient, demanding low memory. It is invariant to diagonal re-scaling of gradients according to [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF].

In neural networks, we aim to minimize error, so, in this scenario, the objective function (in the optimization sense) is referred as a cost function, and the result of this cost function is called loss. Let x denote the error (y predictedy real ), the loss function is defined as the logarithm of the hyperbolic cosine of the prediction error:

log(cosh(x)) = log e x + e -x 2 (8)
The algorithm runs during 1000 epochs. The code was developed using Python, Keras and Tensorflow libraries. It initially defines latitude and longitude limits and separates the dataset into training and test datasets and then starts getting the current day as a copy (datetime object), gets the oldest day, and set as current day. A slice (latitude / longitude) data is captured for the current day, getting the precipitation matrix. For each pair of matrices, two periods (k as t and k + 1 as t + 30min) are registered as X and Y . The batch size was set to 2. Then the training period starts and a model is generated and stored. The training history is transformed into a dataframe and saved in a ".csv" format. Then the algorithm advances to the next training day, setting the training dataset to its original state.

Application of ConvLSTM Models, Second Approach

For the second method, reflectivity data from the X-band radar in Fortaleza, Brazil, obtained from FUNCEME was used to create the dataset.

The selected coordinates were 5°00'00.0"S 40°00'00.0"W, 3°00'00.0"S 37°00'00.0"W, and 184 sample data points separated by a time-step of 7 to 8 minutes were obtained every day for 11 days between February 1st, 2020 and February 11th, 2020.

For pre-processing, we first retrieved radar images of size 360 × 360 from the raw data file and transformed the intensity values into gray-level pixels. To address computational limitations and the absence of data in the outer edges of the images, we decided to crop and cut the images. For each image, we retrieved the inner square of dimensions 256 × 256, which we then cut in four squares of equal dimension.

Patches from the same position were used to create sequences of 10 images. Such a sequence is shown in figure 4. Each sequence was then added to the final dataset, which was saved as a N-dimensional array of shape (808,10,128,128,1).

Figure 3. Sequence of frames obtained from meteorological radar (Fortaleza) images from FUNCEME (accessed on June 15, 2021).

The model was constructed using three 2D convolutional LSTM (ConvLSTM2D) layers with batch normalization, followed by a 3D Convolutional (Conv3D) layer for the spatiotemporal outputs. For the three ConvLSTM2D layers, the number of output filters was 64 and the ReLU activation function was used. The kernel sizes were the following: 5 × 5 for the first layer, 3 × 3 for the second layer, and 1 × 1 for the third layer. The Conv3D layer used one output filter, a kernel size of 3 × 3 × 3, and the sigmoid activation function.

Similarly to the first model, the Adam optimization function was used. The loss function used was binary cross entropy. A cross-entropy loss increases the performance (in terms of learning rate and generalization) of a model with sigmoid and softmax output when compared with the mean squared error loss, which presents problems such slow learning rate and saturation.

Data from the first 8 days were used for training and validation. The trained model was then evaluated on testing data for the following 3 days. The algorithm was run for 15 epochs, using batch sizes of 4 instead of 2 to reduce the computational cost.

RESULTS AND DISCUSSION

The results for both approaches are shown below, including discussions of qualitative and quantitative aspects and comparisons with the performance of some works presented in Section 2.

Results for the First ConvLSTM Approach

Figure 4 shows the comparison between real and predicted frames at the same moment for different time instants (30, 60, 90, 120 minutes).

The losses were under 0.02. As it is possible to see, ConvLSTM can predict with good accuracy the dynamic behavior of the frame sequence.

Results for the Second ConvLSTM Approach

The trained model was able to predict images up to 30 min ahead based on 50 minutes of previous data. Figure 6 shows the comparison between real and predicted frames at the same moment for different time instants (8, 15, 23, 30 minutes).

The predicted images all had grey backgrounds; for this reason, we decided to apply filters to each image, thus matching the observed images and allowing for more accurate comparisons. Figure 7 illustrates such changes.

The losses were the following: training: 0.0538, validation: 0.0577, testing: 0.0542.

This ConvLSTM model was able to get relatively accurate results with little data and training. We found these results to be transposable, as the model did not need to be trained on the same days used to make the predictions.

Discussion

The tables 1, 2 and 3 summarize the architecture and parameters for the 1st and 2nd approach and ConvCast according to discussed in [START_REF] Koutník | Convcast: An embedded convolutional lstm based architecture for precipitation nowcasting using satellite data[END_REF]. As is possible to see, the 1st architecture has the smallest loss, followed by the 2nd approach and ConvCast. It is possible to compare this performance indicator with the present one shown in [START_REF] Shi | Convolutional lstm network: A machine learning approach for precipitation nowcasting[END_REF], where it is used a Mean Square Error (MSE), resulting in 1.4. However, the 1st approach demands a significantly larger number of epochs (1,000 against 15 of 2nd architecture and ConvCast). Is important to note that the prediction of 2nd approach was based on less than one hour of observations. In addition to a reduced loss, the number of inputs necessary is expressively smaller: from 10 up to 14 against the values presented on the table 3 of ConvCast.

CONCLUSION AND FUTURE WORKS

The ConvLSTM algorithms were shown to provide feasible results with a reduced number of inputs for training and testing, with predicted and actual frames presenting similar profiles. In one of the presented configurations, the number of necessary epochs is equal to the presented in [START_REF] Koutník | Convcast: An embedded convolutional lstm based architecture for precipitation nowcasting using satellite data[END_REF]. The method can be used, for example, for hurricane and typhoon monitoring, allowing early warning on regions under the constant presence of these phenomena.

A possible subject to explore is the use of Adversarially Learned Inference (ALI) and Generative Adversarial Networks (GAN). A second alternative would be the use of Physics Informed Neural Networks, discussed in [START_REF] Raissi | Physics-informed neural networks: a deep learning framewoek for solving forward and inverse problems involving nonlinear partial differential equations[END_REF]. As a next step, a second line of study is the use of a stochastic process based-approach to find a generalized model (probably an α-stable like process) -more precisely, a Probability Density Function (PDF) for climatological phenomena that converges to fractal representation on small scale, Gaussian Random Field on a medium scale and Generalized Neyman-Scott Process and Coxer process to model large scale phenomena such as frontogenesis. This PDF can be used to extract important statistical features that can be used for weather prediction as well.

Another possibility is the use of stochastic resonance to improve the quality of radar/satellite images.

  7)A piece of information is forgotten only when a new entry takes its place.

Figure 1 .

 1 Figure 1. LSTM architecture.

  , containing meteorological radar images obtained from Cearense Foundation of Meteorology and Hydric Resources (Fundação Cearense de Meteorologia e Recursos Hídricos -FUNCEME), website accessed on April 23rd, 2021.

Figure 2 .

 2 Figure 2. Sequence of frames obtained from meteorological radar (Quixeramobim) images from FUNCEME (accessed on April 23rd, 2021).

Figure 4 .

 4 Figure 4. Frame prediction of the first ConvLSTM approach: comparison between real and predicted frames. The observed (left) and predicted (right) frames are compared side-by-side for 30 and 60 minutes ahead (first row) and 90 and 120 minutes ahead (second row).

Figure 5 .

 5 Figure 5. Frame prediction of the second ConvLSTM approach: comparison between real and predicted frames. The observed frames are shown on the first row and the respective predicted ones are shown on the second row for 8, 15, 23 and 30 minutes ahead.

Figure 6 .

 6 Figure 6. Frame prediction of the second ConvLSTM approach with filters applied: comparison between real and predicted frames.

Table 1 .

 1 Summary of architecture and parameters for the 1st approach

		1st Approach
		3 ConvLSTM layers
		Activation function: tanh
		Convolutional window: 7 × 7
		28, 64, 64 output filters respectively
	Structure
		1 3DConvLST:
		Activation function: ReLU
		Convolutional window: 1 × 1 × 1
		1 output filter
	Batch size	2
	Optimization function	Adam
	Loss function	logcosh(x)
	Input	14 coloured images
	Number of epochs	1,000
	Loss	0.02

Table 2 .

 2 Summary of architecture and parameters for the 2nd approach

Table 3 .

 3 Summary of ConvCast architecture and parameters according to Kumar et al. (2020).

		ConvCast
		Activation function:
		tanh (ConvLST)
		ReLU (3DConvLSTM)
	Structure	4 hidden layers
		128 input-to-state filters
		64 state-to-state
		kernel size: 7 × 7
	Batch size	2
	Optimization function	Adam
	Loss function	MSE
		1,276 in the training set
	Input	319 in the validation set
		242 for the test set.
	Number of epochs	15
	Loss	0,6049

ACKNOWLEDGEMENT

The authors would like to thank Institut de Recherche pour le Développement for the financial support.