Pablo Arias
email: arias@cmla.ens-cachan.fr

Mila Nikolova
email: nikolova@cmla.ens-cachan.fr

Below the Surface of the Non-Local Bayesian Image Denoising Method

The non-local Bayesian (NLB) patch-based approach of Lebrun, Buades, and Morel [12] is considered as a state-of-the-art method for the restoration of (color) images corrupted by white Gaussian noise. It gave rise to numerous ramifications like e.g., possible improvements, processing of various data sets and video. This article is the first attempt to analyse the method in depth in order to understand the main phenomena underlying its effectiveness. Our analysis, corroborated by numerical tests, shows several unexpected facts. In a variational setting, the firststep Bayesian approach to learn the prior for patches is equivalent to a pseudo-Tikhonov regularisation where the regularisation parameters can be positive or negative. Practically very good results in this step are mainly due to the aggregation stage -whose importance needs to be re-evaluated.

Introduction

In this paper we analyse the Non-local Bayesian (NLB) image denoising algorithm introduced by Lebrun, Buades and Morel in [START_REF] Lebrun | A nonlocal Bayesian image denoising algorithm[END_REF], which is based on the assumption that sets of similar patches are IID samples from a Gaussian distribution.

In recent years several works have proposed Gaussian models, or Gaussian mixture models (GMMs) as priors for image patches, achieving state-of-the-art results in image denoising and other inverse problems. In [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] a GMM with hundreds of components is used as a fixed prior for image patches. The mixture is learnt via an EM algorithm from a database of two million patches. Instead of a fixed prior [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From Gaussian mixture models to structured sparsity[END_REF][START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF] proposed to learn the GMM from the input image. Guillemot et al. [START_REF] Guillemot | Covariance trees for 2D and 3D processing[END_REF] introduced the covariance tree, a hierarchical data structure capturing the covariance of image patches at several scales.

The NLB algorithm differs from the previous approaches in that it does not consider a single model (GMM, covariance tree) for all image patches. Instead, there is a Gaussian model for each patch and its nearest neighbors. This is closely

The work of P. Arias is funded by BPIFrance and Région Ile de France in the FUI 18 Plein Phare project. P. A. is also partly founded by the "IDI 2016" project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02, the Office of Naval research (ONR grant N00014-14-1-0023) and by ANR-14-CE27-001 (MIRIAM). The work of M. Nikolova work was partially funded by the French Research Agency (ANR) under grant No ANR-14-CE27-001 (MIRIAM). related to the approaches that denoise patches by applying principal component analysis (PCA) on groups of patches [START_REF] Muresan | Adaptive principal components and image denoising[END_REF][START_REF] Dabov | BM3D image denoising with shape-adaptive principal component analysis[END_REF][START_REF] Zhang | Two-stage image denoising by principal component analysis with local pixel grouping[END_REF][START_REF] Deledalle | Image denoising with patch based pca: local versus global[END_REF].

Several extensions of NLB were considered. A multiscale version was developed in [START_REF] Lebrun | The noise clinic: a blind image denoising algorithm[END_REF] to handle the structured noise typical of compression artifacts. Extensions to video were proposed in [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF][START_REF] Arias | Video denoising via empirical Bayesian estimation of spacetime patches[END_REF]. NLB was adapted to denoise manifoldsvalued images in [START_REF] Laus | A nonlocal denoising algorithm for manifold-valued images using second order statistics[END_REF]. A Bayesian hyper-prior model on the mean and the covariance in the NLB approach was proposed for inverse problems in [START_REF] Aguerrebere | A hyperprior bayesian approach for solving image inverse problems[END_REF].

In the NLB method, for each set of similar patches an empirical Wiener filter is used to estimate the corresponding clean patches. The denoised image results from the aggregation by averaging of these estimated patches. The Wiener filter requires to estimate the mean and the covariance matrix of the a priori distribution. The method consists of two steps. In the first step the statistics of the prior are estimated from the noisy patches. The second step is similar but uses the output of the first step as an oracle to estimate the a priori statistics.

In this work we carefully dissect the NLB algorithm, focusing on the first step which is responsible for most of the denoising. On the surface, the NLB algorithm appears as a maximum a posteriori (MAP) estimation of each patch. But in fact these empirical MAP estimates are poor, as a consequence of the fact that the prior covariance estimate has in general negative eigenvalues which result in an unstable filter. The success of the algorithm resides in that the aggregation phase averages out the errors of the MAP estimates, achieving similar and often better results than the ones obtained with a more standard non-negative definite estimate of the prior covariance matrix. Although the aggregation has been recognized as a good practice in patch-based denoising [START_REF] Lebrun | Secrets of image denoising cuisine[END_REF], the results presented below show that it plays a crucial role for the success of the NLB.

The theory of the algorithm

For an original gray value image x ∈ R M1×M2 we are given an observation contaminated with centered white normal noise n with known variance σ 2 , independent of x: y = x + n.

The approach proposed in [START_REF] Lebrun | A nonlocal Bayesian image denoising algorithm[END_REF] relies on the fact that natural images exhibit patch-based self-similarity, and assumes that groups of similar patches follow a normal distribution. The NLB algorithm consists in two main denoising steps:

(1) compute an oracle image x using the data y and σ 2 ; (2) obtain the denoised image x using the oracle x, the data y and σ 2 . Each one of these steps consist of two stages: (a) denoise individually all patches using a maximum a posteriori (MAP) estimator; (b) aggregate the results of (a).

An example of the output of each step is given for the 256 × 256 parrot image, shown in Fig. 1. Some noise artifacts reside in the oracle image x. They are restored in the denoised image.

Bayesian model for patches

A patch of y (resp., x) of size √ K × √ K at pixel i = (i 1 , i 2 ) is denoted by y i (resp., x i ).
In what follows, we adopt a vector notation for patches. The likelihood of x i given the observed y i is

P(y i |x i ) ∝ exp - x i -y i 2 2σ 2 . ( 1
)
With each patch x i , there is associated a set of N most similar patches S(i) within a fixed size window. The patch x i is called the reference patch. Similarity is measured with respect to the 2 distance. It is assumed that all patches in S(i) have a normal distribution sharing the same mean m x i ∈ R K and the same covariance C x i , see [START_REF] Lebrun | A nonlocal Bayesian image denoising algorithm[END_REF]. Thus

P(x i ) ∝ exp - 1 2 (x i -m x i ) T C -1 x i (x i -m x i ) . (2)
Using Bayes's rule, P(x i |y i ) ∝ P(y i |x i )P(x i ). The MAP estimate x(i) i of the noisy reference patch y i , x(i) i = arg max

x i P(x i |y i ) = I + σ 2 C -1 x i -1 (y i + σ 2 C -1 x i m x i )
is known to provide some denoising. Using that for an invertible matrix C one has

I + σ 2 C -1 -1 C -1 = C + σ 2 I -1 and σ 2 C + σ 2 I -1 -I = -C C + σ 2 I -1 ,
the MAP estimate reads also as

xi = m x i + C x i C x i + σ 2 I -1 (y i -m x i ) . (3)
This individual MAP patch denoising is considered as the main phase of the NLB algorithm [START_REF] Lebrun | A nonlocal Bayesian image denoising algorithm[END_REF]. Remark 1. Extensive numerical tests reported in [START_REF] Lebrun | Implementation of the "nonlocal Bayes" (nl-bayes) image denoising algorithm[END_REF] have shown that the denoising formula (3) can be applied simultaneously to all of the N most similar patches, i.e. those belonging to S i . As a consequence the same patch y k can be denoised several times, if k ∈ S(i) for different i. We denote the corresponding MAP estimates as

x (i)
k . After all patches are denoised using the local MAP estimator, an aggregation stage is applied by computing for each pixel i the mean over all restored patches that contain this pixel.

The main difference of the NLB method compared to other Gaussian-based models for image patches [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: From Gaussian mixture models to structured sparsity[END_REF][START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF] is the way how the matrix C x i is learned from the image itself. The question is analysed in the following sections.

Learning the a priori model

We have to distinguish between the first and the second steps. The number of the pixels in the square patch window is denoted by K j and the number of similar patches by N j , where j = 1, 2 indicates the step.

First step. The N 1 patches S 1 (i) similar to the reference patch y i are selected using the 2 distance between y i and all patches y k in a neighborhood around y i . Let us fix at patch y i and its similar patches y k ∈ S 1 (i).

Since

x i ∼ N (m x i , C x i ) and n i ∼ N (0, σ 2 I)
and since the noise is independent from the original x, we have that

y k ∼ N (m x i , C y i ), with C y i = C x i + σ 2 I, for y k ∈ S(i).
The maximum likelihood estimates of m x i and C y i are given by the sample mean and sample covariance matrix. The ML estimate of the covariance matrix is biased. An unbiased estimator can be obtained by dividing by N 1 -1 instead of N 1 , yielding

y i = 1 N 1 k∈S i y k , C y i = 1 N 1 -1 k∈S1(i) (y k -y i )(y k -y i ) T , (4)
where the superscript T means transposed. Using that C y i = C x i + σ 2 I, the authors in [START_REF] Lebrun | A nonlocal Bayesian image denoising algorithm[END_REF] estimate C x i as

C x i = C y i -σ 2 I. (5)
This estimate is often not positive definite. There are different ways to avoid this; e.g., clipping the negative eigenvalues to zero has been used on the context of patch-based denoising in [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF]. This may also have undesirable effects, see the discussion in [16, p. 406]. New facts are given in subsection 3.3.

In spite of the theoretical background explained above, extensive experiments in [START_REF] Lebrun | Implementation of the "nonlocal Bayes" (nl-bayes) image denoising algorithm[END_REF] have shown that using the MAP estimator yields artifacts in flat areas.

Remark 2. [Flat areas]

Flat areas are detected by applying a χ 2 Gaussianity test to the set of all pixels of the patches in S(i). If a flat area is detected, the patches in S(i) are estimated as constant patches, where the constant value is the average of all pixel values of all patches in S(i).

Second step. This step is grounded on the presumption that the oracle estimate x is faithful enough to provide two informations on the unknown ground truth image. These are: [START_REF] Aguerrebere | A hyperprior bayesian approach for solving image inverse problems[END_REF] the distances between patches in the oracle are well estimated, so the set of the N 2 patches similar to y i is computed from the oracle estimate x; (2) the covariance matrices for patches provided by the oracle image x are "nearly true", so these matrices are computed using only the oracle image.

An oracle image in the first step

In this section we study in detail the first step of the algorithm.

Statistical model and filtering for 1st step patch denoising

Replacing C x i and m x i in (3) by their estimates (4), we end up with the following Wiener-like filter

x (i) k = y i + C y i -σ 2 I C -1 y i (y k -y i ) . ( 6
)
This filter is applied to all patches k ∈ S 1 (i). After simplification, x

(i) k = y k - σ 2 C -1 y i (y k -y i ).
Since C y i 0 its eigen-decomposition is of the form

C y i = U y i diag(λ y i ) U T y i . (7)
Here U y i is an orthonormal matrix (i.e. it satisfies U T y i U y i = U y i U T y i = I) containing the eigenvectors of C y i and λ y i ∈ R K ≥0 are the corresponding eigenvalues of C y i in decreasing order. The estimate of C x i in (5) reads as

C x i = C y i -σ 2 I = U y i diag λ y i -σ 2 U T y i . (8)
Thus the Wiener-like denoising in ( 6) can be interpreted as a filtering in the basis of the eigenvectors U y i of C y i :

U T y i x (i) k -y i = diag λ y i -σ 2 λ y i U T y i (y k -y i ) . (9)
where division in the diagonal matrix is vector-wise. A filtering of this form is stable and the learnt prior model is consistent if and only if C y i -σ 2 I has positive eigenvalues.

The eigenvalues of the estimated prior covariances C yi -σ 2 I. It appears that C y i -σ 2 I nearly always has negative eigenvalues even for non-flat patches, whatever the level of the noise. Table 1 gives more details for five images and different values of σ ∈ {1, 5, 10, 20, 40}. The parameter choice follows the recommendations in [START_REF] Lebrun | Implementation of the "nonlocal Bayes" (nl-bayes) image denoising algorithm[END_REF]. Patches are square; for σ < 20 one takes (K 1 , N 1 ) = (9, 27) and for 20 ≤ σ < 50 (K 1 , N 1 ) = (25, 75); as recommended, N i = 3K i and the search window is a square of size N i /2 . Flat areas are detected as described in Remark 2. As in [START_REF] Lebrun | A nonlocal Bayesian image denoising algorithm[END_REF], with each index i in the image (except at the boundary) there is associated a reference patch y i . The set of indexes corresponding to non-flat reference patches is denoted by J. Column (a) contains the least eigenvalues of C y i for each image and for each noise level σ over all patches in the relevant set J. Since the eigenvalues are ordered decreasingly, the least one is K 1 . All these least eigenvalues are much smaller that the corresponding σ 2 . Column (b) shows the percentage of matrices C y i -σ 2 I over J with at least one negative eigenvalue, which also reads as {y i , i ∈ J : g(λ K1 ) < 0}. Observe that these percentages are large, going up to 100 %. What is more, the number of negative weights in g(λ y i ) for non-flat reference patches is important: column (c) presents the average number of negative weights in g(λ y i ) taken over J. Consequences. Since C x i = C y i -σ 2 I has negative eigenvalues it is not a covariance matrix, and the learned prior P(x i ) is not a probability.

Table 1 shows that most of the noisy prior covariances estimated by the denoising algorithm have a large number of negative eigenvalues, and thus they have a significant impact in the pseudo-Wiener denoising in ( 6) and ( 9). Pseudo-Tikhonov regularization. Observe that x (i) k in ( 6) is the unique stationary point of the objective F below

y i ,i∈J λy i (K1) (b) {y i : g(λK 1 ) < 0} (c) {k ∈ J : g(λ k ) < 0} boat σ =
F (x) = x -y i 2 + σ 2 (x -y i ) T C y i -σ 2 I -1 (x -y i ). ( 10
)
Using the eigen-decomposition in [START_REF] Dabov | Image denoising by sparse 3d transform-domain collaborative filtering[END_REF] and the fact that U y i is orthonormal, F also reads as

F (x) = U T y i (x -y i ) 2 +σ 2 (x-y i ) T U y i diag(λ y i -σ 2 I) -1 U T y i (x-y i ). ( 11
)
We denote λ y i := (λ 1 , . . . , λ K ). Let u y i ,k stand for the kth column of U y i . Then

F (x) = K k=1 u T y i ,k (x -y i ) 2 + β(λ k ) u T y i ,k (x -y i ) 2 , (12)
where

β(λ) = σ 2 /(λ -σ 2 ). The function λ → β(λ) is continuous on R \ {0} with lim λ 0 β(λ) = -1 lim λ σ 2 β(λ) = -∞ lim λ σ 2 β(λ) = ∞ β(2σ 2 ) = 1 lim λ→+∞ β(λ) = 0.
The objective F is a quadratic convex-concave function. It can be seen as a pseudo-Tikhonov regularization in the basis U y i where the regularization parameters β(λ k ) can be positive or negative. Let us remind that classical Tikhonov regularization with all β > 0 entails oversmoothing; see e.g., [3, sec. 3.2]. Table 1(c) shows that F in (12) has many negative regularization parameters.

Pseudo-Wiener filtering. Equation ( 9) reads also as

U T y i x (i) k = diag λ y i -σ 2 λ y i U T y i (y k -y i ) + U T y i y i . (13)
Since U y i is an orthonormal matrix, the noise statistics of y i and of U T y i y i are the same. The weighting coefficients in the diagonal matrix in [START_REF] Lebrun | Secrets of image denoising cuisine[END_REF] are given by a function g : R >0 → R

g(λ) := (λ -σ 2 )/λ = 1 -σ 2 /λ (14)
which is differentiable, strictly increasing with lim λ 0 g(λ) = -∞, g(σ 2 ) = 0 and lim λ→+∞ g(λ) = 1. The weights g(λ) are negative for λ < σ 2 and rapidly decrease with λ 0. This g is related to β in ( 12) by g(λ) = (β(λ) + 1) -1 . Thus the components of U T y i x

(i) i in (13) satisfy u T y i ,k x (i) k -y i = g(λ k ) u T y i ,k (y k -y i ) , 1 ≤ k ≤ K 1 . (15)
Fig. 2 shows a sample covariance matrix C y i with its eigenvectors and eigenvalues λ y i for a non-flat patch of the noisy parrot image in Fig. 1 where σ = 20. The eigenvectors are plotted in a raster order. The weights g(λ y i ) in ( 14) are depicted in the same figure. More than half of the eigenvalues are smaller than σ 2 , resulting in negative eigenvalues in C x i and negative filter coefficients g(λ). The resulting filter gives large negative weights to the eigenvectors associated to the smallest eigenvalues, which look like noise with no apparent structure. Remark 3. According to the Eckart-Young theorem, the capability of the eigenvectors of C y i to encode important features of the patch y i decreases as the eigenvalue decreases. As far as the negative weights g(λ k ) decrease, insignificant eigenvectors with inverted signs are amplified.

Aggregation

The process of individual patch MAP estimation terminates when all patches have been visited. A patch can be estimated multiple times using [START_REF] Dabov | BM3D image denoising with shape-adaptive principal component analysis[END_REF]. The set of all MAP estimates of a patch x k reads as X k := {i : k ∈ S 1 (i)}. All these patch estimates are aggregated by averaging them on their corresponding locations on the image. We denote by H the local indexes within a

√ K × √ K patch, H = {h = (h 1 , h 2 ) : 0 ≤ h 1 , h 2 ≤ √ K -1}.
A pixel j is contained in all the patches x j-h with h ∈ H. The value of the oracle image x at pixel j results from averaging all estimated values for that pixel:

x(j) = 1 A j h∈H i∈X j-h x (i) j-h (h), with A j = h∈H X j-h . ( 16
)
There is a large variability in the aggregation weights across the image. Some patches have a high degree of self-similarity with their surroundings, and appear in many nearest neighbors sets S 1 . Thus their aggregation weights A i can be quite higher than those of other patches with singular patterns. Several patch based methods have a similar aggregation of the patch estimates [START_REF] Dabov | Image denoising by sparse 3d transform-domain collaborative filtering[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: From Gaussian mixture models to structured sparsity[END_REF], and it was identified in [START_REF] Lebrun | Secrets of image denoising cuisine[END_REF] as one of the key ingredients of state-of-the-art image denoising strategies.

To evaluate the impact of the different stages of the denoising algorithm, we track the non-flat patches during the denoising of the parrot image with σ = 20. In Figure 3 we show the histograms of the RMSE of the non-flat patches in the input noisy image, (1) after the individual MAP estimates but before the aggregation, and (2) after their aggregation in the output oracle image. The RMSE before the aggregation (BA) is computed as

RMSE BA =   1 C i∈J k∈S1(i) x (i) k -x k 2 2   1 2 , (17)
where C = i∈J S 1 (i) for J the set of non-flat patches. The RMSE after the aggregation is computed similarly, with the patch of the oracle x k instead of x

k . The left figure corresponds to the first step of the NLB and the right to the second. It is seen that the MAP formula [START_REF] Dabov | BM3D image denoising with shape-adaptive principal component analysis[END_REF] in the first step gives quite a poor result, with a RMSE worse than that of the noisy input image. This is a consequence of the unreliable estimate of the covariance matrix. The aggregation stage has a crucial role to obtain the first-step denoised patch x i . Most of the denoising work is done during the first step of the algorithm; see also Fig. 1.

Fig. 4 shows different stages of the denoising of a patch y i : the MAP estimate (6), the aggregation stage in ( 16), followed by the second step MAP and the final denoised patch x i . Fig. 2 already showed that eigenvectors with small positive and also negative eigenvalues have a random structure. In Fig. 4, the significant contribution of such eigenvectors can be clearly seen: The step 1 MAP estimate before the aggregation is not better than the noisy patch. However, the aggregation of the MAP estimates yields a good approximation of the clean patch, very similar to the output after the second denoising iteration. Most of the denoising is carried out by the aggregation of the MAP estimates in the first stage. In Fig. 5 we show all the estimates that are aggregated for two pixels in the patch shown in Fig. 4. We refer to these plots as pixel paths. These values correspond to the MAP estimates (Eq. ( 6)) of patches that contain the pixels. It is quite a general rule that a pixel has a large number of updates. The variance of the updates is huge, which results from the negative weights in the pseudo-Wiener filter. The surprising fact is that in spite of this huge variance the final value is reasonable. By manipulating the values -removing outliers based on the histogram -one can get better results. Since the number of updates is not known, an automatic outlier rejection does not seem simple.

3.3

The effect of the thresholding the negative weights

In the light of the results after the first step MAP estimation (Figs. 3 and4) one is curious to better realize the role of the negative eigenvalues in the estimated matrix C x i . A common approach is to clip negative eigenvalues to zero [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF]. Our aim here is to go into the question using a family of estimates of the covariance matrix by discarding small eigenvalues in C y i -σ 2 I (see [START_REF] Arias | Video denoising via empirical Bayesian estimation of spacetime patches[END_REF]):

C τ x i = U y i diag T τ λ y i -σ 2 U T y i , (18)
where τ ∈ R and T τ (s) = s if s ≥ τ σ 2 and T τ (s) = 0 otherwise. We can then investigate the performance of the algorithm under variations of τ . In Fig. 6 we test the estimators C τ x i for different values of τ and evaluate the PSNR on non-flat patches after the individual MAP estimation (dashed lines), and after the aggregation (continuous lines) for both denoising steps. As expected the plots become constant below τ = -1, since all eigenvalues are larger than -σ 2 (C y i is positive definite). This is exactly the first-step MAP estimate [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF] in NLB [START_REF] Lebrun | A nonlocal Bayesian image denoising algorithm[END_REF]. When τ increases, more and more entries in the diagonal matrix in [START_REF] Zhang | Two-stage image denoising by principal component analysis with local pixel grouping[END_REF] are zeroed out, together with the corresponding filter coefficients g(λ y i ). For τ large enough C x i → 0 and the MAP denoising reduces to the x [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF].

(i) k = y i . If τ = 0 we recover estimator (C y i -σ 2 ) + of
The plot confirms the bad quality of the MAP estimates in the first denoising iteration when τ ≤ -1, and the impressive correction caused by the aggregation. Regarding the output of the individual MAP estimates in the first step, the best results are obtained for τ = 0 (1.8dB higher than with τ ≤ -1). However, this is not the case after the aggregation (continuous red curve), where slightly better results (around 0.1dB higher) correspond to τ ≈ 2, and also to τ ≈ -0.7. The PSNRs for the second step of the algorithm are almost independent of τ . The negative eigenvalues in the first step do not have a negative effect on the performance of the whole algorithm. In fact, allowing some small negative eigenvalues seems slightly better than keeping all the positive eigenvalues. Although we show our analysis on a single example, we have found that this is the generally the case.

Second step denoising

The second step denoising is essentially a Wiener filter where the a priori statistics are estimated using the oracle image resulting from the first step. In addition, the oracle image is also used to define the set of similar patches S 2 (i). As in the first step, following Remark 1, all patches in S 2 (i) are denoised as proposed in (3). The a priori covariance matrix is estimated as the sample covariance matrix of the patches from the oracle image C x i , which are assumed noiseless. The coefficients of the resulting Wiener filter matrix on the basis of eigenvectors U x i have the well-known expression f (λ) = λ/(λ + σ 2 ). Fig. 7 shows a sample covariance matrix C x i with its eigenvectors and eigenvalues, and the weights f (λ x i ). The covariance matrix corresponds to the reference patch in Fig. 2. The filter weights concentrate their mass on the largest eigenvalues, which correspond to smooth eigenvalues. This is contrary to what happens in the first step, as shown in Fig. 2. Thus, unlike to the the first step, eigenvectors with a random structure are much less present in the second step MAP estimates. The final denoised image is obtained by aggregation as described in subsection 3.2. Although the second step has a contribution to the quality of the final result, arguably most of the denoising is done in the first step, as can be seen by the Figures 1, 3 and 4.

Discussion

Several patch-based methods for image denoising follow Bayesian approaches with a Gaussian (or GMM) prior for patches. The validity of the Gaussian model is questionable, but these methods produce state-of-the art results and are appealing due to their simplicity. In the case of NLB, the indefinite prior covariance estimate causes poor individual patch MAP estimates. However, they systematically lead to sets of pixel estimates whose mean (aggregate value) is visibly consistent. A good understanding of the NLB algorithm requires models that include the aggregation stage.

Fig. 1 .

 1 Fig. 1. From left to right: noisy image y for σ = 20; first-step denoised x (29.57 PSNR), denoised x (29.95 PSNR) and original image.

  (a) The least eigenvalue of Cy i for all reference patches in non-flat areas. (b) Percentage of reference patches y i in non-flat areas with at least one negative weight in g(λy i ). (c) Average number of negative weights in g(λy i ) for reference patches y i in non-flat areas.

Fig. 2 .

 2 Fig.2. For a non-flat patch in the noisy parrot image in Fig.1, with σ = 20, from left to right: sample covariance matrix Cy i ; its eigenvectors in raster order; the eigenvalues λy i and the resulting weighting function g(λy i ). Notice the large filter weights g(λ) corresponding to λ < σ 2 .

Fig. 3 .Fig. 4 .

 34 Fig. 3. Histograms of the RMSE for non-flat patches, for the parrot image with σ = 20.

Fig. 5 .

 5 Fig.5. Pixel paths for two pixels in the patch y i shown in Fig.4. These are the values that are aggregated according to[START_REF] Yaroslavsky | Theoretical Foundations of Digital Imaging Using MATLAB[END_REF].

  before agg. step 1, after agg. step 2, before agg. step 2, after agg.

Fig. 6 .

 6 Fig.6. PSNR of all non-flat patches after different stages of the algorithm, using the thresholded covariance matrices C τx i , as a function of τ . See text for details.

Fig. 7 .

 7 Fig.7. For the same non-flat patch as in Fig.2we show, from left to right: sample covariance matrix C x i ; its eigenvectors in raster order; the eigenvalues λ x i and the resulting weighting function f (λ x i ). Compare to the weights in Fig.2.

Table 1 .

 1

		1	0.087	10.25 %	1.702 %
	(512 × 512)	σ = 5	0.883	96.23 %	40.75 %
		σ = 10	3.747	99.85 %	53.47 %
		σ = 20	31.88	100 %	55.13 %
		σ = 40	122.7	100 %	62.25 %
	man	σ = 1	0.059	10.37 %	2.584 %
	(1024 × 1024) σ = 5	0.909	95.84 %	41.41 %
		σ = 10	3.833	99.64 %	54.65 %
		σ = 20	29.01	100 %	56.83 %
		σ = 40	134.6	100 %	62.23 %
	parrot	σ = 1	0.051	45.43 %	16.59 %
	(256 × 256)	σ = 5	1.201	80.75 %	29.84 %
		σ = 10	4.811	93.88 %	42.26 %
		σ = 20	35.22	100 %	48.16 %
		σ = 40	134.9	100 %	89.42 %
	peppers	σ = 1	0.126	13.48 %	1.914 %
	(512 × 512)	σ = 5	0.772	97.67 %	46.82 %
		σ = 10	3.309	99.88 %	56.94 %
		σ = 20	34.05	100 %	56.63 %
		σ = 40	127.3	100 %	62.62 %
	stream	σ = 1	0.093	5.194 %	0.856 %
	(512 × 512)	σ = 5	0.942	67.77 %	20.44 %
		σ = 10	4.354	97.14 %	39.19 %
		σ = 20	29.55	100 %	50.38 %
		σ = 40	130.5	100 %	61.69 %


Atis Elsts

Xenofon F Fafoutis

Adeyinka Adeleke

Robert Piechocki

George Oikonomou

Simon Duquennoy

Antonio Liñán

Marc F Abregas

Antonio Li

Ñ Án

Marc F Àbregas Zolertia

Competition: Adaptive Time

Keywords:

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

Background

There is a growing need to make low-power wireless networks more reliable and more predictable in order to open them up to a wider range of applications (e.g., industrial, automotive, e-health applications). Time-Slotted Channel Hopping (TSCH), specified in the IEEE 802.15.4-2015 standard [1], offers these features and has recently attracted attention from both industry and academia.

TSCH makes use of pseudorandom channel hopping to combat external interference and frequency-selective multipath fading. The core observation here is that several major causes of bad performance, such as WiFi interference or deep fading are very unlikely to equally affect all of IEEE 802.15.4 MAC-layer channels. Since a retransmission of an unacknowledged packet in TSCH with high probability takes place on a different channel, having just a few lowquality channels do not render the whole network unusable. In a similar way, changing the channel allows the nodes to get out of frequency-specific deep fading pockets caused by the multipath effect. In this way, the reliability of the net-working plane in increased by making use of frequency diversity. Since the low-power wireless environment is notoriously dynamic and unpredictable, TSCH typically shows better results than adaptive selection of a single channel, as seen, for example, in ZigBee.

OpenWSN [START_REF] Watteyne | OpenWSN: a standards-based low-power wireless development environment[END_REF] is the de facto TSCH open-source reference implementation. In the recent years, Contiki port of TSCH has emerged as a competitor [START_REF] Duquennoy | Orchestra: Robust Mesh Networks Through Autonomously Scheduled TSCH[END_REF]. OpenWSN and Contiki versions of TSCH were demonstrated to be interoperable in Internet Engineering Task Force (IETF) plug-tests in Prague (2015) and in Paris (2016). Similarly to Open-WSN, the Contiki implementation supports multiple hardware platforms: NXP JN516x, Texas Instruments CC2538 and CC2650, Zolertia Z1, and others, including TelosB.

The Contiki implementation of TSCH is characterized by energy-efficient wake up pattern. Before entering sleep mode, the microcontroller calculates the time to the next active slot (transmission or reception), and schedules a wake up. In this way, no energy is wasted waking up at the start of every idle slot.

Our enhancements to TSCH

The objectives of this work are to increase energy efficiency, decrease latency, and keep high reliability of TSCH networks even in harsh conditions. To achieve these objectives, we will experiment with a number of techniques.

Adaptive time synchronization

This technique consists of learning the drift of the local clock compared to the drift of the network's time source and adaptively compensating for this expected drift [START_REF] Stanislowski | Adaptive synchronization in IEEE802.15.4e networks. Industrial Informatics[END_REF]. In this way, the desynchronization between nodes can be kept small; in particular, it is bounded as long as the temperature changes in the environment remain bounded as well.

Adaptive guard time

Idle listening in the Rx guard time is know to have high impact on TSCH energy consumption [START_REF] Elsts | Microsecond-Accuracy Time Synchronization Using the IEEE 802.15.4 TSCH Protocol[END_REF]. By increasing the accuracy of time synchronization, this idle listening can be reduced. We plan to dynamically reduce guard time depending on the state of time synchronization. Our results show that decreasing guard time leads to higher energy efficient compared to the more conventional method of decreasing the rate of network synchronization packets [START_REF] Elsts | Microsecond-Accuracy Time Synchronization Using the IEEE 802.15.4 TSCH Protocol[END_REF].

Adaptive channel selection

To avoid transmissions on low-quality IEEE 802.15.4 channels, either the hopping sequence itself may be changed on the network nodes, or the nodes can pseudorandomly remap the set of low-quality channels to a different set of channels every time a transmission is done (Fig. 1). We have experimented with multiple techniques [START_REF] Elsts | Adaptive Channel Selection in IEEE[END_REF]: upstream-based selection and downstream-based selection, as well as multiple channel quality metrics: PRR and idle-channel energy levels. Assuming a situation where some of the channels are more interfered than other channels, the adaptive selection helps to decrease link-layer packet error rate (Fig. 2). Since the dependability competition requires relatively fast adaptations compared to the number of packets being sent, upstream-based channel selection using the idle-channel energy level metric is expected to perform better than the downstream-based selection using the PRR metric.

Network join process improvements

The robustness of TSCH join process can be improved by increasing the number of physical channels, as this increases the probability that some of the less-interfered channels will be used for sending enhanced beacon (EB) packets. However, using more channels also makes the join process slower in the average case. Similarly, the speed of the process can be improved by sending EB packets more frequently; however, too fast EB scheduling is going to lead to packet collisions and actually decrease performance. We will investigate practical solutions for these optimization problems, i.e., selecting the number of channels to use and the EB frequency. We also will experiment with sending EB packets on multiple channels simultaneously, and sending EB packets in multiple slots in the same slotframe.

Routing and scheduling

We plan to use flooding or multi-path forwarding to additionally exploit the spatial diversity of nodes besides the frequency diversity of channels already naturally provided by TSCH. Link selection will be based on EB packet reception probabilities and other factors. The design of an efficient TSCH schedule for this competition is an open question that we plan to tackle at a later stage.

Improvements in hardware-specific code

The Contiki hardware driver for the CC2420 radio chip can be improved in terms of energy efficiency, for example, by turning the chip off completely during idle periods -such as is done, for example, in the ManOS operating system [START_REF] Strazdins | MansOS: easy to use, portable and resource efficient operating system for networked embedded devices[END_REF], rather that keeping it in idle mode, with its voltage regulator remaning powered. This "radio off" mode compared to the "radio idle" mode has the drawback of a much slower wakeup. However, since the time of the next active TSCH slot is known beforehand, techniques such as pre-emptive wakeup can be used to work around this problem. These techniques are already present in the code of more recent Contiki hardware platforms, for example, CC2650 [START_REF] Elsts | Microsecond-Accuracy Time Synchronization Using the IEEE 802.15.4 TSCH Protocol[END_REF] and JN516x; we will backport them to CC2420.

Similarly, MCU energy usage can be optimized by putting the MCU in sleep during radio packet reception and transmission. This requires changes both in the CC2420 driver and the core TSCH code.

Concluding remarks

We plan to build on the open-source implementation of TSCH in Contiki already known to provide high reliability in terms of packet delivery rate even in challenging environments [START_REF] Duquennoy | Orchestra: Robust Mesh Networks Through Autonomously Scheduled TSCH[END_REF]. We plan to enhance it in several ways to make the implementation more energy-efficient and have lower latency.

Figure 1 :

 1 Figure 1: Overview of the adaptive channel selection process [4].

Figure 2 :

 2 Figure 2: The effect from adaptive channel selection on link-layer packet error rate. The adaptive algorithm is able to select and use the less interfered channels (with 20 % PER), keeping the overall packet error rate low.

Acknowledgments

This work was performed under the SPHERE (a Sensor Platform for Healthcare in a Residential Environment) IRC funded by the UK EPSRC, Grant EP/K031910/1 and a grant from CPER Nord-Pas-de-Calais/FEDER DATA.


S Lewandowski
email: simon.lewandowski@onera.fr.

S Duzelier
email: sophie.duzelier@onera.fr.

J Eck

S Dagras

C Tonon

H Jochem

P Jouanne
email: pierre.jouanne@thalesaleniaspace.com.

E Laurent

J M Desmarres

Particle Flux Effects on Physicochemical Polymer Degradations

Particle flux effects on polymer physicochemical degradations were investigated. Different polymers were studied depending on application: silicone, fluoropolymer, acrylic, epoxy, polyimide, polyurethane, etc. Polymers with surface functional properties (paint, coating, etc.) were dissociated from polymers with bulk properties (potting, cable, etc.). Surface materials were exposed to protons, and bulk materials were exposed to electron radiations. Two fluxes were used: 5 × 10 11 and 7.4 × 10 12 cm -2 • s -1 for electrons and 8.7 × 10 9 and 2.2 × 10 11 cm -2 • s -1 for protons. The total dose was fixed at 500 Mrad for electron radiations and 1 Grad for proton radiations. Functional properties (thermo-optical and mechanical) depending on material application were measured before and after irradiation exposure. Chemical structure modifications were also investigated. From all combined results, no global tendency was evidenced. The study shows that low or high irradiation flux induces more or less degradations according to polymer type, the composition of the polymer, and functional properties of concern. Flux effects can also depend on particle nature (electron or proton). The last point was highlighted by silicone polymer degradations; high proton flux caused more damage than low flux. On the contrary, low electron flux induced more degradations than high flux.

Nomenclature T g = glass-transition temperature, °C α s = solar absorptance ε irr = infrared emissivity

I. Introduction

I N SPACE applications, materials are exposed to specific environ- ments: charged particles (from the radiation belts), the solar spectrum, atomic oxygen (in low orbit, <1000 km), and micrometeorites and debris. Vacuum conditions and temperature are added to this. Among these different environmental parameters, the charged particles are constituted by electrons or protons trapped in the radiation belts, rash protons (also x rays on the sunny side) and solar wind. Materials on external parts of the satellite are directly exposed to this harsh environment. Charged particles mainly generate an ionizing absorbed dose profile into these materials. Because of their energy spectrum and the mode of interaction, the protons are responsible for the surface dose (first micrometers) and electrons for the bulk dose.

The interaction of particle radiation with an organic material causes the phenomena of chain scission and/or crosslinking via ionization mechanisms. The effects are linked to the accumulation over time of the ionizing dose (cumulative effect) causing surface and/or volume changes in the material. The induced macroscopic effects consist mainly of color change (absorbance increase and transmittance decrease) and the decrease of mechanical properties performance (mass loss, brittleness, reduced stress and strain at the break, etc.).

A classification of polymer materials by particle sensitivity and degradation type is given in the literature survey depending on polymer structures (e.g., the presence of the aromatic group leads to better stability) and the chemical composition of its monomer [START_REF] Kinderman | Nuclear Radiation Effects on Structural Plastics and Adhesives[END_REF]. Additives can also have a high influence on degradation. An index radiation indicator is given by Tavler et al. [European Laboratory for Particle Physics (CERN)] for prioritizing materials per family according to their mechanical behavior (ultimate strength and deformation) [START_REF] Tavlet | Compilation of Radiation Damage Test Data Part II[END_REF]. However, these classifications are based on data obtained under specific conditions, generally not representative of the space environment (data from nuclear industry or from neutrons and a gamma accelerator, with exposure to air). They only could be used as indicative.

Therefore, evaluation of material degradation during space missions is estimated by performing ground tests that closer simulate the space environment in order to select and qualify the appropriate materials [START_REF] Tighe | Space Environment Simulation and Materials Testing[END_REF]. For this, test facilities are used to replicate specific degrading components of the space environment. For many reasons, particularly to reduce the aging time and due to specifications of the particle accelerator, these tests are not performed under the exact same conditions as the environment the material will see over the lifetime of the mission. Accelerated tests, especially important for particle irradiation, are performed, according to specific standards [START_REF]Space Systems, Space Environment, Simulation Guidelines for Radiation Exposure of Non-Metallic Materials[END_REF][START_REF]Standard Practice for Combined, Simulated Space Environment Testing of Thermal Control Materials with Electromagnetic and Particulate Radiation[END_REF][START_REF]Space Product Assurance Particle and UV Radiation Testing for Space Materials[END_REF]. The possible acceleration impact on the material degradation is currently not evaluated. Dose rate effects have been observed in the literature, with potentially more degradation observed for low dose rates [START_REF] Tavlet | Behaviour of Organic Materials in Radiation Environment[END_REF][START_REF] Gillen | Time-Temperature-Dose Rate Superposition: A Methodology for Extrapolating Accelerated Radiation Aging Data to Low Dose Rate Conditions[END_REF][START_REF] Wilski | The Radiation Induced Degradation of Polymers[END_REF]. These dose rate effects may be due to the oxygen environment and the oxidation of radiolysis species created during irradiation. Again, many data come from nuclear applications (neutrons/in air) or CERN (γ and neutrons/in air). In this case, because of the oxygen involvement, surface effects are dominant, but low particle flux tests (for long durations) on thin layers can be critical to bulk properties. We can minimize the effect of oxygen by working under a stream of neutral gas (nitrogen and argon) or under vacuum, provided that there is no adsorbed oxygen. The effect of dose rates can be seen in vacuum particle irradiation. Briskman et al. report data on Polytetrafluoroethylene (PTFE) and low-density polyethylene [START_REF] Briskman | Determination of Dose Rate Effects in Polymers Irradiated in Vacuum[END_REF], including the gas emission rate measurement, for which significant differences were observed for different dose rates. An increase of degradations of tensile strength and elongation at rupture was observed for PTFE at low dose rates. In contrast, Hardgrove presents other electrons e -protons p dose data at different acceleration rates that show for white paints that the worst case is observed at high acceleration rates [START_REF] Hardgrove | Space Simulation Test for Thermal Control Materials[END_REF]. A comparison with in-flight data leads him to the conclusion that the dose rate closest to space arrival provides the most realistic degradation value.

Considering these different results, the study presented in this paper aims to evaluate the particles flux effects on the physicochemical properties of different polymers. Indeed, as previously indicated, the irradiation aging tests performed on the ground in order to simulate the materials aging in orbit cannot be fully representative of the in-flight environment. In particular, the particle flux on the ground is highly accelerated compared to that in the space environment. In the first part of this paper, the experimental process is described with a presentation of the materials tested, the different irradiation conditions, and the associated analysis methods.

In the second part, a global synthesis on the flux effects observed on the different polymers tested is presented.

Finally, in the last part, a focus on the flux effects on silicone rubber is discussed.

II. Materials and Experimental Methods

A. Materials

Various polymers used on satellites were studied including surface materials, such as paints and coatings, and bulk materials, such as cables, fixing, and structural bonding.

The following surface materials that were of concern that were tested include room temperature vulcanized (RTV) silicone, polyurethane (PU), polyetherimide (PEI), and polyester film.

The investigated bulk materials were acrylic adhesive polyimide tape, RTV silicone, epoxy, fluorinated ethylene propylene (FEP), and acetal.

Two types of irradiations tests were performed depending on the function of the materials studied; surface materials were exposed to proton flux, and bulk materials were exposed to electron flux.

B. Irradiations 1. Proton Flux Irradiation

Proton irradiation was performed in the SEMIRAMIS facility (ONERA), using a 2.5 MeV Van de Graff accelerator. The proton energy was 240 keV, and the total dose was 1 × 10 7 Gy (penetration at about 3 μm). Two different fluxes were selected: 8.7 × 10 9 p • cm -2 • s -1 and 2.2 × 10 11 p • cm -2 • s -1 . A vacuum lower than 3 × 10 -7 mbar was obtained after a 1 day pump down period. The sample holder temperature was maintained at 25°C for the duration of the test.

Electron Flux Irradiation

The electron irradiation was performed by TRAD using a Van de Graff accelerator able to produce high-energy electrons (1-3 MeV).

A fixed electron current is delivered by the horn. The samples are horizontally disposed on the table. The dose rate is directly correlated to the average amount of electrons received by the samples under the beam.

The irradiation table is cooled down to 20°C. A thermocouple is put on the table next to the samples, to check with a remote controller the eventual temperature increase during the experiment.

A display screen in the control room provides a direct view of the table and the samples.

The irradiation is performed under a nitrogen environment, thanks to a flowing nitrogen setup containing all the samples that are irradiated.

Before every irradiation, a preliminary dosimetry step is performed in order to determine the dose rate.

The dosimetry is performed using thin film radiachromic dosimeters at a given beam current and by measuring the deposited dose. The dose rate is calculated as the needed beam current to perform the irradiation at the target dose rate.

In the frame of this study, irradiations are performed at 3 MeV and at standard and high dose rates, respectively, 5 × 10 11 and 7.4 × 10 12 e -• cm -2 • s -1 .

C. Analysis Techniques 1. Ultraviolet Visible Near-Infrared Spectroscopy

Different Ultraviolet-Visible-Near-Infrared (UV-Vis-NIR) analyses were performed at ONERA: 1) For measurements in air, the UV VIS NIR measurement system was composed of a Perkin Elmer Lambda 1050 spectrophotometer associated with a 150 mm diameter integrating sphere in air. The calibration was performed using a diffuser Spectralon SRS-99-010 standard. The diffuse transmission spectra in the range 200-2500 nm were acquired at intervals of 1 nm, with the lamp change at 319.2 nm.

2) For measurements under vacuum, the SEMIRAMIS facility is equipped with a Perkin Elmer Lambda 1050 spectrophotometer connected to a mobile integrating sphere under vacuum. Reflectance spectra in the range 250-2500 nm can be obtained, without breaking vacuum, before and after irradiation.

X-Ray Photoelectron Spectroscopy

Atomic concentration and chemical bonding information were obtained from the specimen surfaces by x-ray photoelectron spectroscopy (XPS). XPS measurements were realized with a Thermo Scientific K alpha spectrometer with a monochromated Al Kα radiation (hν 1486.6 V). Spectra were collected using an x-ray spot size of 400 μm.

Attenuated Total Reflection Fourier Transform Infrared Spectroscopy

Attenuated total reflection (ATR) Fourier transform infrared (FTIR) spectroscopy was performed at ONERA to characterize surface chemical degradations. The ATR FTIR spectra were recorded using a Nicolet 800 FT-IR spectrometer equipped with a deuterated triglyceride sulfate detector in the range 400 -4000 cm -1 , with a resolution of 4 cm -1 . A diamond ATR accessory was used to investigate the first micrometers of the samples.

Nuclear Magnetic Resonance Spectroscopy

29 Si solid-state nuclear magnetic resonance spectra were recorded on a Brucker Avance 400 spectrometer equipped with a 4 mm probe operating at 9.39 MHz. Samples were spun at 12 kHz at the magic angle using ZrO 2 rotors. 29 Si Magic Angle Spinning (MAS) spectra were obtained under high-power proton decoupling conditions, with a small flip angle (30 deg) and a recycle delay of 60 s.

Differential Scanning Calorimetry

At Thales Alenia Space, Differential Scanning Calorimetry (DSC) measurements were carried out using a NETZSCH DSC 204 apparatus for RTV silicone and a DSC823e from Mettler Toledo for other polymers. The ramp was fixed at 10°C∕ min. The analysis temperature range was fixed depending on polymer characteristics.

Peel Test

Peeling tests were performed in accordance with ASTM D1000 standard at Airbus Industries. Peelings were performed with a 180 deg angle at 400 mm∕ min and at ambient temperature on a Zwick Roell facility.

Adherence Test

Adherence tests performed on coatings according to ISO EN 2409. P224 tape was used at Airbus Industries.

Shear Test

Shear tests were performed according ISO specification (ISO.4587:2003). Tests were performed at room temperature on a Zwick Roell facility at Airbus Industries.

Shear Tests were performed on adhesives (epoxy and silicone) at Thales Alenia Space. Tests were performed at room temperature on a Zwick 100 kN apparatus.

Nanoindentation

At Centre National d'Études Spatiales (CNES), nanoindentation was carried out with a nanoindenter AX (MTS), using a Berkovich point and a depth penetration of 5 μm. Twelve indentation points were performed on each sample. A sample Poisson ratio was fixed at 0.4.

Dynamic Mechanical Analysis

Three-point flexure tests were performed on a Dynamic Mechanical Analysis (DMA) 242 NETZSCH at CNES. The temperature range of measurement was fixed according to polymer characteristics. Measurement was performed with a ramp of 5°C∕ min and a frequency of 1 Hz. The distortion amplitude was 30 μm, dynamic strength limit was fixed at 7.2 N, and the constant part of static strength was fixed at 0 N. Each sample was tested two or three times.

III. Flux Effects on Polymers: General Results

A. Surface Irradiation

Table 1 summarizes all the results obtained on surface materials during the aging tests.

No flux impact was evidenced except for RTV silicone. RTV silicone showed higher degradation at high proton flux.

B. Bulk Irradiation

Table 2 summarizes all the results obtained on bulk materials during the aging tests.

At this dose level and with the considered fluxes, no flux impact was evidenced except for RTV silicone, the second epoxy, and FEP. It has to be noted that in the same material family the degradation could be different (epoxy adhesive). Fig. 2 Reflectance from 250 to 800 nm of pristine and protons irradiated silicone rubber.

IV. Flux Effects on Silicone Rubbers

To illustrate the process used and the results obtained in this study, particular attention is paid to the material considered as a surface and bulk sample during this testing: RTV silicone rubbers.

A. Surface Irradiation

A slight yellowing of samples was observed after proton irradiations. Moreover, inhomogeneous degradation on the high flux sample was noted, as shown on the optical micrograph (Fig. 1). An important surface damage was observed with the formation of many cracks localized in one-half of the sample. This phenomenon was previously observed on silicone under proton irradiation [START_REF] Jochem | Effects of 400 keV Electrons Flux on Two Space Grade Silicone Rubbers[END_REF][START_REF] Jochem | Effects of Electron and Proton Flux on the Chemical Structure of a Silicone Rubber[END_REF].

Thermo-Optical Properties

The SEMIRAMIS facility enables measurement of thermo-optical properties in situ (under vacuum). Thermo-optical properties were measured by reflectance under air and vacuum conditions. Results are presented for the two different irradiation fluxes in Fig. 2.

A UV cutoff shift to higher wavelengths was observed under vacuum after irradiation, showing absorbent species formation. This shift is more important after the high flux irradiation. After the replacement of the samples under ambient air, new measurements showed that there were recoveries of their reflectance spectra. This recovery had been previously observed on different materials [START_REF] Kiefer | Space Environment Effects on Polymeric Materials[END_REF][START_REF] Marco | Eight Years GEO Ground Testing of Thermal Control Coatings[END_REF]. This phenomenon is particularly significant for the high flux irradiated sample. The low and high flux final spectra (measurement in air) were closed.

Thermo-optical coefficients (IR emissivity and solar absorptances) have been calculated. IR emissivity changes were not significant. Solar absorptances are reported in Table 3.

Recovery phenomenon was also observed in the change in α s , particularly for the high-flux irradiated sample. Accordingly, the α s end of test values (measurement in the air environment) for both aged samples were similar, even though the interim α s (measurement under vacuum just after irradiation) was higher when proton flux was higher.

X-Ray Photoelectron Spectroscopy

XPS measurements show that the atomic percentages obtained for the different samples are similar (C, O, and Si). However, a shoulder corresponding to the SIO x>3 bond is observed for irradiated samples. This is more evident for the high-flux sample (Fig. 3). This result indicates an evolution of the surface nature of the silicone rubbers into silica, which is a function of the total dose and the irradiation flux.

Attenuation Total Reflection Fourier Transform Infrared Spectroscopy

ATR FTIR spectra obtained before and after proton irradiation are presented in Fig. 4.

A distinction can be seen on the ATR FTIR measurements between the high-flux cracked sample and the other irradiated samples (lowand high-flux noncracked):

1) For the high-flux cracked sample, ATR FTIR measurement confirms silica formation on the high-flux irradiated sample, with an increase of the SiOH bands absorption (3100 -3600 cm -1 ) [START_REF] Zhang | Damage Effects and Mechanisms of Proton Irradiation on Methyl Silicone Rubber[END_REF][START_REF] Tomozawa | Infrared (IR) Investigation of the Structural Changes of Silica Glasses with Fictive Temperature[END_REF]. The appearance of a new absorption band (1700 cm -1 ), attributed to the CO groups, responsible for the sample yellowing, is also observed.

