Aymeric Mainvis

Vincent Fabbro

Christophe Bourlier

Henri-Jose Mametsa

Spectral Decomposition Method for Large Sea Surface Generation and Radar Backscattering

Keywords: Sea surface, spectral method, small slope approximation, radar scattering, ocean remote sensing

This paper presents an original method to simulate rapidly a realization of a large sea surface and also investigates its impact on the backscattering normalized radar cross section (BNRCS). Indeed, for radar microwaves frequencies an for wind speed between 1 and 10 m/s the area of the sea surfaces must be ranged from 50 to 200 m 2 to take into account all the surface roughness scales, which can contribute to the scattering process. The proposed method consists in splitting up the full sea surface height spectrum into sub-spectra of smaller scales and recombining the resulting sea height generations from different interpolations and recombination techniques. A closed-form expression of the resulting sea surface height spectrum is also derived to interpret the simulation results. Finally, from the firstorder Small Slope Approximation (SSA1) scattering model, the BNRCS computed from sea surfaces generated by the suggested method and the conventional method are compared to show the efficiency of the method in terms of accuracy and memory consumption.

I. INTRODUCTION

Deterministic model -and sea surface generation in particular-is a key-element within the radar sensor design optimization but also within the modeling and understanding chain of measurements. Either from airborne or satellite sensors or from satellite ones, sea surface generation makes it possible to compute the radar backscattering and therefore, to investigate the recorded data, like altimeters data. However, some radar systems (like satellite Synthetic Aperture Radar) require a sea surface generation with large area and high resolution. So, it becomes crucial to offer low-consumption sea surface generation model in terms of computing resources and / or time requirement. Indeed, electromagnetic (EM) wave scattering from random rough surfaces models -like the first-order Small Slope Approximation (SSA1) method [START_REF] Voronovich | Small-slope approximation in wave scattering by rough surfaces[END_REF]-require a small surface sampling interval to obtain accurate results. Commonly, for convergence criteria, this sampling interval is equal to one-tenth of the incident radar wavelength. Furthermore, a large scale of waves within the sea surface have to be taken correctly into account to simulate the sea surface geometry. In particular, the gravity and gravity-capillary waves drive the EM scattering from the illuminated sea surface and so constrain the sea surface area. Therefore, EM scattering of large sea surfaces remains a critical computing constraint.

August 23, 2017 DRAFT

A fast simulation of a sea surface is described in [START_REF] Pinel | Fast simulation of a moving sea surface remotely sensed by radar[END_REF]. Pinel et al. studied the slope probability density function and the slope autocorrelation function after dividing the height sea spectrum in two parts and generating sea surfaces with different spatial resolutions. In [START_REF] Jiang | Spectral Decomposition Modeling Method and Its Application to EM Scattering Calculation of Large Rough Surface With SSA Method[END_REF], Jiang et al. introduced the spectral decomposition method to reduce memory consumption and generate different-scale rough surfaces. The whole spectrum is divided in several parts used to generate these different-scale rough surfaces. This method is applied to perform unified device architecture (CUDA) parallel computation. This approach and the conventional one are tested and compared with the first-order SSA1.

The same method is then studied for sea surface generation in [START_REF] Jiang | EM scattering calculation of large sea surface with SSA method at S , X , Ku , and K bands[END_REF] and tested with SSA1 by simulating the sea surface backscattering normalized radar cross section (BNRCS) and Doppler spectra. This paper provides a quantitative analysis of the spectral decomposition method. Indeed, this particular sea surface generation is analytically described and developed to express its computational complexity and its memory consumption. Furthermore, an in depth study is also performed to highlight the impact of the interpolation process and the two suggested combination techniques on the sea surface characteristics and on the BNRCS.

Section II details the formalism of the spectral decomposition method which describes a split-spectrum process and a reconstructed sea surface generation with an interpolated surface and a combination technique. Section III introduces the SSA1 method and the sea surface BNRCS expression. Section IV presents numerical results for a two-dimensional problem by evaluating the sea surface height spectrum and the height structure function. BNRCS computed from SSA method considering a conventional sea surface generation and the Spectral Decomposition Method (SDM) are described in Section V. The influence of SDM parameters are discussed.

II. SEA SURFACE GENERATION AND SPECTRAL DECOMPOSITION METHOD

A. Sea Surface Model

From a conventional spectral method, the height of sea surface H(r, t) is expressed as [START_REF] Tsang | Scattering of Electromagnetic Waves: Numerical Simulations[END_REF] H

(r, t) = Re R 2 S(k)E(k)e -jω(k)t e jk•r dk , (1)
where r = r(x, y) the Cartesian position, t the time, S(k) the height sea spectrum, k the wavenumber vector, E a gaussian process and ω(k) the pulsation defined through a dispersion relation [START_REF] Elfouhaily | A unified directional spectrum for long and short wind-driven waves[END_REF]. This conventional method is very efficient for sea surface generation, allowing to use Fast Fourier Transform (FFT). However, EM scattering computation using rigorous techniques requires an accurate sampling of the surface, inducing prohibitive computing ressources for high frequency and for high sea state in a three-dimensional problem. For this reason an optimization of the method is proposed by applying a spectral decomposition.

B. Spectral Decomposition Method

To introduce the Spectral Decomposition Method; first, the function Γ is defined as

Γ(k, t) = S(k)E(k)e -jω(k)t . (2)
Then, this function is split up into N functions Γ n such as

Γ n (k, t) =    Γ(k, t) if k n ≤ k < k n+1 0 else, (3)
August 23, 2017 DRAFT with Γ defined in Eq. ( 2), k the wavenumber vector, k n the cutoff-wavenumber, for which k 0 = 0, k N = +∞ and n ∈ [0, N -1]. Consequently, one has to choose N -1 cutoff-wavenumbers k n to define Γ n . Eq. ( 1) can be rewritten as

H(r, t) = Re N -1 n=0 k <kn+1 k =kn Γ(k, t)e jk•r dk (4) = Re N -1 n=0 R 2 Γ n (k, t)e jk•r dk = Re N -1 n=0 h n (r, t) ,
with h n (r, t) the height of sea surface generated from the n-th part Γ n . Here, the full generated sea surface H(r, t)

is seen as the summation of N elementary generated sea surfaces, corresponding to a specific roughness scale.

C. Reconstructed Sea Surface

Computing sea surface implies choosing a surface size L x × L y (or N x × N y sampling points) and sampling intervals (∆ x , ∆ y ). For more clarity, the surface length and the sampling interval to generate the sea surface H are chosen such as L x = L y = L 0 and ∆ x = ∆ y = ∆, respectively. Then, by applying SDM and from Eq. ( 3), the use of FFT to generate each surface -corresponding to the different roughness scales-implies

k n = 2π L n k 0 = 2π L 0 k N = π ∆ , (5)
with n ∈ [0, N -1], L n = N n ×∆ n the length of the sea surface h n linked to the cutoff-wavenumber k n in Eq. ( 4)

with N n sampling points and ∆ n the sampling interval. So, selecting cutoff-wavenumber k n in SDM leads to the parameters of h n

L n = 2π k n ∆ n = 2π N n k n . (6)
Eq. ( 3) and (6) imply

L n > L n+1 , (7)
and Eq. ( 6) leads to

∆ n > ∆ n+1 , (8)
if N n = N , a constant. Consequently, generating sea surfaces by using SDM induces a sea surfaces serie, each sea surface with different sizes and different sampling intervals. However, to be able to sum the different surfaces corresponding to the different roughness scales, the surface meshings must be equal. Then, to solve this problem, two techniques are suggested: an interpolation process and a combination technique. Fig. 1 plots a schematic diagram for the generation of surfaces h n and h n+1 and their respective length, L n and L n+1 , and sampling interval, ∆ n and ∆ n+1 according to the SDM. The interpolation process is a way to reduce the sampling interval from ∆ n to ∆, the smallest sampling interval. In this paper, three kinds of interpolation are studied, namely, linear, quadratic and cubic. The combination technique is a way to increase the surface length from L n to L 0 , the largest surface length. For example, the interpolation process is a mean to get the same size between the two sea surfaces in August 23, 2017 DRAFT Fig. 1 and the combination technique, the mean to get the same sampling interval. Thus, these two surfaces can be finally additionated. For a sake of clarity, a focus on a two-dimensional problem and a spectrum division into two parts is proposed to formalize the combination techniques. Thus, the total sea surface H is composed of a low-frequencies-scale sea surface h LF and a high-frequencies-scale one h HF

H(x) = h LF (x) + h HF (x), (9)
h LF is the interpolated sea surface and h HF the combined one. Two combination techniques are studied: the Repeated Surfaces Technique (RST) and the Combined Surfaces Technique (CST). Then, the RST principle is that the final HF surface is composed of P times the same HF surface. It can be formalized by

h HF,RST (x) = h HF (x) * P -1 n=0 δ(x -nL), (10)
with * the product-convolution operator, h HF the P -times-repeated surface, L its length, h HF,RST the composed surface of length P ×L and δ the Dirac delta function. This combination technique ensures the continuity between each combined surface h HF thanks to the FFT property. The CST is defined by Jeannin et al. [START_REF] Jeannin | A Large-Scale Space-Time Stochastic Simulation Tool of Rain Attenuation for the Design and Optimization of Adaptive Satellite Communication Systems Operating between 10 and 50 GHz[END_REF], a composite surface h comp is defined by

h comp (x) = √ d -xh 1 (x + L -d) + √ xh 2 (x) √ d , (11)
with x ∈ [0; d], h 1 and h 2 two independent rough surfaces with length L. These two surfaces are combined on an interval d. Then,

h HF,CST (x) = P -1 n=0 h HF,int,n (x) * δ[x -n(L -d)], (12)
with

h HF,int,n (x) =    h HF,comp,n-1 (x) if x ∈ [0; d] h HF,n (x) if x ∈]d; L -d], (13)
and h HF,n the rough surface generated from a unique random gaussian process and Thus h HF,int,n is a rough surface of length (L -d). Furthermore,

h HF,comp,n (x) = √ d -xh HF,n (x + L -d) + √ xh HF,n+1 (x) √ d . (14)
h HF,comp,-1 (x) = √ d -xh P -1 (x + L -d) + √ xh 0 (x) √ d , (15)
to ensure the combined sea surface periodicity and so, the continuity. The composed surface h HF,CST length is equal to (L -d)P . This technique keeps continuity between each combined surface h HF,int,n . For a simplicity reason, the interval d is chosen as L/2 in this work. Fig. 2 illustrates a schematic diagram for the generation of the surface h HF with any of the two combination techniques, RST (2a) and CST (2b). To check the sea spectrum integrity, the HF sea height spectrum from RST is expressed. It can be written as

S HF,RST (k) = S HF (k) P sin kP L 2 sin kL 2 2 (16)
with S HF,RST the spectrum from RST, S HF the sea height spectrum used to generate the P combined surfaces of length L and k the wavenumber. The proof is detailed in Appendix A. Then, from Eq. ( 16), it appears that S HF,RST is the conventional sea spectrum S HF modulated by a 2π/L-periodic function. This function has local maxima for

kL 2 = nπ ⇔ k = n2π L , (17)
with n ∈ Z. To quantify the efficiency of the SDM, the computational complexity of the FFT is a relevant tool. This is expressed as

O(N log 2 N ), (18)
August 23, 2017 DRAFT

with N the number of samples computed by FFT. Let us consider a simple 3D case, as previously discussed, the spectrum is divided into two parts like in Eq. ( 9), that is

H(r, t) = h LF (r, t) + h HF (r, t), (19)
with h LF the interpolated sea surface and h HF the reconstructed one. According to the chosen combination technique, the computational complexity C HF linked to the surface generation h HF is

C HF =    O(N 2 HF log 2 N 2 HF ) if RST L L-d P × O(N 2 HF log 2 N 2 HF ) if CST, (20)
with N 2 HF the number of samples of each combined surface of area L 2 , d the CST parameter in Eq. ( 12) and P the number of combined high-frequencies-scale surfaces. So, the computational complexity C T to generate the sea surface H is

C T = O(N 2 LF log 2 N 2 LF ) + C HF , (21)
with N 2 LF the number of samples of h LF before interpolation process. For example, suppose N LF = N HF = N , then,

C T = (1 + α) × O(N 2 log 2 N 2 ), (22)
α = 1 or P L/(L -d) Eq. ( 20). The equivalent computational complexity C ref for a conventional sea surface generation is

C ref = O(P 2 N 2 log 2 P 2 N 2 ). (23)
Indeed, with a given number of samples N 2 and a given sampling interval ∆, the total area of the generated sea surface thanks to SDM is L 2 = (P ×N ×∆) 2 . So, by keeping the same sampling interval, (N ×P ) 2 sampling points are needed to reach the same area. Fig. 3 sets out the computational complexity of sea surface generation versus the number of samples N with P = 8 according to Eq. ( 22) (RST and CST) and Eq. (23) (Reference). This result shows a gain between 5 and 14 (for a number of samples N = 10 4 ) by using SDM rather than a conventional sea surface generation. Fig. 4 plots the computational complexity of sea surface generation versus the number of combined surfaces P with N = 2 13 . This time, the gain is between 10 (for CST) and 200 (for RST) by using SDM and considering a parameter P = 16. These simulations clearly highlight the interest of such a multiscale method.

For memory consumption, by keeping the same notations introduced in Eq. ( 22) and Eq. ( 23), the total memory consumption to store generated sea surface data is

M ref = mP 2 N 2 (24)
M T = m(1 + α)N 2 , (25)
with m the memory allocation of an elementary piece of data, M ref the memory consumption for a conventional sea surface generation and M T the one for SDM with α = 1 for RST or α = P L/(L -d) for CST Eq. (20).

According to Elfouhaily et al. [START_REF] Elfouhaily | A unified directional spectrum for long and short wind-driven waves[END_REF], the minimum surface wavenumber k min should verify k min ≈ 0.3k p with

k p ≈ Ω 2 g/u 2 10 , (26)
August 23, 2017 DRAFT where Ω is the inverse wave age equal to 0.84 in a case of a fully developed sea, g the acceleration of gravity and

III. SIMULATED RADAR BACKSCATTERING

To analyze the impact of SDM on EM scattering, the BNRCS is considered, computed by a local model, the first-order Small Slope Approximation (SSA1).

First-Order Small Slope Approximation

The SSA1 is efficient in the whole range of incidence angles, from 0 • (nadir) to 60 • . The scattering operator is

given by [1]

S(k, k 0 ) = B(k, k 0 ) Q z r e -jQzη(r) e -jQ H •r dr, (27)
with B(k, k 0 ) the first-order small perturbation model (SPM1) kernel [START_REF] Pinel | Electromagnetic Wave Scattering from Random Rough Surfaces: Asymptotic Models[END_REF], a polarization term, Q H and Q z horizontal and vertical components of the vector Q = k-k 0 , respectively, k 0 and k the incidence and observation wave vectors and η(r) the surface elevation. Then, the incoherent normalized radar cross section (NRCS) of an infinite surface σ 0 is expressed as

σ 0 (k, k 0 ) = S(k, k 0 )S * (k, k 0 ) κΣ - S(k, k 0 ) S(k, k 0 ) * κΣ , (28)
with S(k, k 0 ) defined in Eq.( 27), Σ the sea surface area (length in a 2D problem) and κ a constant equal to π for a 3D problem and 4k 0 for a 2D problem with k 0 the radar wavenumber. In this numerical approach, a Thorsos beam [START_REF] Bourlier | Method of Moments for 2D Scattering Problems: Basic Concepts and Applications[END_REF] of parameter g = L/3 (with L the total length of the sea surface) is considered to illuminate the generated sea surface. This beam is a tapered plane wave with a Gaussian shape; the tapering is used to reduce the incident field to near zero at the ends of the generated sea surface and so, to reduce the potential edge effects to a negligible levels. From Eq. ( 28) and for a Gaussian process, an analytical expression of the incoherent NRCS [START_REF] Bourlier | Comparison of asymptotic backscattering models (SSA, WCA, and LCA) from one-dimensional Gaussian ocean-like surfaces[END_REF] can also be performed,

σ 0 (k, k 0 ) = |B(k, k 0 )| 2 κQ 2 z e -Q 2 z σ 2 η r e -jQ H •r e Q 2 z W (r) -1 dr (29) = |B(k, k 0 )| 2 κQ 2 z r e -jQ H •r e -1 2 Q 2 z D(r) -e -Q 2 z σ 2 η dr,
with σ 2 η the height mean square value, W the height autocorrelation function and D the height structure function defined as

D(r) = 2 σ 2 η -W (r) . (30)
August 23, 2017 DRAFT So, the BNRCS is directly linked to the Fourier transform of a function F which is related to the sea surface characteristics,

F(r) = e -1 2 Q 2 z D(r) . (31)
Then, by applying the SDM and studying the function F, one can ensure a correct estimation of the BNRCS.

IV. GENERATED SURFACE CHARACTERISTICS

It is necessary to analyze the generated surface characteristics to check its realism in comparison to that obtained from the conventional method. First, the interpolation process impact (for LF sea surface generation) on sea surface height spectrum is investigated. Secondly, the generated surface characteristics from combination techniques (for HF sea surface generation) introduced in subsection II-C are studied. For a sake of clarity, this study is conducted for 2D problems but the results can be extended to 3D problems.

A. Interpolation Techniques

One scenario is proposed here and the parameters are listed in Tab. II. Before the interpolation process, the cutoff-wavenumber k c is π/(P ∆); after, it becomes equal to π/∆ with P = 8 the interpolation parameter. Fig. 5 illustrates the isotropic part of sea surface height spectrum from the model of Elfouhaily et al. [START_REF] Elfouhaily | A unified directional spectrum for long and short wind-driven waves[END_REF] versus the wavenumber k. Three interpolation techniques are studied: linear, quadratic and cubic. The full sea surface height spectrum is obtained by both a numerical computation (S LF (k) Num) -which consists in a Monte-Carlo process by generating 500 sea surfaces and then by computing their sea surface height spectrum-and the analytical expression [START_REF] Elfouhaily | A unified directional spectrum for long and short wind-driven waves[END_REF] (S LF (k) Theory). First, one can observe that interpolation processes create higher frequency harmonics than the original generation. Also, it can be seen that the quadratic interpolation presents over-occurred harmonics which can severely disturb the BNRCS, especially by using the Small Perturbation Method (SPM), that is directly proportional to high-frequencies sea surface height spectrum. Besides, linear and cubic interpolations seem to be relevant techniques to upgrade the sampling intervals of a given sea surface, creating low energy high frequency components. So, linear interpolation is the best choice by adding time and/or computer resources optimizations. Wind speed is u 10 = 8 m/s. Three interpolation parameters are presented, P = {8; 16; 32}, the interpolation method is linear.

discussed after adding the HF surface. Moreover, an interpolation process -especially when linear-is efficient to reduce the sampling interval at almost no added cost. Besides, the more the interpolation parameter P is, the earlier the oscillations occur in the sea surface height spectrum. This phenomenon is explained by the chosen sampling August 23, 2017 DRAFT interval to reach. Indeed, before the interpolation process, the maximum wavenumber value is k max = π/(P ∆), so, the more the interpolation parameter P is, the shorter k max is and thus, the earlier the oscillations occur.

B. Combination Techniques

The scenario is similar to that one in subsection IV-A, Tab. II. While, P surfaces are now combined thanks to one of those techniques presented in subsection II-C. Before the combination process, the HF surface length L is N ×∆ and after, it will be P ×L with P the number of combined surfaces. Thus, the cutoff-wavenumber before the combination technique is k min = 2π/L. Fig. 7 plots the isotropic part of sea surface height spectrum versus the surface wavenumber. This spectrum is compared to those performed by using a combination technique. Fig. 7a illustrates the RST, the theoretical spectrum of RST previously derived in Eq. ( 16) is also displayed and is in accordance with the numerical one. The RST slightly overestimates the harmonics within the spectrum. Despite the apparition of harmonics lower than that located in k min , the CST seems to get the best accuracy by ensuring continuity and preserving from overestimated harmonics, Fig. 7b. Moreover, the SDM height mean square value (σ 2 HF, X , X the combination technique) is in accordance with the conventional one (σ 2 HF ). Indeed, σ 2 HF = 0.084 m 2 , σ 2 HF, RST = 0.086 m 2 and σ 2 HF, CST = 0.083 m 2 . Fig. 8 plots the combined sea surfaces height spectrum versus the surface wavenumber k by using the CST. Whatever the number P of combined surfaces is (between 8 and 32), the height spectrum is qualitatively similar.

C. Height Spectrum and Height Structure Function

The SDM is applied to create a N ×P -samples sea surface with a sampling interval ∆. Firstly, one sea surface with N samples and a sampling interval P ×∆ is generated and then linearly interpolated to get a new sampling interval ∆. Secondly, one sea surface with N samples and a sampling interval ∆ is generated to perform RST and CST and therefore, to create a combinated sea surface with N ×P samples and a sampling interval ∆. Then, these two surfaces are added to generate the composite two-scales surface. Notice that, to avoid spectral redundancy between the two spectra used to generate these two surfaces, harmonics in the interval I are forced to 0 in the first spectrum -that is the LF part-with