2) For the other irradiated samples, spectra obtained after irradiation are closer to that of the pristine sample. However, an increase of the SiH and a decrease of the SiCH3 and SiPhenyl groups are observed. A SiOSi characteristic peak (1070 cm -1 ) decreases after irradiation, for both high-flux irradiated samples, while it is unchanged for the low-flux irradiated sample.

The formation of a silica-type layer localized on the surface can explain why cracks appear: the density of silica is higher than that of silicone rubber.

29 Si Solid-State Nuclear Magnetic Resonance Spectroscopy

29 Si solid-state nuclear magnetic resonance (NMR) spectra of all silicone rubber samples are presented in Fig. 5 (e -and p irradiation).

Pristine and irradiated spectra show a high signal at -22 ppm attributed to D siloxane units within the chain. Considering proton irradiation, a slight peak corresponding to SiCH 2 Si bonds appears Fig. 3 XPS spectra of silicone rubber before and after proton irradiation.

B. Bulk Irradiation

To assess the bulk properties of the material, different characterization techniques were used such as UV VIS NIR spectroscopy, NMR, DMA, nanoindentation, and DSC to determine the evolution of the physicochemical properties of the material submitted to the electron beam.

The ATR FTIR spectroscopy shows a significant degradation of the samples submitted to a standard dose rate, leading to the apparition of SiOH, CO, and SiH bonds. For high dose rate, the sample modification consists mainly in the creation of SiH bonds and on the decrease of the number of SiPhenyl bonds (Fig. 6). This tendency is confirmed by the NMR analysis, which reveals the formation of crossing points (SiO 3 bonds) coupled with the diminution of the SiPhenyl bonds (Fig. 5). After irradiation at the standard dose rate only, another type of crossing point is also observed (SiCH 2 Si).