I = 2π N ∆ , π P ∆ . ( 32
)
The frequency is 10 GHz, N = 2 13 samples, ∆ = λ 0 /10 with λ 0 the wavelength, P = 8 and the wind speed u 10 is 8 m/s. This generation is repeated over a Monte-Carlo process by generating 500 composite two-scales sea surfaces. Fig. 9 plots the isotropic part of sea surface height spectrum versus the surface wavenumber. This spectrum is compared to those performed by using the SDM. Once again, harmonics lower than that located in k c = 2π/(N ∆) are greater than their theoretical counterparts in CST, Fig. 9b. Indeed, this technique is based on the combination of independent surfaces. Fig. 9a illustrates the RST. As previously depicted in Fig. 7, the latter overestimates the harmonics within the spectrum. Finally and again, the CST gets the best accuracy by ensuring continuity and preserving from overestimated harmonics. After a spectral investigation of SDM, a spatial analysis is necessary through the heigh structure function introduced in Eq. (30). Indeed, this quantity leads to the BNRCS August 23, 2017 DRAFT estimation by using SSA1 (Eq. ( 29)). Fig. 10 plots the theoretical height structure function (D(x) Theory) versus the surface position. This height structure function is compared to the two obtained thanks to SDM. The RST produces oscillations within the height structure function. This phenomenon is produced by the repetition process and so, by the correlation renewal between one surface elevation point and its copy, located every N ×∆ meters.

The CST height structure function is qualitatively in accordance with the theoretical one.

August 23, 2017 DRAFT

D. From sea surface characteristics to BNRCS

The right description of the function F defined in Eq. ( 31) is a crucial step into the BNRCS computation. Indeed, as previously described, this function is one of the key-parameter in the analytical expression of the BNRCS with SSA1 Eq. (29). Therefore, by ensuring a non-impact of combination techniques on this function, SDM becomes a relevant way to compute the BNRCS from sea surfaces. Fig. 11 plots the theoretical F function versus the surface position. Those computed by using SDM and the two different combination techniques, RST and CST, are also displayed. Two incidence angles are considered here, θ = 0 • , which is located in the Geometrical Optics domain, and θ = 60 • in the Bragg scattering domain. Whatever the combination technique is, the SDM is in accordance with the theory and so, will not disturb the BNRCS estimation.

V. SEA SURFACE BNRCS

A two-dimensional problem is considered to compute the sea surface BNRCS. The same parameters introduced in Tab. II are chosen and the SDM parameter P (that is the number of combined surfaces) is 8. The sea surface BNRCS is computed with a monostatic configuration and HH polarization. The sea dielectric permittivity ε is 53.2 + j37.8.

To obtain this BNRCS, a thousand of sea surfaces are generated. Thus, the surface length is N×P ×∆ ≈ 196 m and the gravity waves are correctly taken into account in the sea surface height spectrum, Eq. (26). For this scenario, we study the impact induced by the combination technique -and therefore the SDM-on the sea surface BNRCS. of the suggested combination techniques do not quantitatively disturb the BNRCS estimation, both in VV and HH polarizations. For RST, the maximal error is ±1 dB and for CST, the error is about 0 dB after the incidence angle 15 • and remains inferior to ±1 dB along the incidence angle track. The error function is similar whatever the polarization is. Indeed, only the sea surface generation process is modified and this does not impact on the polarization term within SSA1, Eq. ( 27). Thus, SDM with RST is an efficient mean to perform numerical computation of the sea August 23, 2017 DRAFT to an offset (-10 dB for CST). The results show a same trend for VV or HH polarization, the tested parameter P values have no impact on the result. This observation is confirmed by Fig. 13c and13d. Indeed, the BNRCS ratio between P = {16, 32} and P = 8 is inferior to ±1 dB along the incidence angle track, and so, whatever the investigated combination technique is. Again, the two combination techniques are discriminated thanks to an offset (-2 dB for CST). Like in Fig. 12c and12d, these error functions are similar; the sea surface generation process does not interfere with the polarization term in SSA1, Eq. ( 27).

VI. SUMMARY AND OUTLOOKS

Sea surface generation is a high-demanding process to achieve accurate BNRCS. Indeed, large sea surface areas and high resolution are required. In this context, the Spectral Decomposition Method (SDM) is a relevant tool to make the sea surface generation faster. The former is formalized and the new sea surface generation process is introduced. An interpolation and a combination technique are necessary to complete this generation. Three kinds of interpolation are suggested, namely, linear, quadratic and cubic. Two combination techniques are also investigated, the Repeated Surfaces Technique (RST) and the Combined Surfaces Technique (CST). Then, the SDM computational complexity is compared to the reference one, it appears that the SDM reduces the complexity by a factor between 10

August 23, 2017 DRAFT and 200, depending on the chosen combination technique. Similarly, the memory consumption is drastically reduced by using SDM rather than a conventional method, the ratio is roughly from 10 -2 to 10 -3 . The SDM and these interpolation and combination techniques are studied through the generated sea surface characteristics and through the sea surface BNRCS. About the interpolation process, the linear method is the best choice. Indeed, the linear method is the quickest interpolation process and presents weak distortions within the sea surface height spectrum (crucial parameter since its inverse Fourier transform drives sea surface NRCS from SSA1). Besides, despite different behaviours linked to each of the proposed combination techniques on the sea surface height spectrum and on the height structure function, the sea surface BNRCS computed from the SDM with either the RST or the CST is in agreement with the one from a conventional sea surface generation. Therefore, the SDM is valid from near nadir sensors like altimeters to moderate observation angles. This approach is analytically formalized -both in spatial and frequency domains for RST-and tested for a subdivision in two spectra according to the sea surface characteristics and the BNRCS. This process is a relevant tool to accelerate the radar backscattering computation from large sea surfaces. In future work, it should be coupled with a two-scale electromagnetic model to further speed up the simulation. Moreover, the spectral decomposition method could be used to simulate sea surfaces with range variations of characteristics (wind speed in particular) and then compute composite sea surfaces, closer to the real weather conditions.

August 23, 2017 DRAFT

with S HF (k) the theoretical sea height spectrum. That is the response of a uniform linear array of phased array antenna with S HF the elementary antenna.

Fig. 1 .

 1 Fig. 1. Schematic diagram for the generation of surfaces hn and h n+1 according to the Spectral Decomposition Method.

Fig. 2 .

 2 Fig. 2. Schematic diagram for the generation of surface h HF with any of the two combination techniques, RST (2a) and CST (2b).

Fig. 3 . 8 Fig. 4 .

 384 Fig. 3. Computational complexity of sea surface generation versus the number of samples N with P = 8

u

  10 the wind speed at ten meters above the sea. So, with a sampling interval of one-tenth of the incident radar wavelength and u 10 = 8 m/s; 67, 108, 864 samples are needed to generate a conventional 3D sea surface. That is 536, 870, 912 octets for a float64 (m = 8 octets). For SDM, with α = 1 for RST and P = combined surfaces, M T = 16, 777, 216 octets. The memory consumption ratio is 1/32. Tab. I gives the memory consumption ratio August 23, 2017 DRAFT

Fig. 6 Fig. 5 .

 65 Fig.6plots the isotropic part of interpolated sea surface height spectrum versus the wavenumber k. The wind speed u 10 is 8 m/s. The linear interpolation method is considered. Three values of the interpolation parameter are studied: 8, 16 and 32. This result shows a qualitatively-low impact of the interpolation parameter and has to be

Fig. 6 .

 6 Fig. 6. Isotropic part of the interpolated sea surface height spectrum versus the surface wavenumber k from the model of Elfouhaily et al. [6].

Fig. 7 .

 7 Fig. 7. Isotropic part of the sea surface height spectrum S HF versus the surface wavenumber k from the model of Elfouhaily et al. [6]. Wind speed is u 10 = 8 m/s. The cutoff-wavenumber before the combination process k min is also displayed. The isotropic part of the sea surface height spectrum from the two combination techniques introduced in subsection II-C are also illustrated; RST (7a) and CST (7b).

Fig. 8 .

 8 Fig. 8. Isotropic part of the combined sea surfaces height spectrum versus the surface wavenumber k from the model of Elfouhaily et al. [6]. Wind speed is u 10 = 8 m/s. The inspected combination technique is the CST. Three parameters are shown; 8, 16 and 32.

Fig. 12a andFig. 9 .

 9 Fig.12aand 12b plot the incoherent monostatic BNRCS versus the incidence angle from a conventional sea surface generation -spectral method introduced in subsection II-A-and from SDM with either RST or CST. The ratios between RST / CST and the reference are also showed in Fig.12c and 12d. One can see that the SDM and any

Fig. 10 .

 10 Fig. 10. Height structure function D versus the surface position x from the model of Elfouhaily et al. [6]. Wind speed is u 10 = 8 m/s. The theoretical height structure function D is plotted in full black line. The two combination techniques from subsection II-C are illustrated in dashed-red and discontinuous-blue line; RST and CST respectively.

  Fig. 11. F function versus the surface position x from the model of Elfouhaily et al. [6]. Wind speed is u 10 = 8 m/s for a frequency f = 10 GHz. The theoretical F function is plotted. The two combination techniques from subsection II-C are illustrated; RST and CST.

  Fig. 12. The wind speed is 8 m/s for a frequency f = 10 GHz in VV and HH polarizations. Comparison of the SSA1 from conventional sea surface generation and the SSA1 from SDM, considering the two different combination techniques. 1000 surfaces of length L = 196 m were generated.

  Fig. 13. The wind speed is 8 m/s for a frequency f = 10 GHz in VV and HH polarizations. Comparison of the SSA1 with different P parameters and considering the two combination techniques. 1000 surfaces of length L = 196 m were generated. CST -X dB stands for an offset of X dB to improve the discrimination between the two techniques.