Complementary NMR analysis shows a stiffening of the silicone skeleton after electron irradiation.

The DMA tests focused on the mechanical behavior of the material depending on the temperature (Fig. 7). Results show an evolution of the temperature of the inflexion point for the storage modulus, with a decrease of about 10°C after electron irradiation. Additionally, no significant difference was noticed between the values after standard and high dose rate irradiation.

DSC analysis was also performed. The results reveal that the material was not fully polymerized, leading to the apparition of an exothermic peak during the first scan. These peaks were localized at -118 and -100°C for pristine and high dose rate irradiated samples, respectively. It can be noted that such a peak was not observed for the sample irradiated at the standard dose rate. The evolution of the T g was also monitored. Its value was measured at -110, -103, and -104°C, respectively, for one pristine sample, the one irradiated at the standard dose rate and the one submitted to a high dose rate. Thus, there was an evolution of this temperature after irradiation, but no effect of the dose rate was noticed.

Finally, mechanical properties of the material were assessed through the irradiation of shear test samples. A clear decrease of the tensile strength of around 40% was measured, from 4.1 to 2.3 MPa. The samples irradiated at the standard dose rate seemed slightly more degraded.

V. Conclusions

Different materials families were tested in the frame of this study. They were selected to represent some of the main organic material families used on spacecraft such as epoxies, silicones, polyurethanes, fluoropolymers, etc.