hh 1 P P - 1 n=0 P - 1 m=0hhW 2 S

 1112 HF,RST (x) = h HF (x) * HF (x -nL),(33) with h HF the P -times-repeated surface, L its length, h HF,RST the composed surface of length P ×L and δ the Dirac delta function. Then, the height autocorrelation function W HF,RST from RST is expressed asW HF,RST (r) = HF (x 1 -nL)h * HF (x 1 + r -mL) HF (α n )h * HF (α n + r + (n -m)L) HF (r + (n -m)L),(34)with W HF the theoretical height autocorrelation function, x 1 an abscissa and α n = x 1 -nL. Therefore, by taking the Fourier transform of the height autocorrelation function, one can get S HF,RST (k) = S HF (k) P be simplified by using formulae from geometric series to finally obtain S HF,RST (kHF (k),

  T /M ref versus the number of combined surfaces P and the combination technique. Once again, the SDM is more efficient than the conventional sea surface generation.

	TABLE I
	MEMORY CONSUMPTION RATIO
	Number of combined surfaces	RST	CST
	P = 8	0.031 0.266
	P = 16	0.008 0.129
	P = 32	0.002 0.064

M

ACKNOWLEDGMENT

The authors would like to thank Total company for the fundings and especially Veronique Miegebielle and Philippe Lattes for supporting this work.

August 23, 2017 DRAFT


Benoît Michel
email: benoit.michel@insa-lyon.fr

Marc Clausse

Thermochemical Heat Tranformer for Waste Heat Recovery: consideration of dynamic aspects for design

Keywords: Thermochemical Heat Transformer, Waste Heat, Salt screening

Waste Heat Recovery in the industrial sector often requires to upgrade the temperature of the heat fluxes in order to make them useful for the process. Thermochemical Heat Transformer can play a major role as they can provide high temperature lift when compared to alternative technologies. The aim of the presented work is twofold:  to identify suitable working pairs using water as working fluid taking into account three criteria: COP, temperature lift and equilibrium drop (difference of the operating temperature and pressure to their value at equilibrium)  to run a parametric study, thanks to a dynamic 2-salt/2-reactor model, to identify optimal specific heat power as a function of cycle time, working pairs and reactor heat transfer coefficient. For waste heat temperature at 90 °C and a cold source at 30 °C, the CaCl 2 /H 2 O and Ca(NO 3 ) 2 /H 2 O working pairs were identified as promising. The achieve temperature lift is as high as 90 K for a theoretical COP of 0.42. In dynamic operating, a specific heat power of 325 W/kg was found.

INTRODUCTION

As reported by Ademe, waste heat valorization is a promising solution to dramatically increase the energy efficiency of the industrial sector (ADEME, 2015). According to this report, heat to electricity conversion (e.g. ORC) becomes competitive for temperatures of recovered heat above 150 °C whereas a huge amount of waste heat is available at lower temperature. While sorption systems are industrially available for heat to cold conversion, there is today a great challenge for heat to heat conversion, i.e. to upgrade waste heat temperature from 90-100°C to 150°C and above.

Compared to Electrical Heat Pumps (EHP) and Absorption Heat Transformers (AHT), Thermochemical Heat Transformer (THT) is a far less mature technology [START_REF] Li | A target-oriented solid-gas thermochemical sorption heat transformer for integrated energy storage and energy upgrade[END_REF][START_REF] Visscher | Experimental investigation on a novel solid-gas thermochemical sorption heat transformer for energy upgrade with a large temperature lift[END_REF]. However, THT could potentially tackle the issues of EHP and AHT by achieving high temperature lifts with 1-stage heat driven system without any crystallization or corrosion issues [START_REF] Li | A target-oriented solid-gas thermochemical sorption heat transformer for integrated energy storage and energy upgrade[END_REF]. Furthermore, 2-salt technology (Figure 1a), also called resorption heat pump, allows a reduction of the high-pressure lift occurring in AHT when high temperature lifts are sought [START_REF] Yu | Energy upgrading by solid-gas reaction heat transformer: A critical review[END_REF]. In this case, the evaporator and the condenser are replaced by two reactors filled with a so-called low temperature salt (LT salt). Salt/ammonia pairs have been extensively studied for THT applications while there have been some evaluations of salt/methanol, salt/H 2 O and salt/CO 2 pairs with encouraging results for the latter two [START_REF] Yu | Energy upgrading by solid-gas reaction heat transformer: A critical review[END_REF]. In such systems, a main issue is to select the appropriate salt pair for given operating temperatures (T L , T m & T H ), to achieve high performance.

Hence, the first objective of this paper is to perform a salt screening for industrial waste heat upgrade thanks to a static thermodynamic approach. The second objective is to study the dynamic behavior of the THT for pre-selected pairs. The model developed and validated for this purpose is presented as well. Furthermore, it was chosen to focus on water as working fluid.

STATIC THERMODYNAMIC APPROACH FOR SALT SCREENING

In this section, a static thermodynamic approach is used to perform a rapid screening of available salt-H 2 O pairs [START_REF] N'tsoukpoe | A systematic multi-step screening of numerous salt hydrates for low temperature thermochemical energy storage[END_REF][START_REF] Visscher | Experimental investigation on a novel solid-gas thermochemical sorption heat transformer for energy upgrade with a large temperature lift[END_REF]. It is based on the following set of equations. Eq. 1 states the hydration equation of a salt (S) by a vapour (G). The reaction equilibrium is given by Eq. ( 2), where p eq is the equilibrium pressure at T eq as shown in Eq. (3). Eq. ( 1)

( ) ( ) Eq. (2)
Eq. ( 3)

For the sake of clarity, the work presented in this paper is focused on given heat source temperatures representative of operating conditions suitable for heat upgrading in the industry. The cold source temperature is chosen at T L = 30 °C, the waste heat temperature at T m = 90 °C and a high temperature source at T H = 150 °C.

As shown on Fig. 1a, during the regeneration phase (phase 1), waste heat at T m is supplied to the HT reactor and the LT salt is connected to a heat sink at T L . In order to allow simultaneously the hydration of the LT salt and the dehydration of the HT salt, the operating conditions (vapor pressure and temperature) of the LT reactor must be in the blue zone presented on Fig. 1b. Therefore, this reactor has to be connected to a heat sink at T L < T Lmax (i.e. the equilibrium pressure of the LT salt at T L must be inferior to the equilibrium pressure of the HT salt at T m ). In the case of T L is superior to T L,max , the operating conditions do not allow the hydration of the LT salt and the dehydration of the HT salt. At the end of this phase, the reactive LT salt is fully hydrated and the HT salt is fully dehydrated.

During the temperature upgrade phase (phase 2), the LT reactor is connected to the waste heat at T m . In order to achieve the dehydration of the LT reactor and the hydration of the HT reactor, the operating condition of the HT reactor must be in the orange zone of the Fig. 1b, thus the HT reactor temperature must be inferior to T Hmax (i.e. the equilibrium pressure of the LT salt at T m have to be superior the equilibrium pressure of the HT salt at T H ). Therefore, this phase allows to supply heat at high temperature (T H ). Finally, at the end of this phase, the reactive LT salt is fully dehydrated and the HT salt is fully hydrated, requiring the beginning of a new regeneration phase. From these considerations, a simple numerical tool was developed to perform a first selection for LT and HT salts by checking whether the temperature conditions are satisfied (T L < T Lmax and T H < T Hmax ). 77 pairs were identified from 66 initial reactions.

A second step to discriminate between possible pairs is to consider the achievable performance in terms of COP and maximum temperature lift.

Using the thermodynamic considerations previously presented, neglecting heat accumulation by the salts (heat capacities are low compared to reaction enthalpies) and assuming that there are no heat losses, the theoretical maximal COP of the 2-salt THT is given by:

(T m ) (T m ) T L max T m T H max
T H =150°C T L =30°C ΔT lift Hydration LT sink (Q LT ) ΔT lift max a) b)
Eq. ( 4)

With the reaction enthalpy of the salt n (in J/mol v ). The maximal temperature lift is achieved for adiabatic operation and therefore corresponds to the theoretical maximum lift attainable by the HT salt during the temperature upgrade phase: T Hmax -T m . For the targeted high and waste heat temperatures (150 & 90 °C), a minimum temperature lift of 60 K is then required. To take into account the temperature pinches required in the HT reactor for the heat transfer to occur in finite time, only pairs achieving a minimum temperature lift of 80 K are considered on Fig. 2. From these results, 21 of the 77 promising reaction couples are kept. A large discrepancy between ΔT max values is observed, from 80 K to 125 K (for the CaCl 2 /CuSO 4 /H 2 O pair) with a COP max ranging from 0.39 to 0.52. Besides the static thermodynamic considerations, dynamic aspects have also to be taken into account in order to estimate specific heating power.

DYNAMIC MODEL

Governing equation

Several modeling of solid/gas reactor, with different approaches (local, global) and complexity levels (1/2D, stationary/dynamic, …) have been presented in the literature (Castaing-Lasvignottes and Neveu, 1997; [START_REF] Marias | Thermodynamic analysis and experimental study of solid/gas reactor operating in open mode[END_REF][START_REF] Nagel | Multi-physical continuum models of thermochemical heat storage and transformation in porous media and powder beds -A review[END_REF]. Most of the reported works refer to cooling or heat storage applications and only few models deal with THT applications, especially for those based on 2-salt cycles. Thus, from modeling approaches of a solid/gas reactor developed by [START_REF] Lebrun | Models of heat and mass transfers in solid-gas reactors used as chemical heat pumps[END_REF], a simple and original model of 2-salt THT reactor has been developed in the frame of this work. It should allow to study the dynamic behavior to estimate its performances, mainly heating capacity, and temperature lift. The following assumptions are made: a) Uniform temperature (T SLT and T SHT ) and local thermal equilibrium between the gas and solid phases is assumed in the reactors. b) Heat accumulation in the gas phase and heat losses are neglected. c) The heat capacities are considered constant during the reaction. d) In closed thermochemical systems (working with pure gas), the mass transfer usually does not limit reaction kinetic. Thus, the vapor pressure is assumed uniform inside the reactor. e) The assumptions (a) and (b) imply that the salt reaction advancement (X S ) is uniform in each reactor. X S equals 0 when the salt is hydrated and 1 when is dehydrated f) The kinetics of the solid/gas reactions depends only on the reaction advancement and on the equilibrium drop. This drop corresponds to the difference of the operating temperature and pressure to their value at equilibrium.

Energy balance

Couples HT Salt | LT Salt

According to assumptions (a), (b) and (c), the energy balances of the HT and LT salt reactors are respectively written as:

[ ( ( ) ) ] ̇ ̇ (
) Eq. ( 5)

with: ̇ ̇ Eq. ( 6)

With ̇ the source ( ̇ in dehydration phase) or sink ( ̇ in hydration phase) of vapor due the reaction of the salt n. The outlet temperature of the heat transfer fluid is given by: ̇ Eq. ( 7) 3.1.2 Kinetics of the solid/gas reaction According to the assumption (f), the reaction kinetic depends on the reaction advancement and on the equilibrium drop. [START_REF] Mazet | Analysis and experimental study of the transformation of a non-isothermal solid/gas reacting medium[END_REF] have shown that it can write as a 1st order kinetic law:

( ) Eq. ( 8
)
where: ( ) in hydration phase and in dehydration phase.

Operating vapor pressure

The vapor operating pressure of the interconnected HT and LT salt reactors, p v , is defined such as the vapor mass flow rate caught/released by the HT salt reactor is equal to the vapor mass flow rate respectively released/caught by the LT salt reactor. Thus, using Eq. ( 6) and Eq. ( 8), the vapor pressure writes:

( ) ( ) ( )
Eq. ( 9)

Note that in case of the kinetic coefficients of both reactive salts are equals and if the salts masses are chosen in a balanced way (i.e. in a way that the overall of both reactive salts could react, cf. § 3.2), we have A = -1. Therefore, the operating vapor pressure is equal to the average of the equilibrium pressure of each salt:

( ) ( )
Eq. ( 10)

Model parameters

Reactor heat transfer coefficient (UA)

The global heat transfer coefficient of the reactor has been identified from experimental results of [START_REF] Esaki | Experimental evaluation of the heat output/input and coefficient of performance characteristics of a chemical heat pump in the heat upgrading cycle of CaCl 2 hydration[END_REF] obtained with a prototype of THT using CaCl 2 /H 2 O as working pair. Although UA depends on the reactor geometry and on the reactive material used, we will suppose in this paper, that the change of the reactor scale and reactive material does not affect the heat transfer characteristics. The following specific heat transfer coefficient has been evaluated:

Eq. ( 11)

Kinetic parameter

The kinetic constant, k cin , lumps the heat and mass transfers in the solid pores and the intrinsic kinetics of the reaction. A value of 68.10 -4 s -1 has been identified from an experimental hydration of a THT using CaCl 2 as reactive salt [START_REF] Esaki | Experimental evaluation of the heat output/input and coefficient of performance characteristics of a chemical heat pump in the heat upgrading cycle of CaCl 2 hydration[END_REF].

As shown by [START_REF] Michel | Comparison of closed and open thermochemical processes, for long-term thermal energy storage applications[END_REF], this value indicates that the kinetic is not the main limitation of the thermochemical reaction. In the following, this constant value of k cin has been used for all the studied working pairs.

Mass of salts

In order to avoid that one of the two reactive salts limits the reaction of the other one, the mass of each salts should be chosen in a way that both reactive salts could totally react. Thus, each salt must allow to catch (respectively release) the water vapor released (respectively caught) by the other salt:

Eq. ( 12)

Model validation

The model has been validated using experimental results of a CaCl 2 •2H 2 O THT prototype [START_REF] Esaki | Experimental evaluation of the heat output/input and coefficient of performance characteristics of a chemical heat pump in the heat upgrading cycle of CaCl 2 hydration[END_REF]. Fig. 3 presents simulations that fit adequately to experimental results. The difference between experimental and simulated reaction rates and outlet heat transfer fluid temperature is rather low: at X=0.5, it is respectively of 9% and 0.5%. These results demonstrate that this simple model predicts in a satisfactory way the thermal behavior of a reactor.

SALT SELECTION CONSIDERING DYNAMIC ASPECTS

In a real cycle an equilibrium drop must be considered in order to realize the thermochemical reaction. This parameter influences the heat and mass transfer rates and therefore the reaction kinetics [START_REF] Mazet | Improving the performance of the reactor under unfavourable operating conditions of low pressure[END_REF][START_REF] Obermeier | Thermodynamic analysis of chemical heat pumps[END_REF]. Thus, any change in the operating temperature (T L , T m or T H ) results in a change in the equilibrium drop and then in the reaction kinetic. However, as indicated by [START_REF] Obermeier | Thermodynamic analysis of chemical heat pumps[END_REF], the equilibrium drop is not an independent parameter and equal values do not produce the same heat and mass transfer rates for different reaction systems. This parameter can be related to the chemical potential difference which represents the thermodynamic driving force of the reaction. Usually, the greater the driving force, the greater the reaction kinetic is. Under the assumptions of pure substances and ideal mixtures, the driving force can be described by the Gibbs free energy of the reaction [START_REF] Obermeier | Thermodynamic analysis of chemical heat pumps[END_REF]. Thus according to the eq. ( 2) and (3) this parameter can be writes as:

Eq. ( 13)

In order to take into account the thermal resistances inside the reactors of a 2-salts THT, the reactor temperature must be evaluated for each reactor. Thus, assuming that the heat accumulation is negligible compared to the reaction enthalpy, the eq. ( 6), ( 7) and ( 8) lead:

Eq. ( 14)

Thus, assuming a reaction advancement of 0.5, determining the vapor pressure using eq. ( 10) and imposing a heat transfer fluid pinch (5°C for the HT reactor and 1°C for LT reactor), the temperature of each reactor has been determined and then the reaction driving forces have been evaluated.

CYCLE PERFORMANCES

This part aims at analyzing the dynamic behavior and performances of a 2-salt THT. Three key parameters allowing to define the performance of the system are studied: the system heating power (during the upgrading phase), the COP and the temperature lift. The cyclic behavior of the two LT and HT reactors is simulated with successive heating and cooling. The results are given at cyclic steady state and optimal cycle time (maximum heating power).

Dynamic behavior

The dynamic behavior for CaCl 2 /Ca(NO 3 ) 2 and CaCl 2 /BaCl 2 reactive couples is presented Fig. 5a andb. Maximal specific heat power respectively of 325 W/kg and 300 W/kg have been reached during the temperature upgrade phase, corresponding to average specific heat power of 49 W/kg and 26.9 W/kg. Furthermore, due to the system regulation, the temperature and specific heat power are not stable, especially during the temperature upgrade phase.

In agreement with the driving forces presented in the Fig. 4, the reaction advancement evolution shows that the reaction kinetic of each phase is quite equilibrated for the CaCl 2 /Ca(NO 3 ) 2 THT, while the temperature upgrade phase is much faster than the other phase for the system using the CaCl 2 /BaCl 2 couple (30 min instead of 1.5h). In this last case, the system performances are therefore limited by the regeneration phase, especially due to the adverse operating conditions.

Couples HT Salt | LT Salt

Sensitivity to operating temperatures

Simulations have been performed to study the sensitivity of the operating temperatures. Their influence to the average specific heat power and the reaction driving forces is presented on Fig. 7.

The impact of the operating temperatures clearly depends to the driving force values. Thus, although the diminution of T L increases the driving force associated to the generation phase, no significant change of the specific heat power is observed. As contrary, the increase of T H (and then of the driving force associated to the temperature upgrade phase) allows in a first time to significantly increase the system performances. This is because in the reference case, the temperature upgrade phase corresponds to the lower driving force value. This phase is therefore the limiting phase, which is subject to the most unfavourable conditions on which it is necessary to play in priority.

Other limitations, as the heat transfer or the kinetic reaction can impact the system performances. However, this is out the scope of this paper. Note that in the studied cases, these parameters are weakly limiting.

CONCLUSION

A screening of salt couples was performed to design a 2-salt Thermochemical Heat Transformer for waste heat recovery. For the sake of clarity a representative case was considered with a low temperature source at 30 °C, an intermediate at 90 °C and a high temperature source at 150 °C.

A first tool was developed to discriminate pairs using a static thermodynamic approach insuring cycle feasibility. From the selected 77 pairs, 21 were kept considering a minimum requirement in temperature lift. The last step was based on dynamic consideration (driving force) and allows to identify three candidates: CaCl 2 /Ca(NO 3 ) 2 , CaCl 2 /SrBr 2 and CaCl 2 /K(OH). Cycle operation was simulated for the most promising pair CaCl 2 /Ca(NO 3 ) 2 . An average specific heating power as high as 49 W/kg is achieved.

ACKNOWLEDGEMENT

The authors wish to thank the INSA LYON (BQR FATAL RECUP) for funding this study.

NOMENCLATURE

Figure 1 :

 1 Figure 1: a) Working principle of a two-salt THT. b) Clapeyron diagram representation of the operating areas during hydration and dehydration for the CaCl 2 /SrBr 2 pairs.

Figure 2 :

 2 Figure 2: Maximal COP versus theoretical maximal temperature lift, for T LT =30°C and T m =90°C.

Figure 3 :

 3 Figure 3: Experimental and simulated results of the hydration of the THT (A type reactor, (Esaki et al., 2017)) a) Reaction advancement vs. time during. b) Outlet heat transfer fluid temperature vs. time.

Fig. 4

 4 presents a comparison of the driving force of the hydration of the LT and HT reactors for different reactive couples (i.e. a comparison of the driving force during the regeneration and the upgrading temperature phases). A large discrepancy between F m values is observed, from 0.001 to 0.35. However, promising reactive pairs are identified, as CaCl 2 /Ca(NO 3 ) 2 , CaCl 2 /SrBr 2 and CaCl 2 /K(OH). They present high driving forces whose values are close between each cycle phase. This allows to avoid large reaction times gap between these phases.

Figure 4 :

 4 Figure 4: Driving force for the hydration of the LT salt versus driving force for the hydration of the HT salt (for X Sn =0.5, T L =30°C, T m =90°C, T H =150°C, ΔT htf,HT =5°C and ΔT htf,LT =1°C).

Figure 5 :

 5 Figure 5: Specific heat power and reaction advancement versus time of 2-salt THT operating with: a) CaCl 2 / Ca(NO 3 ) 2 working pairs. b) CaCl 2 / BaCl 2 working pairs.

Figure 6 :

 6 Figure 6: Impact of the operating temperatures. a) Evolution of the driving forces. b) Specific heat power release during the upgrading temperature phase. Case of the CaCl 2 /Ca(NO 3 ) 2 reactive couples.


C Bernay-Angeletti

J-M Favreau

Technical Report : 3D mesh Segmentation Application to glasses

Object segmentation in meaningfull parts is an important research topic as these have become a part of many mesh and object manipulation algorithms in computer graphics, geometric modelling and computer aided design. In particular, tasks that were before realized by experts are more and more assumed by people without particular formation in geometric or computer graphic problems. So there is a need to have some ergonomic tools that have expert know-how while user keep control on specific constraints to its particular application. Here, we present how for segmentation purpose we model expert know-how and two algorithms that can be used to make segmentation task. One is generic and the other is specialized for glasses segmentations.

Introduction

In reason of 3D printers democratisation and easy access to 3D scanners, there is more and more 3D meshes available on internet. So use of 3D mesh becomes commonplace, increasing needs in automatic processings and in interfaces for computer-aided manipulations. In particular, tasks that were before realized by experts are more and more assumed by people without particular formation in geometric or computer graphic problems. So there is a need to have some ergonomic tools that have expert know-how while user keep control on specific constraints to its particular application. Mesh segmentation has become a key point in many computer graphics tasks and applications and also in geometric modelling. Segmentation can help for parametrization, mesh editing or animation and many other tasks [START_REF] Chaudhuri | Probabilistic reasoning for assembly-based 3d modeling[END_REF][START_REF] Chen | A Benchmark for 3d Mesh Segmentation[END_REF]. So segmentation is a crucial point when we want to offer tools for automatic processing.