The aim of this study was to determine the influence of the irradiation flux on material degradation. It has to be noticed that, due to facility limitations, the flux as applied in the frame of this study, even the lowest ones, was still largely above real fluxes in orbit.

During this study, various analysis tests were performed on the materials in order to be able to determine the evolution of the functional properties of the tested materials but also to determine some physicochemical changes.

The general following observations can be reported: 1) There is no evidence of a similar effect of the irradiation flux on the different materials families. Indeed, for some materials families, the degradation is higher with a higher irradiation flux, while it could be the contrary for another material family.

2) The effects can also be different for two materials in the same family (different behavior observed on two epoxies).

3) The irradiation flux effect can be different on the same material depending on the type of particle irradiation (electrons or protons).

This study demonstrated that ground-testing parameters are important to the evaluation of material behavior against radiation. The current study focused on the irradiation flux effect. The investigation of other irradiation parameters such as the type of particles, the atmosphere, and the irradiation temperature is also relevant and will be investigated in future work.

1 ) 1 )Fig. 1

 111 Fig. 1 Optical micrograph of the high-flux proton irradiated sample (cracks on the left side).

Fig. 4 ATR

 4 Fig. 4 ATR FTIR silicone rubber before and after proton irradiation.

Fig. 5

 5 Fig.529 Si NMR spectra of silicone rubber before and after electron and proton irradiation.

Fig. 6 ATR

 6 Fig.6ATR FTIR silicone rubber before and after electron irradiation.

Table 1

 1 Flux effects on surface polymers10 9 part • cm -2 • s -1 2 × 10 11 part • cm -2 • s -1 10 11 part • cm -2 • s -1 7 × 10 12 part • cm -2 • s -1

				Table 2 Flux effects on bulk polymers
		1 Grad (p )	1 Grad (p )		500 Mrad (e -)	500 Mrad (e -)
	Materials 6 × RTV silicone 1) Increase of α s , decrease of stress at break and	Materials 5 × Acrylic tape 1) High adherence degradation due to outgassing
		modification of the polymer structure			under irradiation for high and low flux: bubble
		2) Higher degradation for High flux			creation (adherence test: -95% after irradiation)
		3) (2 × 10 11 part • cm -2 • s -1 ) (developed in Sec. IV.A)		2) No significant flux effect at this dose level
	PU	1) No significant degradation (variation of α s , ε irr , and	RTV silicone	1) Degradation of mechanical properties
		adherence test lower than the uncertainty of the		2) Increase of the thermo-optical properties:
		measurement)			yellowing (visual inspection)
		2) No flux significant impact			3) Higher degradation for low flux (developed in
	PEI	1) No significant degradation (variation of α s , ε irr lower		Sec. IV.B)
		than the uncertainty of the measurement)	1st epoxy adhesive	1) Low modification of mechanical properties
		2) No flux significant impact		bicomponent, glass beads: 1 wt%	(stress at break: -15% after irradiation) even if physicochemical property changes observed [T g
					decrease (-10°C after irradiation) and observation of
					a new pic in ATR]
					2) No significant flux impact
				2nd Epoxy adhesive
				bicomponent,
				glass beads: 12 wt%