Mesh segmentation partitions the surface into a set of clusters under some criteria often depending of the final application. Comparative survey on segmentation technics [START_REF] Attene | Mesh segmentation-a comparative study[END_REF][START_REF] Shamir | A survey on Mesh Segmentation Techniques[END_REF][START_REF] Theologou | A comprehensive overview of methodologies and performance evaluation frameworks in 3d mesh segmentation[END_REF] classify them in two categories: geometric segmentation and semantic segmentation. We are only interesting in the second one: our aim is to extract meaningful parts (having experts knowledge about what that mean) of an object in reasonable time, as it will ba a pretreatment in a more complete application. Different parts should be well identify because different treatment will be allowed depending on their nature. Automatically identify meaningful parts of a 3D mesh is challenging. A common practice is to build some specific features by extracting some geometric or volumetric properties like approximate convexity anlysis [START_REF] Van Kaick | Shape segmentation by approximate convexity analysis[END_REF], concavity-aware field [START_REF] Oscar | Mesh segmentation with concavity-aware fields[END_REF], spectral clustering [START_REF] Liu | Segmentation of 3d meshes through spectral clustering[END_REF], K-means [START_REF] Shlafman | Metamorphosis of polyhedral surfaces using decomposition[END_REF], core-extraction [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF], graph-cuts and randomized cuts [START_REF] Golovinskiy | Randomized cuts for 3d mesh analysis[END_REF], randomwalk [START_REF] Lai | Fast mesh segmentation using random walks[END_REF] and then use it for segmentation purpose. Part boundaries or frontiers are faces which have at least one neighbor with a different label (that is which belongs to another part than theirs). Inducing part boundaries along negative curvature minima is often used even for semantic segmentation because they often offer natural separation between parts like legs and torso. Volumetric considerations can also be used and are more related to skeletal shape representation like shape diameter function [START_REF] Shapira | Consistent mesh partitioning and skeletonisation using the shape diameter function[END_REF].

Rathers than trying to segment a single object some recent researches are done on co-segmentation [START_REF] Hu | Co-Segmentation of 3d Shapes via Subspace Clustering[END_REF][START_REF] Meng | Unsupervised cosegmentation for 3d shapes using iterative multi-label optimization[END_REF][START_REF] Shu | Unsupervised 3d shape segmentation and co-segmentation via deep learning[END_REF][START_REF] Sidi | Unsupervised co-segmentation of a set of shapes via descriptor-space spectral clustering[END_REF][START_REF] Yang | Unsupervised Co-segmentation of 3d Shapes via Functional Maps[END_REF]: the aim is to segment in the same time a family of objects (for example different glasses) because more semantic knowledge can be extracted by simultaneously analyzing a set of shapes instead of a single shape. But their final objective is shlightly different of ours as they only guarantee consistent segmentation over the set whereas we want to use expert knowledge about what is for a given problem a good segmentation. Indeed, often saying what is meaningful parts or good segmentation is not a so easy task. Good exemples of that can be found in [START_REF] Chen | A Benchmark for 3d Mesh Segmentation[END_REF] where users are asked to segment different objects. We can see that different persons segments same object differently (see figure 1). Why one should be consider better than another?

That why, we consider that we have at least an exemple of good segmentation in the sense that it has be done by an expert and so it can be considered wellsuited for its specific problem. In consequence, here, expert knowledge will be represent by at least one shape and its good segmentation and our aim will be : given another shape of the same family, how can we segment it so that extracted meaningful parts are coherent with example shape. Figure 1: An exemple of segmentations of the same objet, here glasses, from different persons. Segmentations can be found in http://segeval.cs. princeton.edu/, details of segmentations process are described in [START_REF] Chen | A Benchmark for 3d Mesh Segmentation[END_REF].

Strategy

Context

Our project implies at the same time a local company Polymagine and a research laboratory LIMOS. Company needs directed areas of research so that proposed solutions can be used in industry in the long run. In particular, we will take care about time computations and will try to have few constraints for meshes (we will not fixe number of vertices or faces for example). Polymagine is specialized in personnalisation of objects. Polymagine is interesting in segmentation as a pre-processus before personnalisation step : deformation and transformation cannot be apply in the same way on all the object. They represent the expert knowledge and provide us some segmentations examples for glasses (see figure 2). So, glasses are for us composed of tenon (in red) on which will come to hang on the legs of glasses, circles (in blue) which hold glasses and finally bridge (in green) which connects both circles. We assume that we have at least one labelled example of glasses but not a lot.

Learning with (deep) neural network

Deep learning is widely used in image segmentation process. AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] is one of the more known deep neural network. It is one of the pioneer in deep CNN. Visual Geometry Group (VGG) [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] and GoogLeNet [START_REF] Szegedy | Going deeper with convolutions[END_REF] are two others very well known deep neural network for image segmentation. They give good results. Even if CNN (Convolutionnal Neural Network) cannot be used for mesh segmentation because of the complexity of 3D geometry (a vertice don't have all the time the same number of neigboor whereas in an image the number of pixels surrounded one pixel is fixed), it seems to us that neural network can be interesting in our case, as we want to transfer knowledge from given examples to new examples of same family.

Features

When using neural network, what are entries is one of first question to answer. As previously stated, vertices (or faces) of meshes are not good candidates because they can vary in number and it is difficult to order them : what can be consider as first vertice and how we traversed the others ? So we make the choice to work as in more classical mesh segmentation technics with particular features.

We find that the work of [START_REF] Yang | Unsupervised Co-segmentation of 3d Shapes via Functional Maps[END_REF] is interesting. Their aim is not exactly the same as they want to make co-segmentation and so do not have a model to respect in segmentation. They compute presegmentations (rough segmentations) of a family of meshes using K-means in space of eigenfunctions of Laplace-Beltrami [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF][START_REF] Meyer | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF]. Then they use these presegmentation and properties of Laplace-Beltrami operator to compute Functional Maps [START_REF] Ovsjanikov | Functional maps: a flexible representation of maps between shapes[END_REF] and to make correspondance between parts of differents objects. Then, they are able to make segmentation of their family of meshes. Results on differents family of shapes are impressive.

We make some preliminary tests using the code for Functional Maps provided on http://www.lix.polytechnique.fr/~maks/fmaps_course/publications. html. To use this code, we need to have a pre-segmentation of meshes, that is a segmentation in some coarse regions and we need meshes with the same number of vertices. This two conditions are for us only for tests, because we want to be able to treat meshes from our industrial partner that will not respond to this criteria (meshes tend to have the same order of vertices but not exactly the same number).

On figure 3, tree pre-segmentations are showns. On figure 4, results of Functionals Maps given the previous pre-segmentation are shown. Colors indicates similar regions. On the first sub image, the second glasses (the upper one) has been created by deformation of the first glasses (on bottom). Matching between the two shapes is not perfect in particular near the bridge but it seem to be promising. At contrary, on second subimage, the two shapes have independant meshes and matching is really bad.

Functional Maps not seem to be well-suited for our problematic were meshes are not near isometric deformation one from another. Nevertheless, we keep being interesting in the Laplace-Beltrami Operator as features. Indeed Laplace-Beltrami is a local descriptor where different feature information are represented by different feature functions. Since eigenfunctions of Laplace-Beltrami are order from low-frequencies to high-frequencies. Yet low-frequencies contain more information than high-frequencies functions. So a few of Laplace-Beltrami eigenfunctions can be used as features descriptors. Advantage of Laplace-Beltrami eigenfunction are that they are independant from shape spacial localisation and that it describe well shape geometry. In figure 5, the 4 first eigenfunctions of Laplace-Beltrami are shown.

On experiments, there can be some problems of flips in Laplace-Beltrami eigenfunctions : sometimes sign of similar regions are inverted. An exemple is given on figure 6. On the first glasses, high value (in yellow) are on bottom-left circle and upper-right circle whereas on the second glasses it is on bottom-right circle and upper-left circle. To overcome this problem of flip signs, we simply consider the absolute value of the Laplace-Beltrami eigenfunctions. We will also normalize eigenfunction to be between 0 and 1 :

φ * i = |φ i | -min j |φ i (x j )| max j |φ i (x j )| -min j |φ i (x j )| (1)
where :

• φ * i : is the ith feature used. It is the ith eigenfunction of Laplace-Beltrami in absolute value and normalized.

• φ i is the ith eigenfunction of Laplace-Beltrami.

• x j is a vertice of the mesh

• | -| represents the absolute value.
In fact, we are more interested in having one label per face than having one label per vertice so we convert this vertice features in face features. To make that, for a face f , we consider the 3 vertices which compose it (as we work on triangulated mesh) and we say that the feature face value is the mean of the 3 features vertices.

φ * i (f k ) = 1 3 3 j=1 φ * i (v k j ) (2)
where :

• φ * i (f k ) is feature vector computed for face f k and we will refer to it by absolute-normalized eigenfunctions of Laplace-Beltrami for face f k .

• φ * i (v k j )
is the absolute-normalized eigenfunctions of Laplace-Beltrami for a vertice v k j that is one of the 3 vertice of the face f k .

By default, when we will refer to this face feature, we will named these absolute normalized eigenfunctions of Laplace-Beltrami (and we will ommitted for face).

Another problem with the Laplace-Beltrami eigenfunction is an ordering problem : even if all regions can be well characterized by these functions, sometimes it is not the same function that represent correspondant regions. Whereas sign flip can be easily removed by taking the absolute value of eigenfunction, there is no solution for this ordering problem. Or order has big importance when we want to use neural network. A node should always represent the same component. Or if there is change in order, it will not be possible. An exemple of this problem is illustrate in figure 7. For illustration purpose, we have change objects and show images on animals but problem is the same on glasses. Here, leg regions are represented one time by the 3rd eigenfunction of Laplace-Beltrami and on other on the 4th eigenfunction.

Construction of our neural network

To construct our neural network, we have to take into account to the previous remarks and to industrial time constraint : we want that segmentation process can be used in some application so that should not take to much time. As we are only developping some prototype, a time of few minutes (5 minutes should be the maximum and if possible less than one minute) is the goal fixed by Polymagine, considering this computation time could be improved if technics give good results. Another important constraint is that we want to treat any mesh given so number of vertices or faces cannot be fixed and our neural network must be able to work with. Often, when size of example if not fixed, RNN (Recurrant Neural Network) are used, in particular for text recognition. But, we think it was not suitable for us, as a sequence has natural order : place of words in sentence (or place of letters in words) are important : verb is after subject and so on. In meshes, ordering is difficult, there is no natural beginning or end of meshes. That why, we decide to take in entry only feature for one face. Doing this, all size of meshes can be treated, if a mesh has n faces, we simply need to use n times our neural network only changing entries. As we give only one face after another, we cannot ensure topological consistence from neural network but we think that absolute-normalized eigenfunctions of Laplace-Beltrami contain enough information about face neighbourhood to ensure coherent segmentation.

The main problem that remains is an ordering problem. As discussed previously, absolute-normalized eigenfunctions of Laplace-Beltrami can have different order in different meshes. So we cannot make training on some glasses and then used trained neural network on new glasses. Our main idea, that is really different from usual use of neural network, was that we do not need to learn on a lot of faces and that only few faces can be used to make training. Idea is to have for one particular mesh few labelled faces per class that are well positionned and to transfert this face (like some initial source seeds) to another glasses. Once, we have for new glasses, some labelled faces we can make training with our features faces for new glasses on this few labelled faces (on new glasses) and then we only need to used the trained neural network on all faces of new glasses. As training is on features from new glasses, there is no problem of order (in one mesh, order is always consistent).