Table 3

 3 Evolution of solar absorptance α s at each step of proton irradiation

			Pristine	Irradiated 1 Grad
			Measurement
	Materials	In air	Under vacuum	Under vacuum	In air
	Low flux	0.19	0.18	0.20	0.21
	High flux	0.17	0.17	0.27	0.21

Article in Advance / LEWANDOWSKI ET AL. Downloaded by UNIVERSITA DEGLI STUDI DI MILANO on December 22, 2016 | http://arc.aiaa.org | DOI: 10.2514/1.A33500

Article in Advance / LEWANDOWSKI ET AL.

Downloaded by UNIVERSITA DEGLI STUDI DI MILANO on December 22, 2016 | http://arc.aiaa.org | DOI: 10.2514/1.A33500


Rémy Mingant
email: remy.mingant@ifpen.fr

Julien Bernard

Valérie Sauvant-Moynot

Novel state-of-health diagnostic method for Li-ion battery in service

Introduction

The development of improved State-of-Health (SoH) diagnostic methods is a current research topic for batterypowered applications [START_REF] Ozkurt | Integration of sampling based battery state of health estimation method in electric vehicles[END_REF][START_REF] Omar | Lithium iron phosphate based battery -Assessment of the aging parameters and development of cycle life model[END_REF][START_REF] Hu | Method for estimating capacity and predicting remaining useful life of lithium-ion battery[END_REF]. Due to the rapid development and deployment of Electric Vehicles (EV) there is a strong demand for an accurate and reliable technique for monitoring SoH during operation. Such a technique requires a measurable indicator that can be processed in an embedded system to give a quantitative measure of SoH. This is a key parameter required to optimize battery energy management and to track the degradation of the system performance. In this study, we investigate a novel "quasi electrochemical impedance spectroscopy" (QEIS) technique and apply it to determine SoH of a battery aged with a current profile which mimics the electrical needs of an electric school bus.

Electrochemical impedance spectroscopy (EIS) is a well-known and powerful tool used to analyze [START_REF] Dolle | Development of reliable three-electrode impedance measurements in plastic Li-ion batteries[END_REF][START_REF] Zhai | The Application of the EIS in Li-ion Batteries Measurement[END_REF][START_REF] Andre | Characterization of high-power lithium-ion batteries by electrochemical impedance spectroscopy. I. Experimental investigation[END_REF] and model [START_REF] Buller | Modeling the dynamic behavior of supercapacitors using impedance spectroscopy[END_REF][START_REF] Buller | Impedance-based non-linear dynamic battery modeling for automotive applications[END_REF][START_REF] Montaru | Frequency and Temporal Identification of a Li-ion Polymer Battery Model using Fractional Impedance[END_REF] electrochemical systems such as batteries. It is a non-invasive technique relying on a small perturbation signal being applied over a wide frequency range, and allows the investigation of internal phenomena with different time constants in a single measurement.

The general principle of the EIS method is to apply to a cell a sinusoidal signal and measure the characteristic response from the cell [START_REF] Andre | Characterization of high-power lithium-ion batteries by electrochemical impedance spectroscopy. I. Experimental investigation[END_REF]. The applied signal can either be current or voltage. The response 𝑢𝑢 �(𝑓𝑓) (or 𝚤𝚤 �(𝑓𝑓))

(complex numbers) are recorded. The impedance 𝑍𝑍 ̅ (𝑓𝑓) is then calculated by using the following equation:

𝑍𝑍 ̅ (𝑓𝑓) = 𝑢𝑢 �(𝑓𝑓) 𝚤𝚤̅ (𝑓𝑓) (1)
Finally, 𝑍𝑍 ̅ is a complex vector which can be displayed in a Nyquist plot in order to visualize the influence of the parameters. This frequency-domain analysis requires that the system under investigation is linear or at least can be linearized with respect to signal amplitudes in its working point and is stationary or at least quasi-stationary for the time of measurement. Fitting experimental EIS with an electrical equivalent circuit is a well-known analysis method [START_REF] Troltzsch | Characterizing aging effects of lithium ion batteries by impedance spectroscopy[END_REF][START_REF] Andre | Characterization of high-power lithiumion batteries by electrochemical impedance spectroscopy. II: Modelling. Selected papers presented at the 12th Ulm ElectroChemical Talks (UECT)[END_REF][START_REF] Kim | An innovative approach for characteristic analysis and state-of-health diagnosis for a Li-ion cell based on the discrete wavelet transform[END_REF][START_REF] Mingant | In-situ battery giagnosis method using electrochemical impedance spectroscopy[END_REF][START_REF] Bernard | Methode de diagnostic in-situ de batteries par spectroscopie d'impedance electrochimique[END_REF][START_REF] Galeotti | Synthetic methods for the evaluation of the State of Health (SOH) of nickel-metal hydride (NiMH) batteries[END_REF][START_REF] Monem | Lithium-ion batteries: Evaluation study of different charging methodologies based on aging process[END_REF][START_REF] Waag | Experimental investigation of the lithium-ion battery impedance characteristic at various conditions and aging states and its influence on the application[END_REF]. This fit is generally done with a Levenberg-Marquardt or a simplex method minimizing least-squares between the model and experimental data. The best-fit parameters obtained by this method can then be used to determine an appropriate SoH indicator by mathematical combination to obtain a linear correlation with SoH.

Another method [START_REF] Mingant | In-situ battery giagnosis method using electrochemical impedance spectroscopy[END_REF] directly uses an electrical equivalent circuit model that depends only on State-of-Charge (SoC) and SoH. Using this method, SoC and SoH can be directly determined by fitting the measured electrochemical impedance spectrum. Several unconventional methods have been developed to measure an EIS spectrum without requiring the often time-consuming application of a single sinusoidal signal at successive frequencies. For example, a so-called "multisine" method comprising several frequencies superposed on a single composite signal can be used to determine an entire spectrum in the time required for the lowest frequency point, giving an almost instantaneous EIS response [START_REF] Breugelmans | Odd random phase multisine electrochemical impedance spectroscopy to quantify a non-stationary behaviour: Theory and validation by calculating an instantaneous impedance value[END_REF]. This technique requires mathematical tools such as Fourier transform and a sophisticated signal generator for the input signal. A faster technique called dynamic EIS (DEIS) measures the impedance response of lithium-ion batteries during charge/discharge at finite DC currents [19][20-22]. In this technique, the impedance measurement during charge/discharge, the lithium-ion batteries are subjected to a current consisting of a constant direct current and a small alternating current (AC) perturbation.

Furthermore, methods employing a single current step have been studied [START_REF] Gabrielli | Measurement Time Versus Accuracy Trade-Off Analyzed for Electrochemical Impedance Measurements by Means of Sine, White Noise and Step Signals[END_REF], and another technique using white noise has been tested by Gabrielli et al. [START_REF] Gabrielli | Comparison of sine wave and white noise analysis for electrochemical impedance measurements[END_REF]. These methods require an input signal which need a battery to be in a rest period, and apply a perturbation on laboratory conditions. Another method, developed by Alavi et al [START_REF] Alavi | Time-domain fitting of battery electrochemical impedance models[END_REF] consisted of estimating simple equivalent circuit model parameters directly from time-domain data. However, application of the method to estimate the EIS model parameters using more complex models is challenging.

In order to avoid the additional burden of an input signal, we have proposed to make use of the free signals of current and voltage from acceleration and regenerative braking that occur during a vehicle road trip [START_REF] Mingant | Procédé non intrusif de determination de l'impedance électrique d'une batterie[END_REF]. These solicitations of the battery can be very complex and indeed can be considered as an almost random profile.

Hence, the current and voltage battery signals could be used to determine a "quasi-EIS" (QEIS) battery spectrum using the method described in previous work [START_REF] Mingant | Procédé non intrusif de determination de l'impedance électrique d'une batterie[END_REF]. It is worth noting that batteries are usually operated at currents so large that the nonlinearity of transport and reaction processes cannot be neglected. Secondly, batteries change their structure while being discharged or charged, and they are typically operated until at least one species of reactants depletes non-stationarity. For this reason this technique is not considered like EIS. It is expected that this method could be used as an onboard diagnostic tool.

This paper first introduces the vehicle application considered in this work, followed by a specific battery ageing protocol that has been developed to simulate ageing for two vehicle architectures. The QEIS data processing including the electrical modeling with an equivalent circuit is described afterwards and applied to the ageing data generated in our laboratories. Finally a correlation between the capacity losses and QEIS model parameters is established to discuss the practical performance of this novel method as an in-service SoH indicator.

Experimental

Case study

The target application is an electric school bus whose mission is to pick up children in the morning on the way to school, return to the bus station in order to charge the traction battery, and then repeat its mission in the opposite direction in the afternoon. Hence the battery charges and discharges twice a day. The first task in defining the bus demonstrator architecture was to size the Li-ion battery, considering two different cases:

-An 85 kWh battery (so called "one pack") which is sized to discharge until a state of charge of 14 % is reached;

-A 170 kWh battery (so called "two packs") which is oversized for the application in order to compensate the capacity loss during ageing.

For these two cases, in-service current profiles were created using an electric vehicle simulator [START_REF] Badin | Evaluation of EVs energy consumption influencing factors, driving conditions, auxiliaries use, driver's aggressiveness[END_REF] in order to simulate the morning and afternoon school bus tours. These current profiles were then used to define the ageing protocol.

Materials

Commercial Li-ion batteries (1.4 Ah cylindrical LiFePO 4 /graphite Valence Technology cells) were selected at the beginning of the project to reversibly store electrical energy in the electric school bus. The nominal cell voltage is 3.3 V and voltage limits are 2.5 V and 3.6 V. Two modules of 66 cells in parallel were available for the study.

Ageing method

Ageing campaigns were launched at module level to collect battery voltage and current data with time, at various SoH. The ageing protocol is intended to provoke both capacity fading and an increase in internal resistance in the batteries. Thus, regular measurements of battery capacity and internal resistance, so called check-up, are mandatory alongside ageing campaigns. The check-up procedure involves significant perturbations on the battery, so in order to limit the impact of the check-up phases on the battery aging, 3 cells were removed from each module to perform a full check-up at cell level and get a statistical determination of SoH from three samples. Then, the ageing protocol continued on the remaining cells of the battery pack with the current profile adapted for the remaining number of cells.

Two electrical test benches were used at IFPEN facilities: a Digatron power test bench 200 A/50 V (+/-10 kW) associated with a large climatic chamber was used to age battery modules submitted to dynamic current cycles under controlled ambient temperature. Module voltage and temperature were measured by the test bench and recorded on a computer. The periodic check-up experiments were carried out on individual cells using a Biologic VMP3 multi-channel potentiostatic-galvanostatic system equipped with an impedance module (BioLogic Science Instruments, France).

Figure 1 presents the current profiles simulating the morning tour for the "one pack" and "two packs" cases.

Each module under test was placed into the climatic chamber set to 45 °C and connected to the battery test bench system. The figure 2 represents the cycling procedure of the two modules. On this figure, The full charge is done at C/3 and C/6 current rate for respectively the "one pack" and "two packs" cases until the voltage reaches 3.6 V.

Then, the battery is charged at current voltage until the complete time of charging of 2h35. In this figure, the cycle profile 1 is the morning cycle, and the cycle profile 2 the afternoon profile The cycling procedure was looped for a total of 7 months for each module. Every month, both modules were submitted to a check-up, and 3 cells were removed from both modules and individually submitted to a check-up.

QEIS determination

The principle of EIS is to measure the complex electrochemical impedance of a cell as a function of frequency.

The conventional method consists of the excitation of the electrochemical cell with a sinusoidal current (or voltage) and the measurement of the voltage (or current) response. By applying a small amplitude input signal a linear relationship can be assumed for the response. This procedure is repeated for each frequency studied. Thus, the complex impedance Z is calculated with the formula (1) :

𝑍𝑍(𝑓𝑓) = |𝑍𝑍|. 𝑒𝑒 𝑖𝑖𝜑𝜑 𝑧𝑧 = � 𝑈𝑈(𝑓𝑓) 𝐼𝐼(𝑓𝑓) � 𝑒𝑒 𝑖𝑖𝜑𝜑 𝑧𝑧 (2)
Where f is the frequency, U is the voltage amplitude, I, the current amplitude, and φ z the phase shift between both U and I.

The power spectral density (PSD) is a mathematical function generally used to represent the different spectral components of a signal. It is equal to the square of the modulus of the Fourier transform of X(t), denoted X(f) reported in half the acquisition time (2):