With this idea, we choose to construct a neural network with only one hidden layer considering that as we want to train with few data, we will not have enough data to learn a lot of weighs in a lot of hidden layers. The other good point, in having a small neural network is that training step will be fast (as testing step). We have a fully connected neural network with 50 entries (the 50 first absolute normalized eigenfunctions of Laplace-Beltrami) , 3 outputs (one per class) and 250 nodes on hidden layer. Activation functions from entries to hidden layer are sigmoïds and activation functions from hidden layer to output are softmax functions.

Method 1 : Seed Transfert

In this section, we will present in more details our method using the previous idea.

First Tests

During the analysis phase, we made hypothesis and thus first step is to verify them. We begin by trying to learn segmentation on mesh when we manually select the seeds. Good selection of seeds imply :

• every class should have the same number of seeds. Seeds must be balanced in number.

• in every class, seeds must be spacially distributed. Seeds should be as far away from other seeds of same class as possible.

• in every class, some seeds must be near boundaries : if we segment every face in mesh, some seed must be near a face of different label (this face is not necessary a seed even for the other class).

We compute for all face of meshes, absolute-noramlized eigenfunctions of Laplace-Beltrami. We use 10 seeds per class (for a total of 30 seeds) to train our neural network and after we use the trained neural network for testing segmentation on all faces. Figure 8 shows results for one glasses. Segmentation is coherent with what we want and so it seems that our idea can work.

Seeds Transfer

To have an automatic process, we want to have an algorithm that is able to tranfer seeds from source mesh (we know source seeds labels) to target mesh. Figure 8: Segmentation from manual seeds using our neural network We want to label the same number of faces in target mesh than in source mesh. Labelled faces in taget mesh will become seeds for target mesh.

Distance Transfer

Firstly, we consider we have aligned glasses. Under this hypothesis, tranferring seeds can be seen as a research of k-nearest neighbour where k = 1, tree is constructed using target mesh face (the center of face is taken as 3D point) , query points are center of seeds in source mesh. The target face that is selected become a seed for target mesh with same label as source seed.

Distance Transfer + eigenfunctions

As glasses have not all exactly the same shape, distance transfer can be improved by taking into account more topological information. For every seeds in source mesh, we will find the k-nearest neighbour where k = 250. This ensure a spacial coherence. But after to select the final target face, we compute distance between the 3 first eigenfunctions of Laplace-Beltrami, the face with the lower distance is kept as target seed. The idea is that eigenfunction of Laplace-Beltrami will probably be more similar on place that have similar function in glasses.

Results of the two transfer can be seen in 5.1.

Method 2 : Training

We are also interesting to see if we can use our neural network to make some more usual training as the context change a little bit and more labelled glasses are available. Like we previously discuss, using absolute-normalized eigenfunction of Laplace-Beltrami is not suitable for this purpose so we have to change our features.

We choose to use HKS (Heat Kernel Signature) [START_REF] Sun | A Concise and Provably Informative Multi-Scale Signature Based on Heat Diffusion[END_REF]. As previously, we will take it in absolute value, normalize it and compute a mean to have features for face. So we will named it absolute-normalized HKS. We believe theses features will be less sensitive to the ordering problem. We choose to take 50 features.

We used k glasses to make training. On every glasses we randomly select n faces for every class. So we have kn examples in training set. We train our neural network (the one described in 2.4) only one time with these examples (training can be long). After, we compute for every new glasses, its absolutenormalized HKS features and use the trained neural network on every face by changing entries. Results are shown in 5.2.

Results

In this section, we will present some results obtained from the two methods. To implement our algorithmes, we use Matlab to compute generealized eigenvectors and values, Tensorflow in Python to construct, train and use neural network and C++ to compute features on mesh and make file readable by Matlab or Python. For a mesh with near 23914 vertices and 47832 faces, we need to :

• compute mass matrix and cotangent matrix : 1.4 s

• compute generalized eigenvectors and eigenvalues : 21.6 s

• select seeds 4.9s using distance only or 7.6 s using distance + eigenfunctions

• compute normalized-absolute eigenfunctions of Laplace-Beltrami using previous steps and writing file test : 12s

• train our neural network : 3.4s

• use trained neural network on all faces : 6.5s

Result for the first method

Conditions : on one glasses, 10 faces (seeds) per class are well selected. Selected face are well balanced in every class region and somme of them are near boudaries. All meshes are aligned. As describes in section 3.2, we use two methods for transfert. The first one take into account only euclidean distance between faces, the second one select for each source face, the 100 nearest target faces and among this 100 target faces, the one with the most similar 3 first eigenfunction is chosen. Figure 9 shows one case of mislabelling in transfert in both methods. On the left, it is the distance transfert and on the right it is the distance + eigenfunction transfert.. We can see error in labelling in left tenon, more visible in the second line.

In figure 10, we compare the results on both transfert on 12 glasses. First column reprensents results using transfert with distance only (and after training with seeds and testing) and second column represents results using transfert with distance and eigenfunction (and after training with seeds and testing). Results with distance and eigenfunction give a little bit better results. But two are not perfect segmentation and some glasses are bad segmented. The problem comes from error in labelling as first cause, and position of labelled seeds : they are not as well balanced as in source mesh (they can be less distant or it can lack of seed near boundaries).

Result for the second method

Conditions : we use one glasses as training glasses. We consider we have a good segmentation for all the faces of the mesh. Among this faces, we randomly select 750 in every class, no particular care is done about the repartition of selected faces. Absolute-normalized HKS is computed for each face of all the meshes we want to test. Computing it take 72.3 s. Neural network is the one described in 2.4. Training take 196 s.

In figure 11, the result on the training glasses is shown. Mesh has near 25000 faces. Note that, segmentation is good even if we have train with only 2250 faces.

In figure 12, results on 12 others glasses are shown. Even if segmentation are not perfect, they are interesting and should be improved by modification on training : we can use more than one glasses and we also can select faces in a better way, maybe trying to have a good repartition on class as in the first method.

Conclusions

In this report, we presented two algorithms for mesh segmentation using Laplace-Beltrami Operator to compute interesting features and neural network to learn from examples.

First Method : Seed Transfert

The first method transfers few seeds (a dozen per class is enough) from a source glasses to a target glasses. Seeds must be distributed among faces of similar class and some must be close to faces of others class (so near border of semantic region). Two methods of transfert are proposed. Both consider that glasses are aligned. The first use only euclidean distance between source face and target faces to affect a label. The other assume that there is no problem of order between the three first eigenfunction of Laplace-Beltrami. The 100 nearest facesin target mesh of a source face are selected and in these pre selected face, the one that has the most similar eigenfunctions is chosen.

Once seeds have been transferred, a neural network is trained. In entries, we used the 50 first absolute-normalized eigenfunctions of Laplace-Beltrami, and we have 3 outputs (one per class) that are used to affect label. There is only one hidden layer as we train we very few data (only 10 examples per class). Once the neural network is trained, for each face, we compute tthe absolute-normalized eigenfunctions of Laplace-Beltrami and use the neural network with this entry.

Si if a mesh have n faces, the neural network will be used n times. All this process take near one minute and computed time can be optimized. So it can be used for industrial purpose.

Having good initial seeds give good segmentation. For the moment, results are promising but need to be improved. The limitation step is in the process of transferring seeds. For industrial application, it is possible to imagine to have some graphical interface and let the user select seeds. The main advantage of this manual selection is that the quality of the seed should be better and so will lead to better segmentations. Another advantage is thatwith this manual selection no more alignment will be requirred. The main inconvenient is that ask more work from users.

Second Method : Big Train then Test

The second method use absolute-normalized HKS as features. The neural network is the same as previously but it is trained only one time with more examples (we have tests with 750 examples per class chosen on only one glasses but taking more examples and more glasses should improve results). The training step can be long but it is done only one time at the early beginning. After the training step, all we need is to compute absolute-normalized HKS for each face and use the neural network with these features in entry. We obtain segmentation for others glasses. The process (without training) take less than two minutes and computed time can also be improved. So this method can also be used for industrial purpose.

For glasses, it seems to be very promising but there is some restrictions to keep in mind. The first is that HKS is sensitive to size of meshes because of the use of eigenvalues, so object should be at the same order of size. The second restriction is the application to others objects : using HKS decrease order problem in eigenfunction of Laplace-Beltrami but do not solve all the problems. So, this method will not necessary give results for other shapes.

Conclusion and future works

We have developped two methods that can make segmentation using expert knwoledge in particular case of glasses segmentation. Both methods can be used for industrial purpose as the time they need is low and can be improved. Both methods are promising but results are not exploitable for the moment and need more developpments. One method should be generalized to other shape if we can improve the seed transfert. The other method will maybe give better results with less research but is probably suitable only for particular case.

Aknowledgement

This work has been sponsored by the Auvergne-Rhône Alpes region and by the European Union through the program regional FEDER 2014-2020.

Figure 2 :

 2 Figure 2: Good segmentation of glasses. Tenons are in red, circles in blue and bridge in green.

Figure 3 :

 3 Figure 3: Pre-segmentations for tree differents glasses : glasses a), glasses b) and glasses c).

Figure 4 :

 4 Figure 4: Functional Maps used on our glasses. First column shows front of glasses and second columns the back. Fisrt rows is between glasses b) and a). Second rows is between glasses a) and c). Same colors should indicate correspondant regions.

Figure 5 :

 5 Figure 5: On first row : 1st, 2nd eigenfunctions of Laplace-Betrami. On second rpw : 3rd and 4th eigenfunctions of Laplace-Beltrami.

Figure 6 :

 6 Figure 6: the 4rd eigenfunction of Laplace-Beltrami for two differents glasses. We can observed change in sign.

Figure 7 :

 7 Figure 7: eigenfunction of Laplace-Beltrami that highlights the back legs of the two animals. For the pig, it is the 4th eigenfunction whereas for the veal it is the 3rd.

Figure 9 :

 9 Figure 9: Example of transfert labelling. In light grey it is not labelled face, in blue, it is tenon faces, in green bridge faces and in red circle faces. First line show labelled faces on full mesh and second line is a zoom on the left tenon to show error in labelling process. Fisrt column is labelling using transfert with only distance. Second colum is labelling using transfert with distance and eigenfunctions.

Figure 10 :

 10 Figure 10: First column represents results after transfert with distance only, training and testing. Second columns represent results after transfert with distance and eigenfunction, training and testing.

Figure 11 :

 11 Figure 11: Result of segmentation on training glasses for the second method (general training + HKS).