2

) ( 2 ) ( f X T f x = Ψ ( 3
)
Where T is the integration time.

The PSD represents the distribution of frequencies present in the signal, hence a perfect white noise signal would constitute a constant PSD at all frequencies. There are also cross-power spectral densities which are composed of a combination of two Fourier transforms X(f) and Y(f) (3):

) ( ) ( 2 ) ( * f Y f X T f xy = Ψ (4)
An example of PSD curves measured during road cycling of a battery is presented in figure 4. As shown in this figure, the PSD of this road-cycling protocol consists of two parts: a low-frequency plateau ( f < 0.01 Hz), followed by a decrease at increasing frequencies.

The impedance Z is classically calculated using the formula (4), rearranged for its calculation from PSD(?):

) ( ) ( ) ( ) ( 2 ) ( ) ( 2 ) ( ) ( ) ( * * f f f I f I T f I f V T f I f V f Z I IV Ψ Ψ = = = (5)
However, this formula does not practically provide a good accuracy of the impedance, because of the noise of the measurement. It is therefore necessary to divide the signal into N blocks, then to take the average before calculating the impedance [START_REF] Zhai | The Application of the EIS in Li-ion Batteries Measurement[END_REF].

∑ ∑ ∑ ∑ = = = = Ψ Ψ = = N 1 j I N 1 j UI N 1 j * j j N 1 j * j j ) ( N 1 ) ( N 1 ) ( I ) ( I T 2 N 1 ) ( I ) ( V T 2 N 1 ) f ( Z f f f f f f (6)
This division of blocks causes a reduction of the number of studied frequencies, because of the decrease in the amount of information processed at any time. Thus, for a sample of n values, measured at constant time step, ∆t, the frequency f is a vector comprising the numbers from 1 to n, divided by the duration of the test. Thus, a reduction in sample size also results in a reduction of the number of frequencies by decreasing the number of low frequencies.

In addition, according to the Nyquist-Shannon theorem, higher frequencies must be removed to ensure accurate rendering of the signal by limiting aliasing effects (which reduces even more the number of frequencies of study).Then, the number of studied frequencies is reduced by at least a factor of 2.

There is an important source of uncertainty to consider arising from the signal. This analysis method ideally requires a white noise input signal, but the real-world power cycle in this application is clearly not an ideal white noise, and thus the power spectral density is not constant with frequency (as can be seen in Figure 4). The implication of this is that the impedance spectrum calculated from the ratio U (f) / I (f) contains a significant random uncertainty. This can be seen in figure 5a which shows the result of a direct calculation of Z from the 80 000 values of Ψ(U)/Ψ(I) in figure 4. The basic shape of a semicircular EIS spectrum can be imagined in this figure but there is an unacceptable level of noise in the pattern. Thus, to improve the correlation with a conventional EIS result, a filter system of the power spectral density has been developed. This filter consists to select frequencies of which the current is higher, so dissociated from the noise (figure 5b). Hence, the filter selects the impedance of which the index correspond to a PSD current (PSDI) greater than (PSDI max) / 10. In order to increase even more the relative accuracy of the electrochemical impedance, it was then proposed to cut the signal to reduce the number of points and the electrochemical noise. One solution is to treat the same signal several times, but with a "cutting" becoming smaller . As an example for the voltage signal (x values) on figure 6 : the first cutting presents 2 sections of x/2 values each(and then processing with the previous algorithm), the second cutting presents 4 sections of x/4, the third one, 8 sections of x/8, etc, the end of the cutting represent the one in which the size of the section is under 32 values. It is worth noting that this type of treatment does not allow obtaining impedance spectra on the same frequency range as a standard measurement made from sinusoids (figure 7).

Accuracy study

In order to develop a diagnosis tool, it is necessary to quantify the measurement accuracy to assess the performance of the method. In our study, given the measurement accuracy of both the input current and the output voltage, the measurement error can be considered as a parasitic noise. Thus, according to Gabrielli et al. [START_REF] Gabrielli | Comparison of sine wave and white noise analysis for electrochemical impedance measurements[END_REF] , "the presence of parasitic noise for a low-amplitude signal either increase the measurement time or reduce the measurement accuracy". In this paper, the relative estimation errors ϵ(�Z � �), which will simply be called errors, are defined by :

ϵ²(�Z � �) ≤ 1 2N � 1 η i + 1 η v � (7)
Where η i and η v are respectively the signal-to-noise ratios (current) at the input and at the output, N is the number of blocks. The accuracy measurement is ϵ v =10 mV for the voltage and ϵ i = 0.5 A for the current. The signal-to-noise ratios can be defined by :

η i (f) = |I(f)| ϵ i and η v (f) = |U(f)| ϵ v (8)
Where U(f) is the voltage amplitude and I(f), the current amplitude.

Figure 9 shows an impedance spectrum of an LFP/C 1.4 Ah cell obtained from the QEIS. The calculated error bars according to equation [START_REF] Buller | Modeling the dynamic behavior of supercapacitors using impedance spectroscopy[END_REF] are depicted in red. On this plot, one can observe that the error on impedance measurement is larger at high frequency than at lower frequency, ranging respectively from ± 4 mΩ to 2 mΩ.

Indeed, according to the PSD plot on figure 4, the PSD is higher at low frequency leading to a higher signal-to-noise ratio.

In conclusion of this part, on the whole range of QEIS measurement, a maximal error of 4 mΩ represents a maximal relative error of 10%. It permits us however to fit correctly the shape of the impedance with an equivalent circuit based model.

QEIS data processing (Equivalent circuit based model)

A generic electrical equivalent circuit was used to model Li-ion battery EIS spectra and is shown in figure 9.

This circuit includes one ohmic series resistance R0 which represents electrolyte resistance . The charge transfer resistance is represented by a parallel branch consisting of a capacitance C1 and a resistance R1. Finally, diffusion is represented by a constant phase element (CPE) in which α and Q1 are additional fit parameters. In this way the dynamic electrical characteristics of the battery can be modelled using basic electrical components, like resistors and capacitors, and the fitting parameters can be empirically correlated to the battery SoH.

Results

Ageing campaign

Figure 10 shows the voltage values of modules measured during the ageing cycle in both pack cases under study.

In these plots, lines can be observed every 30 days which correspond to full discharge and charge between check-ups. The first plot relates to the "two packs" case study, and exhibits almost no voltage evolution between the first and the last cycling, implying that there is no significant aging. By contrast, the second plot obtained for the "one pack" case study shows that the minimum discharged voltage of the module decreases to reach the minimal rated voltage of the cell, namely 2.5 V, after only 55 ageing days. At this point, the full battery would not be able to fulfil the theoretical bus mission, but for the sake of this study the ageing test was continued. For this test the cycling procedure was the same, but the battery discharge stopped at a battery voltage of 2.5 V.

figure 11 shows the capacity fade for the modules in the two different sizing scenarios. After an initial increase in capacity, both modules follow an approximately linear capacity fade with time. Consistent with the previous figure, the module representing the "one pack" case exhibits the fastest degradation, with 20% of capacity fade after 210 ageing days, while the "two packs" module undergoes a less severe degradation with 10% of capacity fade after 210 ageing days.

QEIS compared to EIS during ageing

In order to discuss the impact of SoH on electrochemical impedance spectra, EIS and QEIS spectra have been graphically compared in figure 12 and figure 13 using the Nyquist representation on the same frequency range.

The EIS measurements were taken from the individual cells removed for full check-up every 30 days in figure 12, while the QEIS spectra of the modules from three stages of the aging test are shown in figure 13. One can observe that the spectrum moves to the right side as the battery ages, which is indicative of an increase in the ohmic resistance R0 as reported in the literature for classical Li-ion aging characteristics [START_REF] Gyan | Experimental Assessment of Battery Cycle Life Within the SIMSTOCK Research Program[END_REF]. According to this previous work, during ageing a surface film irreversibly grows on the graphite electrode (known as the Solid Electrolyte Interphase). This film is electrically insulating, which accounts for the increase in ohmic resistance.

A comparison between two study cases shows a faster increase of the ohmic resistance with the first configuration than for the two pack case. This reflects a more severe ageing for the one pack case.

Similarly, the QEIS measurements given in figure 13 show that the spectrum moves to the right side as the battery ages, which indicates that a qualitative correlation can be directly observed between the two techniques.

QEIS spectrum modelling

Since there is clearly some inherent uncertainty remaining in the QEIS results, fitting the electrical equivalent circuit model could lead to erroneous results if all the parameters were left free to vary in the fitting algorithm.

For this reason, some parameters were set constant as a function of SoH (table 1), based on the fitting results obtained from all the QEIS measurements. R1 and Q1 parameters have been allowed to vary with SoH, this low frequency impedance is related to a combined effect of the diffusion of lithium ions on the electrode-electrolyte interfaces, which corresponds to the straight sloping line observed in the Nyquist plot towards low frequencies [START_REF] Zhang | Electrochemical impedance study on the low temperature of Li-ion batteries[END_REF].

SoH diagnostic

SoH is considered as the ratio between the real capacity and the initial capacity. The next step of this study is to compare fitted parameters for both ageing campaigns as a function of capacity fade. 5

The "one pack" case has obviously induced a more severe aging effect on the module than the "two packs" case, which is expected due to the substantially higher mean current and larger SoC range encountered in the "one pack" case study, and as already noted in Figure 11. figure 15 shows the increase of the ohmic resistance R0 for the two aged modules as a function of capacity fade, which is thought to be related to SEI growth as already 10 discussed. figure 16 shows the decrease in the Q1 coefficient of the CPE element as a function of capacity fade, Fitted parameters R0 and Q1 both show a monotonic evolution with capacity fade but opposite trends. However, these fitted parameters cannot satisfactorily be directly used to predict capacity fade: for R0, the two curves diverge between 0 % and 5 % of capacity fade; and for Q1, the curves are parallel, but not superposed. In order to have an indicator, the combination, 𝑅𝑅 0 𝑄𝑄 1

, was considered as a function of capacity fade as illustrated in figure 17. According to this figure, the superposition is reasonably effective for both current profiles, revealing that a unique mechanism governs the aging features under consideration, such as capacity fade, in "one pack" and "two packs" cases. When the "one pack and two packs"

Conclusions

In this study, a new in operando method for estimating the SoH of a Li-ion battery was assessed on a commercial module. The method uses the real current and voltage values of the module during an in-use load cycling scenario as input data to calculate an impedance spectrum, so-called QEIS, which could be fitted to an electrical equivalent circuit model. SoH prediction algorithms have been obtained from a mathematical analysis of the model parameters.

This method has been applied to ageing data collected from two modules undergoing different ageing protocols over 7 months. QEIS spectra calculated using our method were fitted to an electrical equivalent circuit model. It was found that the ohmic resistance (R0) and capacitance (C1) parameters of the model evolved monotonically with capacity fade of the battery, for the two different ageing protocols under consideration. By combining the two parameters an indicator of battery SoH was obtained, measurable using only the free output signals of the battery (current and voltage). This correlation could further be applied to other Li-ion battery systems, provided there is a similar ageing mechanism triggered. This method is a promising solution for an embedded system to monitor battery SoH and also for lifetime prediction. To be used for another Li-ion technology it is however necessary to process to a measurement campaign in order to evaluate adapted fitted parameters of other battery chemistries.

Figure 1 :

 1 Figure 1: ageing current profiles for 1 pack (a), 2 packs case study (b) (generated by IFPEN simulation platform [27])

Figure 3 :

 3 Figure 3 : check-up procedures of individual cells

Figure 4 :Figure 5 :

 45 Figure 4: PSD of during road cycling phase of the one pack ageing

Figure 6 :

 6 Figure 6: Representation of the cutting phases of a signal

Figure 7 :

 7 Figure 7: Impedance spectra of an LFP/C 1.4 Ah cell obtained from the QEIS method with "successive cutting" from a signal [U, I] = f (t) road (frequencies of 0.78 to 0.015 Hz), compared with Nyquist spectrum measured on a similar cell at the same temperature (frequencies from 10 000 to 0.005 Hz).

Figure 8 :

 8 Figure 8: Electrical Impedance spectrum of the LFP/C 1.4 Ah cell from an EV profile, blue (*), QEIS, red (-), error bars calculated with equation 7.

Figure 9 :

 9 Figure 9: Electrical equivalent circuit for QEIS spectra fitting

Figure 10 :Figure 11 :

 1011 Figure 10: Module voltage evolution with time during ageing tests for different current profiles corresponding to both case studies under consideration (left, 2 packs case, right 1 pack case)

Figure 12 :

 12 Figure 12: Conventional EIS spectra of the LFP/C 1.4 Ah cell obtained during check-up at 25 °C with 2 packs and 1 pack cases

Figure 14

 14 Figure 14 illustrates an example of a QEIS spectrum fitted with the equivalent circuit model. One can observe a qualitative satisfying adjustment of the QEIS spectrum with the electrical equivalent circuit as previously defined. Fitted parameters obtained on the module at 25 °C before aging have been gathered in table 2. Later in this paper we will use these fitted parameters in order to determine the module's SoH.

Figure 14 :

 14 Figure 14: Example of Nyquist curve fitted with equivalent circuit model

Figure 15 :

 15 Figure 15: R0 vs. capacity fade (modules).

Figure 16 :

 16 Figure 16: Q1 vs. capacity fade (modules).

  , the proposed indicator could predict module capacity fade with a 2 % uncertainty (calculated as the largest gap between the two curves). It is worth noting that the effect of temperature has not yet been taken into account in this study, but this could be considered using an Arrhenius law to extend the predictability of the R0/Q1 correlation with capacity fade as SoH indicator.

Figure 17 :

 17 Figure 17: R0/Q1 vs. SoH.

Table 1 :

 1 Constant values used on the model to fit interest parameters

	Parameter	R 0 [mΩ]	C1 [F]	α
	Value	22.2	11.4	0.5

Table 2 :

 2

	1 pack	22.2	17.8	297

QEIS results obtained on the module @25 °C before ageing

R 0 [mΩ] R 1 [mΩ] Q1 [F 0.5 ]

Acknowledgments

Authors wish to gratefully acknowledge the French FUI for supporting the SCOLELEC project.