Figure 12 :

 12 Figure 12: Results of second method (general train and HKS on different meshes of glasses


domain: spi phys

text:
TEST OF A LWIR CRYOGENIC MULTISPECTRAL INFRARED CAMERA FOR METHANE GAS LEAK DETECTION AND QUANTIFICATION IN TADI PLATFORM AND ERSPERCE SITE

OECD CONFERENCE CENTER, PARIS, FRANCE / 28 -30 JANUARY 2020

Guillaume Druart [START_REF] Abdel-Moati | Advanced Techniques for Autonomous Detection of Gas Releases[END_REF] , Pierre-Yves Foucher [START_REF] Abdel-Moati | Advanced Techniques for Autonomous Detection of Gas Releases[END_REF] , Stéphanie Doz [START_REF] Abdel-Moati | Advanced Techniques for Autonomous Detection of Gas Releases[END_REF] , Xavier Watremez [START_REF] Singer | Design of a cryogenic IR detector with integrated optics[END_REF] , Sophie Jourdan [START_REF] Gibbons | System and method for viewing an area using an optical system positioned inside of a dewar[END_REF] , Emmanuel Vanneau [START_REF] Druart | Towards infrared DDCA with an imaging function[END_REF] , Sylvain Favier

(5) [START_REF] Abdel-Moati | Advanced Techniques for Autonomous Detection of Gas Releases[END_REF] ONERA, BP 80100, FR-91123 PALAISEAU CEDEX (FRANCE), Email: guillaume.druart@onera.fr [START_REF] Singer | Design of a cryogenic IR detector with integrated optics[END_REF] TOTAL, PAU (FRANCE) [START_REF] Gibbons | System and method for viewing an area using an optical system positioned inside of a dewar[END_REF] LYNRED, VEUREY-VOROISE (FRANCE) [START_REF] Druart | Towards infrared DDCA with an imaging function[END_REF] NOXANT, PALAISEAU (FRANCE) [START_REF] De La Barrière | Compact multichannel infrared camera integrated in an operational detector dewar cooler assembly[END_REF] BERTIN, AIX-EN-PROVENCE (FRANCE) KEYWORDS: Infrared, multispectral, cryogenic camera, methane, detection, quantification

ABSTRACT:

Cryogenic cameras are an innovative alternative in the design of miniaturized infrared cameras using cryogenic detectors. We apply this technology to design a snapshot multispectral camera for gas leak detection. A UAV compatible demonstrator in a commercial Detector Dewar Cooler Assembly (DDCA) has been made and tested in the TOTAL Anomaly Detection Initiatives (TADI) platform and the Esperce site of ONERA. We demonstrate the ability to detect and quantify in real time the origin of methane gas leak, the flowrate and the volume of the plume with a four spectral channel LWIR camera on ground or from a UAV. During the test campaigns, methane gas leak flowrates greater than 0.5g/s were detected and quantified. Moreover, an image processing that doesn't need reference images without gas leak was used, facilitating mobile applications.

INTRODUCTION

Methane leaks in petrochemical facilities present a risk to health and safety as well as a negative impact on environment. Major leak can potentially lead to explosive environments with a risk of major industrial accident. Existing autonomous leak detection systems mainly rely on point and path detection but these technologies require that the gas contacts physically the point sensor or moves between two path detectors. The drawback of these methods is that they don't permit to detect and localize early the small gas leak rates and to quantify any leak rates. Manual leak checks are also performed by monitoring one component at a time using a flame ionization detector or an optical gas imager but are time-consuming in regards with the complexity of the facilities. Moreover, with this method, the inspection is limited in time and small leaks may go undetected for extended period. Infrared gas cameras have thus become a key technology capable of visualizing gas plumes and are able to provide a continuous monitoring of large areas [1].

TOTAL is looking for a versatile remote sensing technology that can be applied in different ranges of methane leak flowrates in three main scenarios:

• Major leaks in crisis management with flowrate greater than 100g/s. • Medium size leaks in safety monitoring with flowrate between 10g/s and 100g/s. • Small leaks in environmental monitoring with flowrate smaller than 10g/s. Moreover the remote sensor should allow autonomous detection by replacing the operator with a computer vision algorithm. It should be able to detect and quantify in real time the origin of methane gas leak, the flowrate and the volume of the plume. Finally, the camera should be lightweight to be integrated in robots or drones for autonomous inspection.

Infrared cameras are able to detect methane gas since its absorption bands appear in the Mid-Wave InfraRed (MWIR) bandwidth around 3.2µm and in Long-Wave InfraRed (LWIR) bandwidth around 7.6µm. However state-of-art cameras and commercial cameras don't fully answer all the specifications of TOTAL. Some are monochromatic and thus sensitive to false alarms leading autonomous detection difficult. Indeed the methane plume can be difficult to distinguish by a computer vision algorithm from common interferences such as people, cars and steam.

Other cameras are able to produce hyperspectral or multispectral images but the spectral content of the scene need to be scanned over time with a rotating filter wheel or moving mirrors. The frame rate is low and inadequate for mobile acquisition.

The FUI consortium IMAGAZ has developed SIMAGAZ, a camera for detection and quantification methane gas leak with the following specifications:

• Compact and lightweight to be integrated on a UAV • With a high sensitivity to guarantee quick acquisition time • Multispectral to quantify methane gas leak, to limit false alarm and to allow an autonomous detection with a computer vision algorithm.

To meet all these features, we took benefit of 10 years of R&D between ONERA and LYNRED on cryogenic infrared cameras. Numerous works have shown the advantages of these cameras in terms of miniaturization [2,3,4,5]. Moreover, by using a cryogenic infrared detector, we guarantee the high sensitivity of the camera.

SIMAGAZ is the pre-industrial version (TRL 1 5) of a previous laboratory demonstrator studied in a project between ONERA and TOTAL called NAOMI (New Advanced Observation Method Integration).

FEATURES OF SIMAGAZ

The SIMAGAZ camera relies on the technology of cryogenic optics directly integrated in a standard Detector Dewar Cooler Assembly (DDCA). SIMAGAZ uses the LYNRED's state of the art architecture based on his cooled LW VGA format with a pixel pitch of 15µm and a RM3 cooler. So the outside view of the camera looks like a standard DDCA. Only the cold shield has been modified to integrate the optics. The optical part is simply made of a single array of lenses to produce 4 images of the scene on a single Focal Plane Array (FPA) and an array of filters to give the multispectral property of the camera, each filter being associated to a different optical channel. The MCT technology of the infrared detector array offers high speed operation with a high sensitivity. Indeed, a Noise Equivalent Temperature Difference (NETD) has a mean value of 40mK for all the spectral bandwidths and with an integration time of 5ms and for a scene at around 20°C. Therefore, the camera combines the advantages of being simple (only two major optical components are used), cost effective (a single detection block is used), compact (the camera has the size of the 1 Technology Readiness Level (TRL) detection block), sensitive (a cryogenic infrared detector is used), as well of having a high radiometric stability (the optics are at cryogenic temperature, thus limiting the background signal). The specifications of SIMAGAZ are summarized in Table 1 The core camera uses the processing unit for UAV airborn application of the NoxEngine technology platform from NOXANT. The NoxEngine combines innovative software architecture with specialized processors, so images can be acquired and processed in real-time, in a reduced footprint and energy consumption. A picture of the core camera is given in Figure 1.

FIRST QUANTIFICATION RESULTS OF SIMAGAZ DURING TEST CAMPAIGNS

The first results of SIMAGAZ were made on the TADI (TOTAL Anomaly Detection Initiatives) platform [6] in 2018. Total's Lacq Pilot Platform is a test area for qualifying cost-effective systems designed to complement the gas detection system of a plant and provide valuable information should a gas leak incident occurs. New methodologies for the early detection of anomalies using remote observation systems including drones, robots and artificial intelligence data processing systems are currently being investigated there. The TADI infrastructure manages monitored gas leaks at flowrates from 0.3 g/s to 300 g/s, making this site unique in the world.

First SIMAGAZ was mounted on a 2-axis platinum turret as illustrated in Figure 2 and placed on a pod so that the camera can be elevated at 10m. SIMAGAZ monitored at a fixed position the test area at a distance of 80m.

Figure 3 shows some gas plume detection and quantification at different flowrates: 200g/s, 10g/s and 1g/s of methane. The quantification is given in false colour in ppm.m. Figure 3 illustrates the ability of SIMAGAZ to provide detection for large methane leak but also early leak detection of small leaks and even fugitive emissions, validating the three scenarios defined by TOTAL. For these tests, the quantification algorithms used a reference image without gas leak. That's why, the acquisition of the camera started just before the gas emission.

The measured flowrates were consistent with the controlled flowrates.

Figure 2. SIMAGAZ integrated on a turret Subsequently, SIMAGAZ was integrated in a small hexarotor UAV as illustrated in Figure 4. Test campaigns were first performed in Esperce site of ONERA and then in TADI platform. Figure 5 shows detection and quantification of small methane gas leak (under 1g/s) from a gas cylinder in Esperce site using SIMAGAZ integrated on the small hexarotor UAV at different altitude heights. Figure 6 shows detection and quantification of methane gas leak in TADI platform. These detections are important results since we demonstrate that SIMAGAZ can perform a continuous monitoring on a mobile vehicle thanks to algorithms working image by image (without reference images).

CONCLUSION

We demonstrate the ability of this four spectral LWIR camera on ground or on UAV, to detect in real time a methane gas leak with a quantification of its flow rate. In addition, the system allows an estimation of the volume of the plume and helps to investigate the origin of the leak. The core camera weights around 1kg and has a footprint of around 1L and power consumption at the cooling steady state of around 10W. Results from TADI and Esperce campaigns were presented. SIMAGAZ has proved its operability since it endured changing weather conditions without degradation of its radiometric performance. The quality of the images was stable in time and no additional calibrations were required during one-day test campaigns. The effectiveness of the quantification algorithms was provided through a co-design with the optical design that resulted in an optimal choice of filters and a cryogenic infrared camera with outstanding radiometric performances. SIMAGAZ validates the three scenarios considered by TOTAL: detection of large methane leaks, early leak detection of small leaks and fugitive emissions. In the following, three configurations of SIMAGAZ will be investigated: a fixed monitoring camera, a handheld camera for operators and a mobile camera integrated in an autonomous system.

Figure 1 .

 1 Figure 1. Core camera A Human-Machine Interface (HMI) developed during the NAOMI project has been applied to SIMAGAZ to validate the ability to detect and quantify in real-time methane gas plumes.

Figure 3 .

 3 Figure 3. Detection and quantification of methane gas plumes at different flowrates: A) 200g/s, B) 10g/s, C) 1g/s. SIMAGAZ monitored the test area at a distance of 80m.

Figure 4 .Figure 5 .

 45 Figure 4. Integration of SIMAGAZ in a hexarotor UAV

Figure 6 .

 6 Figure 6. Detection and quantification of a methane plume from a UAV in TADI platform.

Table 1 .

 1 . Specifications of SIMAGAZ

	Detector type		MultiSpectral MCT
	Spectral Range		7.2-8.5µm
	Image Size for a sub-	320x256 pixels with a
	image		pixel pitch of 15µm
	Measurement	Rate	Up to 20Hz
	(with image processing)
	Volume of the core	1L
	camera
	Weight of the core	1.3kg
	camera
	Focal length		7mm
	F-number		4
	Horizontal Field Of	40°
	View (HFOV)
	Measurement rate		Up to 20Hz
	Sensitivity of methane	<0.5g/s
	leak flowrates