domain: spi info phys

text:
Tom Druet
email: tom.druet@cea.fr

Jean-Loup Tastet

Bastien Chapuis
email: bastien.chapuis@cea.fr

Emmanuel Moulin
email: emmanuel.moulin@univ-valenciennes.fr

Autocalibration Method for Guided Wave Tomography with Undersampled Data

Keywords: Structural health monitoring, Guided waves, Ultrasound tomography, Piezoelectric transducers, Baseline-free 15

This paper presents a baseline-free quantitative method for the imaging of corrosion flaws present in thin plates with under-sampled data. This method is based on the Hybrid Algorithm for Robust Breast Ultrasound Tomography (HARBUT) which is itself inherently baseline-free. However, in order to ensure that the incident field component becomes negligible in the reconstruction, a calibration step is necessary. Indeed, it is essential to rescale the data with respect to the acoustic model whether it be simulation data or experimental data. This calibration is usually performed by manually choosing a ray for which the domain of propagation is assumed sound. This can be problematic because this method is not automatic. Moreover, if the chosen ray happens to pass through a flaw, the resulting image will be of poor quality. This paper proposes an autocalibration method for the rescaling step. The field of application is Structural Health Monitoring (SHM) of critical structures with heavy constraints on both sensor intrusiveness and diagnostic reliability. In order to limit intrusiveness, a sub-sampled array of embedded guided waves sensors within the structure is used. Extensions to HARBUT are introduced to compensate for the aliasing caused by the undersampling. The benefits of these extensions are then assessed with numerical simulations and experimental datasets measured by a PZT network.

Introduction

In most industries, corrosion can be disastrous if not correctly monitored. Structural Health Monitoring (SHM) is "the process of acquiring and analyzing data from on-board sensors to evaluate the health of a structure." (see [START_REF]Guidelines for implementation of structural health monitoring on fixed wing aircraft[END_REF] page 4). Guided wave tomography is a promising SHM solution to monitor corrosion. Actually, it is an imaging method studied since many decades, allowing to estimate the remaining wall thickness of corrosion patches in plates, 5 pipes, or more complicated structures. Several algorithms, based on different physical phenomena, have been developed, the older being straight-ray tomography [START_REF] Jansen | Lamb wave tomography[END_REF][START_REF] Mckeon | Parallel projection and crosshole Lamb wave contact scanning tomography[END_REF][START_REF] Malyarenko | Fan beam and double crosshole Lamb wave tomography for mapping flaws in aging aircraft structures[END_REF][START_REF] Belanger | Feasibility of low frequency straight-ray guided wave tomography[END_REF]. This class of algorithms assumes the velocity field to be very weakly non-uniform, thus neglecting refraction and diffraction. Another ray-tomography algorithm, named bent-ray, takes into account the wave field refracted by the potential defect [START_REF] Li | Fast marching method to correct for refraction in ultrasound computed tomography[END_REF][START_REF] Hormati | Robust ultrasound travel-time tomography using the bent ray model[END_REF][START_REF] Li | Refraction corrected transmission ultrasound computed tomography for application in breast imaging[END_REF][START_REF] Willey | Guided wave tomography of pipes with high-order helical modes[END_REF]. This means that very large, non-diffracting defects can be imaged. A different approach is proposed by diffraction tomography (DT) algorithms [START_REF] Wang | Plate-wave diffraction tomography for structural health monitoring[END_REF][START_REF] Simonetti | From beamforming to diffraction tomography[END_REF][START_REF] Rohde | A computer simulation study of imaging flexural inhomogeneities using plate-wave diffraction tomography[END_REF][START_REF] Belanger | Guided wave diffraction tomography within the Born approximation[END_REF][START_REF] Müller | The theory of diffraction tomography[END_REF]. These 10 algorithms are based on the Born approximation and allow to image small defects. More recently, an hybrid algorithm named Hybrid Algorithm for Robust Breast Ultrasound Tomography (HARBUT) has been proposed by Huthwaite and Simonetti [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF] in 2013 and was then followed by many developments [START_REF] Huthwaite | Evaluation of inversion approaches for guided wave thickness mapping[END_REF][START_REF] Huthwaite | Improving accuracy through density correction in guided wave tomography[END_REF][START_REF] Huthwaite | Eliminating incident subtraction in diffraction tomography[END_REF][START_REF] Huthwaite | A new regularization technique for limited-view sound-speed imaging[END_REF][START_REF] Huthwaite | Guided wave tomography with an improved scattering model[END_REF].

The work presented in this paper proposes several extensions to HARBUT in order to satisfy the heavy constraints imposed by the SHM framework, such as sensor intrusiveness and diagnostic reliability. To satisfy the former con-source f i (t) response ϕ i j (t) Defect involving a variation of the parameters p ∈ Υ Ξ : Laws governing the medium Ω Direct problem: find ϕ s.t.

ϕ = Ξ ( f, p, q)
Inverse problem: given ϕ i j , find p s.t.

p = Ξ-1 f i , ϕ i j , q
The medium Ω All the couples ( f i (t), ϕ i j (t)) Positions of the sensors Fig. 1. Inverse problem to solve with guided wave tomography.

preserving the robustness of the algorithm as the number of sensors decreases is presented in this paper. Concerning the latter constraint (the diagnostic reliability), it is well known in SHM that baseline-free methods could decrease tremendously the number of false alarms, typically caused by temperature changes or the ageing of sensors between two measurements. This paper presents an autocalibration method coupled to HARBUT. HARBUT is itself inherently baseline-free but the calibration step is usually performed either manually by choosing a ray for which the domain of propagation is assumed sound or by using a baseline measurement which increases the risk of false alarms. The calibration method presented in this paper is automatic and thus adapted to the SHM context. This paper is organized as follows. Section 2 presents the extensions of the original iterative HARBUT described in [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF], such as the convergence criterion and several added or modified regularization methods (Gaussian filter, threshold based on physical considerations and variable relaxation). Section 3 presents tomographic reconstructions based on the original and modified algorithms with simulated data which take into account the PZT transducers behavior. Three defects involving either refraction or diffraction phenomena are considered. Finally, in Section 4, the method is validated on experimental data measured on a real 2 mm thick aluminum plate. Note that throughout this paper all tomographies are performed with the A 0 mode.

Method

A direct problem consists in finding the response ϕ(x, t) for a given source f (x, t), ∀(x, t) ∈ Ω × R given the laws Ξ governing the medium Ω (Fundamental principle of dynamics, Hooke's law) and the parameters of the problem (q, p) ∈ Υ 0 × Υ for some parameter spaces Υ 0 and Υ (characteristics of the material, geometry). Guided wave tomography consists in finding some of the parameters p ∈ Υ (here, the thickness of the plate) given the laws Ξ-1 governing the medium, the other parameters q ∈ Υ 0 (characteristics of the material) and a partial response ϕ i j (t) for a given source f i (t); it is an inverse problem. Note that Ξ is an invertible regularization of Ξ. Fig. 1 is a representation of the problem considered in this article. Sensors are positioned on a circle. i and j correspond to the sensor indices for an emission i and a measurement j. The partial responses ϕ i j (t) are the inputs to the guided wave tomography. These partial responses are measured for all the couples represented with thin blue lines on Fig. 1. Guided wave tomography consists in inverting these measurements in order to obtain a quantitative image of the plate thickness.

In order to solve the inverse problem, guided wave tomography is usually based on a simplified model of guided wave propagation in order to avoid a prohibitive computation time. An acoustic model is generally used, which is characterized by a scalar wave equation and a dispersion relation ω(k) for the chosen guided wave mode, A 0 in this paper. The original iterative HARBUT [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF] is used here to solve the inverse problem.

The original iterative HARBUT [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF] uses a ray-tomography as a first estimation of the plate thickness, in order to respect the Born approximation in the subsequent DT step. In this section, the differences between the original

v G (m•s -1 ) S 0 A 0 A 1 SH 0 SH 1
Fig. 2. Dispersion curves of the group velocity v G for S 0 , A 0 , A 1 , SH 0 and SH 1 , and for a 2 mm thick aluminum plate (V L = 6360 m s -1 , V T = 3140 m s -1 and ρ = 2.7 g cm -3 ).

HARBUT and our method are presented. Some additional elements are also introduced, such as a regularization method called "variable relaxation", which helps to limit aliasing effects with strongly undersampled data.

Algorithmic differences with the original iterative HARBUT

This section includes two parts; first are presented some differences for the Time-of-Flight (ToF) tomography step and then for the DT one.

Time-of-Flight Tomography step

Whether it be straight-ray or bent-ray, in order to solve the propagation problem, one uses the SART algorithm [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF] (or FMM-SART for bent-ray [START_REF] Li | Fast marching method to correct for refraction in ultrasound computed tomography[END_REF]). The structure to be imaged is represented at the position x by the scalar slowness field s(x) = 1 v G (x) (isotropic material), where v G is the group velocity. Fig. 2 presents the dispersion curves of the group velocity for the configuration presented in this paper (the working product frequency × thickness considered in this paper is very low: 40.8 kHz mm). To represent the slowness s(x), we project it on a finite basis. In the original bent ray algorithm [START_REF] Li | Fast marching method to correct for refraction in ultrasound computed tomography[END_REF], a basis of cubic or spherical pixels / voxels is used for the reconstruction, and a non-interpolating B-spline filter [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF] is applied a posteriori in order to accurately estimate the gradient when computing the rays. Here, we have instead expressed the SART algorithm in terms of B-splines from the very beginning. The projection s(x) of the scalar slowness thus reads:

s(x) = i, j S i j β i j (x), (1)
where β is a 2D B-spline basis and S i j is the coefficient associated to the B-spline β i j .

55

With this decomposition, it is now possible to express the ToF τ k→l along the trajectory (straight or bended) γ k→l between the emitter k and the receiver l:

τ k→l = γ k→l s(γ) dγ γ k→l s(γ) dγ = γ k→l i, j S i j β i j (γ) dγ = i, j S i j A i jkl , (2)
where A i jkl = γ k→l β i j (γ) dγ is the contribution of the basis element β i j to the ToF τ k→l . The use of B-splines allows us to obtain smooth images, which is important for the following DT iterations. In practice, fourth order splines are used. This is a good compromise to get a smooth-enough image while keeping a sufficiently small support for the defect reconstruction.

To perform baseline-free tomography, only ToFs corresponding to the defect data are available: τ defect k→l . Therefore, an absolute estimator needs to be used, and is described here. The n th iteration of the ray-tomography step reads:

S (n) i j = S (n-1) i j + δ k S (n) i j
, where δ k S (n) i j is the correction applied to the coefficient S i j and is calculated according to the following formulation:

δ k S (n) i j = 1 N i j l        W i jkl τ defect k→l -τ (n-1) k→l λ k→l        , (3)
where λ k→l is the trajectory length, W i jkl is equal to A i jkl in the simplest case but usually includes a Hamming window to avoid side effects, and the normalization factor N i j reads: N i j = l A i jkl .

Unlike the absolute estimator, the relative estimator uses the ToF corresponding to baseline data: τ baseline k→l . The relative estimator allows to eliminate some systematic errors due to experimental uncertainties and works on relative ToF. The idea is to subtract independently the experimental ToF (τ defect k→l and τ baseline k→l ) and the calculated ToF (τ k→l and τ healthy k→l ). The relative estimator is then calculated according to the following formulation:

δ k S (n) i j = 1 N i j l        W i jkl (τ defect k→l -τ baseline k→l ) -(τ (n-1) k→l -τ healthy k→l ) λ k→l        . (4)
Those two estimators are compared on experimental data in Section 4.2.1.

Finally, concerning the ray-tomography, a regularization method -presented in Section 2.2 -is applied to the map at every iteration in order to obtain a "cleaner" ray-tomography image for the subsequent DT iterations.

Diffraction Tomography step

Concerning the DT step, the algorithm used in this paper is the same as the original HARBUT. However, several extensions to HARBUT are proposed in order to satisfy the heavy constraints imposed by the SHM framework, such as sensor intrusiveness and diagnostic reliability. Specifically, a Gaussian blur adapted to the undersampled data is applied between iterations and an autocalibration method is used to automate the imaging procedure.

Prerequisite. The same object function O(x) as in [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF] is used to represent the image:

O(x) = k(x) 2 -k 2 0 = k 2 0        v 0 v(x) 2 -1        , (5)
where k 0 is the background wavenumber, v 0 the background phase velocity and v(x) the phase velocity at the position x, which is linked to the thickness thanks to the dispersion curves when the algorithm has converged. Iterative HARBUT consists in calculating iteratively the correction O δ (x), which is added to the previous estimation O b (x) so that O(x) = O b (x) + O δ (x). The correction O δ (x) is calculated by means of a beamforming image for which a filter in the spatial frequency domain (see [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF] for details) is applied. In the HARBUT theory, the input data is the diffracted field but it is shown in [START_REF] Huthwaite | Eliminating incident subtraction in diffraction tomography[END_REF] how it is possible to work with the total field instead of the diffracted field. This is the key step for HARBUT to be free of any baseline.

Autocalibration method. In order to perform HARBUT, it is essential to rescale the data with respect to the theoretical model in order for the incident field component to become negligible in the reconstruction, whether it be simulation data or experimental data. This calibration could be done thanks to a baseline but as explained before, this is unwanted in the case of a SHM application because of the false alarms which could occur. This is why the calibration is usually made with the current measurements, by manually choosing a ray for which the domain of propagation is assumed sound. This is not ideal either in SHM because it is not automatic. Moreover, if the chosen ray passes through a flaw, the resulting image will be of poor quality. Here, an autocalibration method is proposed. This method makes the global imaging process automatic and baseline-free while preserving its robustness.

The theoretical incident field is given by the non-perturbed Green's function: where H (1) 0 is the Hankel function of the first kind and of zeroth order. The idea is to calculate the calibration factors c i j :

G 0 (x; x ) = - i 4 H (1) 0 (k 0 x -x ), (6)
c i j = G 0 (x i ; x j ) ϕ healthy i (x j ) , (7)
where the index i denotes the emitter, the index j the receiver and ϕ healthy the measured incident field. Then these calibration factors are multiplied to the measured data, which provides the input calibrated fields for HARBUT. In the case of tomography relying on baseline signals, the construction of ϕ healthy is straightforward, but for a baseline-free method it is challenging to do it in a robust way. Naively, we could select a couple (i, j) for which the ray crosses no defect and then proceed to the calibration with only this couple. This leads to two problems: the first one is that when experimental uncertainties are too strong this approach could introduce noise in the resulting image; the second one is that it is difficult in practice -even nearly impossible -to be sure that there is no defect on the path of the ray, and this would strongly disrupt the rescaling process. Actually, even if a first ray-tomography using the absolute estimator is performed to approximately locate potential defects, this kind of algorithm is not able to detect diffracting defects. Therefore, it has been decided to develop a method to construct the c i j in a robust way without any baseline. The c i j are constructed in three steps:

1. A ray-tomography with absolute estimator is performed with enough regularization (presented in Section 2.2) to obtain a background without noise. Fig. 3a presents a result of this kind of tomography with all the rays corresponding to all possible couples (i, j). Fig. 4 shows all the potential c i j of Eq. ( 7) for baseline data (black triangles) and for defect data (yellow squares). Obviously, the baseline calibration factors are in close proximity to each other, and concerning the potential calibration factors corresponding to the defect state, it is possible to see some of them which diverge from the others. They correspond to the rays which cross the defect. Moreover, the problem is axisymmetric, which implies that if n sensors is the number of sensors in the distribution, for each visible point there are n sensors points superimposed. The non axisymmetric case is treated in Fig. 19 for experimental data. 2. For each ray, we verify that the ray is not crossing any visible defect on the ray-tomography along its path (in practice, a threshold of 1 % relative to the healthy thickness has been used). If a ray crosses a defect, it is not used in the following step. This means that m rays are kept with m ≤ n rays where n rays = n sensors (n sensors -1)

2 the total number of rays. Fig. 3b shows the m remaining rays. Fig. 4 shows the corresponding c i j (red circle). This time, the factors kept are in close proximity to each other. It is possible to notice that only 2 n sensors factors were discarded (two red circles are not superimposed on yellow squares) despite being in close proximity, which suggests that the threshold may be too strict, but this is preferable to the contrary.

3. Finally, a confidence ellipse is computed based on the factors corresponding to the rays within the healthy area (red circles). As noticed before, the ray-tomography algorithms are blind to diffraction. This means that it is possible for a small defect to be within the healthy area, intersecting some of the kept trajectories. Fortunately, the associated calibration factors should diverge from the "proximity area" and thus be outside the confidence ellipse.

In this paper, a 95 % confidence ellipse has been chosen but we have tested smaller values, down to 5 %, and the global error of the resulting image does not change significantly. Only the p factors c i j are kept (p < 2 m because one ray is composed of two sensors) which are inside the computed confidence ellipse. These p factors are named c ellipse and satisfy c ellipse ∈ C p ⊂ C 2 n rays . They are represented by small green circles in Fig. 4. The 2 n rays -p remaining factors to be calculated are named c defect and are estimated as follows:

c defect = p k=1 c ellipse k p . (8)
These three stages allow to compute all the needed calibration factors to perform baseline-free guided wave imaging with iterative HARBUT.

Gaussian filter. In the usual iterative HARBUT, a Gaussian blur is already applied between DT iterations. It is a kind of regularization. In this paper, a Gaussian blur is applied not only between DT iterations, but also for ray tomography.

Here the importance of this filter when the number of sensors decreases is discussed.

Whether it be for ray-tomography or the DT iterations, it has been decided to adjust the Gaussian filter standard deviation to the algorithm resolution. For ray-tomography algorithms it is well known that the resolution corresponds to the first Fresnel zone √ λL where λ is the wavelength and L is the propagation distance between the sensors, so for our configuration L = 2R where R is the radius of the sensor distribution. The effect of the defect size on the forward propagation when considering ray-theory is presented in Fig. 5 for three defects studied in Section 3. The √ λL length is also plotted for scale. Regarding the HARBUT iterations (DT steps), [START_REF] Huthwaite | Evaluation of inversion approaches for guided wave thickness mapping[END_REF] explains that the resolution of HARBUT is around 1.5-2λ and [START_REF] Huthwaite | Eliminating incident subtraction in diffraction tomography[END_REF] demonstrates that the spatial sampling must be carried out such that ∆ lim = λ 2 where ∆ lim is the distance between two sensors. If ∆ lim is smaller (more sensors), the quality of the results remains the same but the computation time is increased, while if ∆ lim is longer (less sensors), artifacts appear caused by aliasing.

In order to eliminate non-physical oscillations (above the theoretical resolution), the Gaussian filter standard deviation is adapted as follows. A Gaussian filter is the convolution between the image data and a Gaussian of standard deviation σ. If we consider the 1-D case, the Gaussian reads

f (x) ∝ exp -(x-µ) 2 2σ 2
and its Fourier transform

f (k) ∝ σ exp (-ikµ) exp -σ 2 k 2 2 . The normalized Gaussian is f (k) = exp -k 2 2 σ2 = exp -(2π) 2 2 σ2 λ 2 with k = 2π λ and σ = 1
σ . It has been decided to fix the relative remaining spectral density to 1 % for the spatial wavelength equal to half the theoretical resolution. In practice, this means that 29.1 % of the spectral density remain for the spatial wavelength equal to the theoretical resolution. For our configuration, if λ w is the working wavelength, the theoretical ray-tomography resolution is is 2λ w . After some basic calculations, we find that σ ray-tomo = λ w and σ HARBUT = λ w 2 . When the number of sensors decreases, σ HARBUT has to be adapted according to the new value of ∆ lim . For instance, if half the optimal theoretical number of sensors is used, the standard deviation should be σ HARBUT = λ w . It is shown in Section 3.2.2 that it is also possible to keep the standard deviation of λ w 2 if other kinds of regularization are used. It is worth noting that it would have also been possible to filter the disruptive spatial frequency components in the dual space (K-space) as explained in [START_REF] Huthwaite | Eliminating incident subtraction in diffraction tomography[END_REF], but the very smooth Gaussian filter behavior is also interesting because it acts as a relaxation in the iterative process.

√ λ w L = √ 31
We have also tried to decrease the standard deviation at every iteration. The achieved resolution is better but the image is noisier. The regularization method called "variable relaxation" presented in Section 2.2 is then needed. For the sake of simplicity, we have decided to keep a constant standard deviation for the Gaussian blur in this paper.

Convergence criterion.

Contrary to what was used in [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF], we have decided to formulate the convergence criterion using the total map (not only the damage area) to avoid problems which could occur when the noise present in the map becomes bigger than the threshold chosen to identify the defect area. The Q parameter, allowing to define the convergence, then reads:

Q (n) = e (n) (x) -e (n-1) (x) dx e plt dx , (9)
where e plt is the healthy plate thickness and e (n) (x) is the thickness at position x for iteration n. Convergence is considered to be achieved when Q ≤ 2 × 10 -4 , which is different from the criterion chosen in [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF], and this is justified because [START_REF] Huthwaite | High-resolution guided wave tomography[END_REF] considers the defect area only. Moreover, in order to possibly work with very noisy data, another criterion is added. If Q is increasing more than M consecutive iterations, the algorithm stops and convergence is considered to be reached. From experience, we have observed that M can be chosen to be 2 or 3.

Regularization

In order to perform baseline-free guided wave tomography with an undersampled sensor array, it is necessary to add prior information on the set of possible parameters p, otherwise the inverse problem Ξ -1 would be ill-defined, leading to aliasing and bad convergence. This is accomplished through the use of several regularization methods [START_REF] Tikhonov | Solution of incorrectly formulated problems and the regularization method[END_REF]. Regularization is classically implemented by adding a penalty term to the Lagrangian [START_REF] Wysoczański | Performance analysis of regularization algorithms used for image reconstruction in computed tomography[END_REF]. However, this requires to explicitly formulate the tomography as a minimization problem, which adds complexity to bent ray [START_REF] Hormati | Robust ultrasound travel-time tomography using the bent ray model[END_REF] and would be impractical -if not impossible -for HARBUT. Instead, here we implement the regularization implicitly on top of the usual HARBUT algorithm, as additional steps inserted between iterations. This allows us to use the same regularization methods both for the ray and diffraction tomography steps of HARBUT.

Each regularization step operates exclusively on the partial image (weights e i j or map e (n) (x)), independently of the input ϕ i j , in an attempt to enforce some of its properties. If the prior information introduced by the implicit regularization method accurately represents the properties of the SHM setup (such as material properties, possible defects or noise profile), it can enable us to obtain a more accurate reconstruction, or alternatively to maintain the same accuracy while using less sensors, in the spirit of sparse tomography methods [START_REF] Wang | Méthodes de reconstruction d'images à partir d'un faible nombre de projections en tomographie par rayons x[END_REF]. On the other hand, if the prior does not faithfully represent the properties of the setup, then the regularization method can introduce a bias. It is therefore necessary to verify that the prior information is accurate for the specific SHM setup. This is done by validating the model using physical or simulated test defects, which must be representative of the expected defects when the system is later used in production. In principle, any prior can therefore be used provided that it is successfully validated. However, for the present work, we will use our knowledge of the setup to guide our choice of regularization method.

The Gaussian blur has already been presented in Section 2.1.2. Two others methods are used in this paper. The first one is based on a physical constraint: we assume that the thickness can never increase because of corrosion. This is implemented as a constant threshold equal to e plt . The second one, which we call "variable relaxation", is inspired by Tikhonov regularization [START_REF] Tikhonov | Solution of incorrectly formulated problems and the regularization method[END_REF] (which limits the effective number of weights) and weight elimination [START_REF] Weigend | Generalization by weight-elimination with application to forecasting[END_REF] (which further reduces small weights, hence limiting noise at the cost of a lower sensitivity to small defects), and is applied to the thickness map for ray-tomography and to the object function for HARBUT after each iteration. It limits the amplification of spatial fluctuations at each iteration due to noise and the lack of sensors, while still allowing us to accurately reconstruct localized defects above the noise threshold. If we call e i j an element of the image, the "variable relaxation" is defined as follows:

e i j =                ẽ + e i j -ẽ 1 + 1 z 2 i j α 2 if 0 < e i j -ẽ < γβẽ, e i j otherwise, (10)
with e i j the new value of the element, ẽ the reference of the regularization (for example ẽ = e plt for the ray-tomography step), β ∈ [0, 1] a threshold which has to be fixed depending on the amplitude fluctuation of the noise present in the image, γ ∈ [1, +∞[ a constant fixing the relaxation limit, α ∈ R + a constant defining the strength of the regularization and z i j reads:

z i j =                      e i j -ẽ βẽ if e i j -ẽ < βẽ, e i j -ẽ 1 2 βẽ 1 -cos π |eij-ẽ|-γβẽ (1-γ)βẽ if βẽ ≤ e i j -ẽ < γβẽ. (11)
In order to clarify the effects of β, γ and α, several values are tested in Fig. 6 for a plate of thickness e plt = 2 mm, in the case of time-of-flight tomography.

Domain of effect of the regularization:

The effect of the regularization is limited. Indeed, e i j e i j , ∀e i j ∈ (1γβ)ẽ, ẽ and e i j = e i j elsewhere. This can be verified with the red curves in Fig. 6: for the solid red line e i j e i j , ∀e i j ∈ ]0, 2[ and e i j = e i j elsewhere whereas for the dashed red line e i j e i j , ∀e i j ∈ ]1, 2[ and e i j = e i j elsewhere. Thereby, the domain of effect of the regularization is twice larger for the solid red line than for the dashed one.

Effect of the threshold β: The effect of the regularization is low for e i j ∈ (1γβ)ẽ, (1β)ẽ and high for e i j ∈

(1β)ẽ, ẽ . To illustrate that, the solid blue line (scenario 2) and the solid red line (scenario 4) are compared.

Scenarios 2 and 4 have the same regularization domain, i.e. γ 2 β 2 = γ 4 β 4 , the same strength, i.e. α 2 = α 4 , but two different thresholds, i.e. β 2 β 4 . We can see that the effect of the regularization is high in a larger domain for the scenario 4 than for the scenario 2. β has to be adapted to the amplitude fluctuations of the noise present in the image. In real applications, it is possible to estimate the amplitude of the unwanted fluctuations present in the image. This noise can be due to several factors such as an error in the positioning of the sensors, a lack of sensors which causes aliasing, noise in the data or an approximation in the model. In this paper, the unwanted fluctuations have been estimated to be around 10 % of the reference thickness. So it has been decided to fix β = 0.1. Effect of α: Finally, the parameter α is studied thanks to the blue curves in Fig. 6. This parameter is the strength of the regularization: the higher is α, the stronger is the effect of the regularization.

These three hyperparameters have been studied in several configurations and we have noticed that the optimal behavior is obtained when γβ = 1. It is interesting to note that this result allows us to fix only two parameters instead of three.

We have fixed β = 0.1 as explained before. Concerning the strength, we decided to fix α = 4 for the ray tomography because we need a very clean map for the autocalibration step, and α = 2 for the DT step of HARBUT. We expect that hyperparameters will need to be adjusted to each specific SHM setup for optimal performance. However, as long as they only depend on the general properties of the setup, which are stable in time, our regularization method remains compatible with baseline-free tomography.

Global relative error

In order to quantitatively compare tomographic images, a global relative error is used. Let fi be the i th pixel of the reference map and f i the i th pixel of the tomographic image. The global relative error E global reads:

E global = N pix i=1 fi -f i 2 N pix i=1 f 2 i , (12)
where N pix is the number of pixel in the image. The methods used to construct the reference maps for the numerical results and the experimental results are explained in Sections 3.1 and 4.1, respectively.

Numerical results

Numerical model

In order to test the algorithm, numerical 3D elastodynamic simulations based on the Finite Element Method (FEM) were performed. The FEM code is internally developed in the NDE department of CEA-LIST [START_REF] Imperiale | UT Simulation of Embedded Parametric Defects Using a Hybrid Model Based Upon Spectral Finite Element and Domain Decomposition Methods[END_REF]. The configuration is presented in Fig. 7. An aluminum plate of density ρ = 2.7 g cm -3 , longitudinal velocity V L = 6360 m s -1 , transversal velocity V T = 3140 m s -1 and dimensions 1400 mm × 1250 mm × 2 mm is considered. The sensors are positioned on a circle of radius r distrib = 300 mm, as shown in Fig. 7. To simulate the PZT transducer behavior, we use a second order numeric scheme both in time and in space, as well as a mesh adapted to the emitter geometry (see Fig. 7). In-plane forces normal to the circle of same diameter as the PZT disc are applied to the mesh nodes (see [START_REF] Fribourg-Blanc | Actionneurs en couches minces pour le contrôle santé de structures aéronautiques -étude de films de PZT et de PMNT[END_REF] for the modeling of PZT transducers and [START_REF] Chapuis | Contrôle santé intégré par méthode ultrasonore des réparations composites collées sur des structures métalliques[END_REF][START_REF] Etaix | Thickness or phase velocity measurements using the Green's function comparison method[END_REF] for experimental validations). The size of the elements is set to 2 mm (because the working frequency is 20.4 kHz and the associated wavelength 30.7 mm) along the two in-plane dimensions (excepted for the emitter zone). An axisymmetrical model is used in order to minimize the computation time.

The defect is parametrized using a segment Υ bot of length r bot corresponding to the defect bottom radius, the defect depth e def and a transition zone Υ trans of length l trans defined in terms of Bernstein polynomials b (m) i :

       Υ trans (x trans ) = e plt + e def -1 + 3 i=0 β i b (3) i x trans -r bot l trans ∀x trans ∈ [r bot , r bot + l trans ] , b (m) i (x) = m! i!(m-i)! x i (1 -x) m-i , (13)
where β 0 = 0, β 1 = 0, β 2 = 1 and β 3 = 1. With a revolution of Υ bot and Υ trans with respect to the normal axis to the plate surface and crossing the center of the sensors distribution, the defect geometry is constructed as presented in Fig. 8. The three defects studied in this paper are presented in Fig. 5.

The solution is extracted within a ring and then processed to reconstruct either the PZT response or the normal response. To obtain the normal response, we only need to interpolate the ring solution for each time-step and then extract the out-of-plane displacement for the points of interest. To compute the PZT response f rec PZT , it is necessary to integrate the in-plane displacements u 1 and u 2 of the active part of the PZT disc along the closed curve γ PZT :

f rec PZT (t) = k PZT γ PZT (u(t) • n) dγ(x),
where u ≡ (u 1 , u 2 ) and x ≡ (x 1 , x 2 ). (

) 14
where n is the outer-pointing normal of γ PZT and k PZT is a constant depending on the piezoelectric material parameters (see [START_REF] Fribourg-Blanc | Actionneurs en couches minces pour le contrôle santé de structures aéronautiques -étude de films de PZT et de PMNT[END_REF] for more details). k PZT is ignored here because we only need to consider relative amplitude signal variations for this work. The frequency behavior of a PZT transducer can be described as follows (see [START_REF] Chapuis | Contrôle santé intégré par méthode ultrasonore des réparations composites collées sur des structures métalliques[END_REF][START_REF] Giurgiutiu | Tuned Lamb wave excitation and detection with piezoelectric wafer active sensors for structural health monitoring[END_REF] for more details):

f rec PZT (ω) ∝ J 1 π ∅ PZT λ(ω) , (15)
where ∅ PZT is the diameter of the active part of the PZT disc and J 1 is the Bessel function of the first kind and of first order. Fig. 9a shows such a behavior for our setup.

Finally, before studying the algorithm behavior, experimental signals and simulated signals are compared in Fig. 9b. The excitation is a 5 cycles toneburst of central frequency 20.4 kHz. It may be noted that there is an electromagnetic coupling between the emission and the experimental measurement. The simulated normal response does not contain the S 0 mode. This is due to S 0 having practically no out-of-plane displacement in comparison with the in-plane displacement for our setup (the out-of-plane displacement is around 1 % of the in-plane displacement here). However, S 0 is visible in the simulated PZT response as well as in the experimental signal, both of which use the same PZT geometry. This means that, because of reflexions on the plate sides, S 0 will interfere with A 0 (the packet of interest in this paper).

Reconstructions

Ray-tomography

In this section, we explain why straight-ray is preferred to bent-ray as a first iteration for HARBUT -which seems paradoxical given that bent-ray uses a more accurate model of wave propagation. It should be mentioned that all the ray-tomographies of Section 3 use the absolute estimator, because our final aim is to perform baseline-free guided wave tomography. Moreover, in order to avoid S 0 perturbations, as a first step we will only consider the normal response. Furthermore, no regularization (not even the Gaussian blur) is used in this Section.

First we consider the defect of Fig. 5a, which presents no diffraction but features strong refraction. In this particular case, Fig. 10 shows, as expected, that bent-ray gives better results than straight-ray. However, when the defect of Fig. 5c is considered, Fig. 11 shows that straight-ray has a better behavior when strong diffraction is present. This is because straight-ray is better conditioned than bent-ray. Consequently, it is sometimes better to use straight-ray instead of bent-ray, depending on the kind of defect searched for. In SHM, this problem can be overcome. Indeed, it is possible to record the reconstruction each time an inspection is performed. As the defect is growing, inspections are performed frequently enough such that the n th one can use the (n -1) th reconstructions as background, so that the Born approximation is never violated.

HARBUT

Gaussian blur effect. The defect of Fig. 5b, presenting both refraction and diffraction, is considered. Is should be mentioned that no Gaussian blur is used for the bent-ray tomographies.

Concerning HARBUT, Fig. 12a is a reconstruction without any regularization nor Gaussian blur. We used much more sensors than the optimum number of 120 which is necessary to adequately sample the Ewald circle [START_REF] Huthwaite | Eliminating incident subtraction in diffraction tomography[END_REF]. The algorithm has been stopped at the 7 th iteration as it would otherwise diverge. It is possible to see the 2k 0 spatial frequency oscillations. If a Gaussian blur of standard deviation 1 2 λ (explained in Section 2.1.2) is used, the reconstruction is practically perfect, as seen in Fig. 12b. Then, if the number of used sensors is decreased to the optimal number 120, Fig. 12c shows that the reconstruction has the same quality but the computation time is naturally lower (roughly four times faster). If half the optimal number of sensors are exploited (60) and the same Gaussian blur is applied, Fig. 12d shows again some parasitic spatial oscillations but this time of spatial wavelength longer than the ones present in Fig. 12a. If the standard deviation is increased to 1λ, as explained in Section 2.1.2, these oscillations disappear but the defect reconstruction decreases in quality: see Fig. 12e. Finally, it is possible to see in Fig. 12f that by adding regularization (here, the physical threshold) and by keeping the standard deviation to 1 2 λ for 60 sensors, we can still obtain a reconstruction of good quality.

Baseline-free HARBUT with autocalibration method. For the following reconstructions in this paper, only 30 sensors are used. In order to obtain the best reconstructions, we use a standard deviation of 1λ (and not 2λ) with regularization, i.e. the threshold and/or the variable relaxation. Moreover, the Gaussian filter is always used for the ray-tomography reconstructions.

Baseline-free HARBUT combined with the autocalibration method and HARBUT combined with a classic calibration and with the relative estimator for the time-of-flight tomography step are compared in Fig. 13 for defect 3. Only 30 sensors are used, which corresponds to a 75 % reduction compared to optimal sampling. It should be underlined that this time the PZT response is used. Moreover, only the threshold is used (no variable relaxation). It is possible to see some small artifacts on both reconstructions due to the lack of sensors. If variable relaxation had been employed, those artifacts would have disappeared (this will be shown in the experimental Section 4.2.2). The reconstructions of Fig. 13 correspond to the calibration factors presented in Fig. 4.

Table 1 allows us to quantitatively estimate the quality of the reconstructions. The global relative error (see Eq. ( 12)) is very similar between the autocalibration and classic calibration methods. This study demonstrates the good performance of the autocalibration method on simulated data. The next section is dedicated to experimental studies.

Experimental results

Experimental setup

Two pictures of the experimental setup are showed in Fig. 14. We consider an aluminum plate of density ρ = 2.7 g cm -3 , longitudinal velocity V L = 6360 m s -1 , transversal velocity V T = 3140 m s -1 and dimensions 1400 mm × 1250 mm × 2 mm. The sensors distribution is a circle of radius r distrib = 300 mm. The PZT transducers, of diameter ∅ PZT = 18 mm, are glued to the aluminum plate.

In order to have a reference to calculate the global relative error, an ultrasonic scan is performed in water immersion and the thickness map of Fig. 15 is obtained from the ToFs. The first Fresnel zone and the wavelength are superimposed on it. We can see that the defect to be imaged looks like the simulated one presented in Fig. 5c.

Reconstructions 4.2.1. Ray-tomography

Here we compare the absolute estimator (Eq. ( 3)) to the relative estimator (Eq. ( 4)) in the case of experimental data.

The straight-ray reconstruction using the relative estimator is presented in Fig. 17a. The healthy part of the map is relatively "clean". When the absolute estimator is used (see Fig. 17b) the healthy part is noisier which was expected because the experimental uncertainties cannot be compensated with this estimator. However, in order to limit the noise, it is possible to apply the variable relaxation. As explained in Section 2.2, it has been decided to fix β = 0.1, the product γβ = 1 and α = 4 for the ray tomography step in order to get a very clean background image for the autocalibration step. Fig. 16 shows the curve of Eq. ( 10) for these parameters. The resulting reconstruction is presented in Fig. 17c. This time, the healthy part contains no noise, which is essential for the following HARBUT steps.

HARBUT

In order to perform baseline-free HARBUT with autocalibration, the approach proposed in Section 2.1.2 to calibrate the field is applied to experimental data. First, the baseline-free straight-ray reconstruction is used to find all the healthy rays. This is presented in Fig. 18a for the full ray distribution and in Fig. 18b for the healthy rays. All the calibration factors are presented in Fig. 19a. This time, contrary to the axisymmetrical simulation cases, no calibration factors are superimposed on each others. Moreover, the distinction between the factors corresponding to rays crossing the defect and those travelling through the healthy part is not as clear as with simulated data. It is possible to see factors considered healthy which diverge from the proximity area. This means that the confidence ellipse visible in Fig. 19b plays a full role to remove the factors viewed as healthy but which would cause a poor calibration.

Finally, we present in Fig. 21 the results for HARBUT with the use of a baseline (all the calibration factors are calculated thanks to a baseline) and in Fig. 22 those for baseline-free HARBUT with the autocalibration method.

The cross-sections plotted on the reference (see Fig. 15) are presented in Fig. 21b, 21d, 22b and 22d. Fig. 21a and 22a present the reconstructions without the variable relaxation for the DT steps. On both images, the healthy part presents artifacts because of aliasing due to the lack of sensors. By comparing cross-sections (Fig. 21b and22b), ones can observe that the artifacts are stronger on baseline-free HARBUT with autocalibration. Fig. 21c and Fig. 22c present the reconstructions with the use of variable relaxation presented in Fig. 20a and20b. Again, as explained in Section 2.2, it has been decided to fix β = 0.1, the product γβ = 1 and α = 2 for the DT step of HARBUT. α does not have to be as strong as in the case of ray tomography. The artifacts disappear almost entirely. In addition, Fig. 21d and22d show that while small variations present in the defects (see reference) are not reconstructed because of the too large value of λ and the lack of sensors, the global geometry is, however, well identified.

Finally, the global relative errors are presented in Table 2. This allows us to quantitatively compare the reconstructions. We can see that baseline-free HARBUT with autocalibration presents relative errors of the same order of magnitude as HARBUT with the use of a baseline for the calibration.

Discussion

The two main contributions of this paper are the autocalibration method and a way to deal with undersampled data.

Concerning the autocalibration method, one could think of a situation which could lead to poor calibration. This situation is discussed below. The autocalibration method is based on an image obtained with ray-tomography. If a transducer or a group of transducers is mispositioned, this leads to a local change in the reconstructed velocity very close to these transducers, and thus the concerned rays are not used for the autocalibration. What about the consequences of this scenario? First, in an industrial use of this method, the mispositioning is expected to be sufficiently small that it would not cause any false values within the initial image. The experimental results showed in this paper are of good quality although no particular attention has been paid to the positioning of the sensors (i.e. positioning accuracy ≈ 1 mm). But even if one or a group of transducers is mispositioned on purpose (or have errors), the rays concerned by the velocity change will not be used for the calibration. This is not a problem as long as there are enough rays to perform the calibration. As explained in the paragraph concerning the autocalibration, we even tried to keep only 5 % of the healthy rays (the ones which do not intersect the flaws in the ray tomography image) but the resulting image does not change significantly. Finally, the sensor positioning accuracy required to obtain a high-quality reconstruction would be enough to ensure a good autocalibration. Consequently, the autocalibration method seems to be robust enough to be used in a SHM context.

Concerning the undersampling issue, the idea was to use less sensors than the optimal number in order to be less intrusive in the structure. It has been shown in this paper that we could still obtain an image of good quality with four times less sensors than the optimal number. However, for some applications (such as aeronautics), the density of sensors required remains too high to be of practical relevance. Further studies are being conducted in parallel by the authors (and other teams) to further reduce the required number of sensors. Another promising approach is to reduce intrusiveness thanks to Fiber Bragg Gratings (FBGs). Using this technology, all the sensors can be replaced by one fiber with Bragg gratings written in the fiber. This could be a far less intrusive way to perform SHM. However, FBGs can only measure received guided waves and cannot generate them. In order to overcome this problematic the idea is to use so-called "passive methods" to retrieve the signals that are needed to perform tomography [START_REF] Druet | Passive guided waves measurements using fiber Bragg gratings sensors[END_REF]. A work is conducted in parallel on passive tomography [START_REF] Druet | Passive guided wave tomography for structural health monitoring applications[END_REF], and more generally on passive imaging with FBGs; it will be discussed in a future paper.

Conclusion

This paper has presented extensions to the original iterative HARBUT to deal with SHM constraints. Indeed, these extensions allow diagnostic reliability despite low sensor intrusiveness. This work focused on the use of PZT transducers which can be very intrusive when used in large numbers. It has been necessary to find a way to keep the algorithm robust when the number of sensors decreases below the optimal number which is typically required for proper sampling. In order to study a large number of configurations, a simulation campaign, which can take into account either PZT emission and reception or out-of-plane displacement, has been realized. Then, to confirm the results, experimental tomographies have been studied on real PZT data.

A method for working with a baseline-free algorithm while preserving the robustness of the calibration has been presented. This is done by using the absolute estimator for ray-tomography and an autocalibration method which is automatic, and is thus adapted to the SHM context. The absolute estimator does not compensate some systematic experimental uncertainties, which may result in a noisy map. This could make HARBUT diverge locally if this noise is too high. Moreover, it is crucial to have a clean map for the autocalibration step. In order to address this serious issue, a regularization method called variable relaxation has been developed.

Finally, because of the very low number of sensors used (four times less than the optimal number needed to reach the theoretical resolution of HARBUT) several regularizations methods have been presented. The Gaussian blur, which smoothly filters the spatial frequency components, is adjusted to the number of used sensors by tuning its standard deviation. The threshold, based on a physical constraint, allows to stabilize the reconstruction at each iteration. Lastly, the variable relaxation allows to compensate for the lack of sensors by reducing the effect of aliasing as well as the potential noise inherent to real data.

This paper has shown that it is possible to perform baseline-free HARBUT with a very low number of sensors compared to the optimal number, but this is at the expense of resolution. Future work could be performed to maintain a relatively good resolution even with a low number of sensors. Compressed sensing seems to be a promising tool for that purpose [START_REF] Wang | Méthodes de reconstruction d'images à partir d'un faible nombre de projections en tomographie par rayons x[END_REF]. Moreover, recently Shi and Huthwaite worked on the undersampled diffraction tomography problematic in the case of a 50 % reduction of the number of sensors. The idea is to use a forward model to get virtual transducers measurements in an iterative approach and so improve the imaging resolution [START_REF] Shi | Ultrasonic wave-speed diffraction tomography with undersampled data using virtual transducers[END_REF].

Another approach would be to use less intrusive sensors. For instance, Fiber Bragg Gratings (FBGs) is a promising solution to increase the number of measurement points without increasing the intrusiveness of the system. However, unlike PZT transducers generally used in SHM, FBGs cannot emit elastic waves. This shortcoming can be circumvented by employing passive methods in order to retrieve the Green's function from elastic diffuse fields -naturally present in structures -measured simultaneously between two sensors [START_REF] Druet | Passive guided waves measurements using fiber Bragg gratings sensors[END_REF][START_REF] Druet | Passive guided wave tomography for structural health monitoring applications[END_REF][START_REF] Druet | Feasibility of passive SHM for corrosion detection by guided wave tomography[END_REF].

Work is being conducted by the authors in parallel on passive tomography and will be discussed in a future paper. A longer-term program consists in performing passive tomography with FBGs only, in order to use a sufficient number of sensors to reach the maximum achievable resolution by the HARBUT algorithm.

  50

  Full ray distribution usually used on baseline data in order to calibrate the field in classic calibration. All the rays used in order to calibrate the field with the autocalibration method presented un this paper.

Fig. 3 .

 3 Fig. 3. Rays crossing the defect area after straight-ray tomography reconstruction.

Fig. 4 .

 4 Fig. 4. Calibration factors plotted in the complex plane.

  mm 600 mm ≈ 136 mm ≈ 4λ w and the HARBUT theoretical resolution λ Defect 1 -Only refraction phenomenon. e def = 0.7 mm, r bot = 0 cm and l trans = 25 cm (see Fig.8and Eq. (13) for the definitions).

  Defect 2 -Weak diffraction. e def = 0.3 mm, r bot = 5.6 cm and l trans = 6.4 cm (see Fig.8and Eq. (13) for the definitions).

  Defect 3 -Strong diffraction. e def = 0.65 mm, r bot = 3.5 cm and l trans = 3 cm (see Fig.8and Eq. (13) for the definitions).

Fig. 5 .

 5 Fig. 5. Forward propagation (ray-tracing) for different kinds of simulated defects.

γ 1 4 Fig. 6 .

 146 Fig. 6. Explanation of the regularization called "variable relaxation" for the case of time-of-flight tomography with ẽ = 2 mm and different values of β, γ and α.

Fig. 7 .

 7 Fig. 7. General FEM setup.

Fig. 8 .

 8 Fig. 8. Defect FEM setup.

  Frequency filtering for a PZT transducer (∅ PZT = 18 mm) coupled to a 2 mm thickness alu-Comparison between experimental measurements and signals from the simulation.

Fig. 9 .

 9 Fig. 9. Validation of the simulation.

Fig. 10 .

 10 Fig. 10. Ray-tomography with refraction only (defect 1).

Fig. 11 .

 11 Fig. 11. Ray-tomography with strong diffraction (defect 3).

  HARBUT behavior for 60 sensors and with a Gaussian blur of standard deviation 1λ. HARBUT behavior for 60 sensors, with a Gaussian blur of standard deviation 1 2 λ and a 2 mm threshold.

Fig. 12 .

 12 Fig. 12. HARBUT behavior (defect 2) as the sensor number decreases.

Fig. 13 .

 13 Fig.[START_REF] Belanger | Guided wave diffraction tomography within the Born approximation[END_REF]. Comparison between HARBUT with autocalibration and HARBUT with classic calibration (defect 3).

  (a) Aluminum plate and PZT transducer positions. (b) Defect position.

Fig. 14 .

 14 Fig. 14. Picture of the plate studied.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Reference measured with an immersion ultrasonic scan.

Fig. 17 .

 17 Fig. 17. Experimental straight-ray tomographies with a 2 mm threshold.

  Full ray distribution usually used on baseline data in order to calibrate the field for HARBUT. All the rays used in order to autocalibrate the field for baselinefree HARBUT.

Fig. 18 .

 18 Fig. 18. Rays crossing the defect area after straight-ray tomography reconstruction.

Fig. 19 .

 19 Fig. 19. Calibration factors plotted in the complex plane.

  Effect of the variable relaxation for the real part of O(x). Effect of the variable relaxation for the imaginary part of O(x).

Fig. 20 .

 20 Fig. 20. Effect of the variable relaxation for the DT step of HARBUT (β = 0.1, γβ = 1 and α = 2).

  No variable relaxation is used for the DT steps. Variable relaxation (β = 0.1, γβ = 1 and α = 2). Cross-section of Fig.21c.

Fig. 21 .

 21 Fig.[START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. Experimental HARBUT tomographies with a 2 mm threshold and using baseline data to calibrate the fields. The initial image used for iterative HARBUT is a straight-ray (relative estimator) tomography with a 2 mm threshold and with variable relaxation (β = 0.1, γβ = 1 and α = 4).

Fig. 22 .

 22 Fig.22. Baseline-free experimental HARBUT tomographies with autocalibration and a 2 mm threshold. The initial image used for iterative HARBUT is a straight-ray (absolute estimator) tomography with a 2 mm threshold and with variable relaxation (β = 0.1, γβ = 1 and α = 4).

Table 1 :

 1 Global relative error E global of Fig.13tomographies (data from simulation).

Table 2 :

 2 Global relative error E global of the tomographies of Fig. 21 and 22 (data obtained from the experience).


L Le Priol

A Dagmey

S Morandat

K Saleh

K El Kirat

A Nesterenko
email: alla.nesterenko@utc.fr

Comparative study of plant protein extracts as wall materials for the improvement of the oxidative stability of sunflower oil by microencapsulation

Keywords: oxidative stability, microencapsulation, vegetable oil, plant proteins, spray drying

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

Introduction

Polyunsaturated fatty acids (PUFAs), also called essential fatty acids, are divided in two categories, Ω-3 and Ω-6 PUFAs, depending on the position of the first double bond on their carbon chain. PUFAs, which are only produced by plants and phyto-planktons, are essential to enable normal growth and maintain good health of all higher organisms, including mammals and fishes. Unfortunately, they cannot be synthesized and need to be provided by the diet [START_REF] Rustan | Fatty Acids: Structures and Properties[END_REF]. This diet should respect the proper ratio between Ω-3 and Ω-6 PUFAs to meet the nutritional needs [START_REF] Dunbar | Omega 3 and omega 6 fatty acids in human and animal health: An African perspective[END_REF]. Ω-3 and Ω-6 PUFAs can be integrated in cell membranes and released on demand to serve as precursors of eicosanoid molecules [START_REF] Larsson | Dietary long-chain n-3 fatty acids for the prevention of cancer: a review of potential mechanisms[END_REF]. Eicosanoids have different biological effects on blood pressure regulation, modulation of inflammation or even immune responses [START_REF] Deckelbaum | Dietary n-3 and n-6 fatty acids: are there 'bad' polyunsaturated fatty acids?[END_REF]. Due to the highly unsaturated nature of PUFAs, they are sensitive to oxidation and thermic degradations leading to the production of hydroperoxides and unpleasant flavors and smells. For many years, microencapsulation of oils in polymeric matrices has been used to protect them from oxidative degradation [START_REF] Lewandowski | Morphology and microencapsulation efficiency of foamed spray-dried sunflower oil[END_REF]. Microencapsulation by spray drying is a relatively inexpensive, fast model for this study as its fatty acid composition contains a large amount of Ω-6 PUFAs. The effect of protein solubility on microparticles properties was discussed. Indeed, the emulsion stability index (ESI), the droplet size distribution and the viscosity of the corresponding oil/water emulsions were evaluated. The microparticles were characterized in terms of oxidative stability, encapsulation efficiency, water activity, moisture content, as well as morphology. This study showed that pea and soybean protein extracts are suitable wall materials for the protection of sunflower oil by microencapsulation.

Materials & methods

2.1. Materials Sunflower oil was kindly donated by the SAS PIVERT (Compiègne, France) and stored at room temperature. Commercial pea protein isolate (75% w/w protein), soybean protein isolate (90% w/w protein), brown rice protein (78% w/w protein), hemp protein (54% w/w protein) (MyProtein, UK) and sunflower protein (55% w/w protein) (SaludViva, ES) used as wall materials were purchased online and stored at room temperature. All others chemical were of analytical grade.

2.2. Protein characterizations 2.2.1. Solubility The protein solubility was measured according to the method described by [START_REF] Guimarães | Study of the proteins in the defatted flour and protein concentrate of baru nuts (Dipteryx alata Vog)[END_REF]. Protein suspensions were prepared at 1% w/v in distilled water. Protein solubility was tuned by modifying the pH in the range from 1 to 13 by adding NaOH or HCl. After stirring at room temperature for 1 h, the solutions were centrifuged at 6,000 rpm for 20 min at 20 °C (MR1812 centrifuge, Jouan, Saint-Nazaire, France). The supernatant was collected and the soluble protein content was analyzed using a Bradford test [START_REF] Bradford | A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding[END_REF] on an ultraviolet-visible spectrophotometer (Lamba 12, Perkin Elmer, San Jose, CA, US) at 595 nm. The solubility was defined as follows:

S % protein concentration in the supernatant initial protein concentration 100 2.2.2. Moisture content and water activity

The moisture content of the protein powders was measured gravimetrically. Briefly, 1 g of protein sample was put in an air oven at 120 °C for 6 h. The moisture content was calculated by weighting the sample before (W0) and after (W1) drying. The following equation was used:

Moisture % W -W W 100 The water activity was determined using a water activity meter (Aqualab 3TE instrument, Decagon, Pullman, WA, US). The sample was placed in a plastic cup to cover the entire surface. It was then placed in a closed chamber with temperature maintained at 25 ± 2 °C. Water activity measurements were performed after 10 min of sample equilibration in the instrument (Association of Official Analytical Chemists. & Cunniff, 1995).

Emulsion preparation

The protein powders used as wall materials were dispersed in distilled water to form 10% w/v solutions and mixed with a high-speed disperser (Ultra-Turrax T25, IKA-Labortechnik, Staufen, Germany) at 5,000 rpm for 5 min at room temperature to ensure protein hydration. The concentration of wall material was optimized at 10% w/v, based on a preliminary study conducted to determine the maximum concentration that can be incorporated in the solution. The pH values of the native protein solutions were 5.5, 6.3, 4.9, 7.8 and 6.3 for brown rice, hemp, pea, soybean and sunflower protein extracts, respectively. The pH values of these solutions were adjusted to 7.8 with 0.1 M NaOH. This pH value was chosen to ensure a proper protein solubility and remain in an acceptable pH range for food application. The emulsion was prepared by adding 10% w/v of sunflower oil (core/wall materials ratio 1:1) and mixed again with the high-speed disperser at 10,000 rpm for 5 min. The premixed emulsion was then stabilized by passing through a high pressure homogenization (HPH) device (Panda Plus 2000, GEA Niro Soavi, Parma, Italy) operated at 400 bars for two passes. In order to evaluate the impact of the optimization treatments, emulsions have also been prepared by keeping the initial pH values of the native protein solutions and by UT emulsification. The stability of these emulsions has been evaluated and compared to the emulsions prepared with the optimization treatments (i.e. pH adjustment and HPH). Emulsions formulated with brown rice, hemp, pea, soybean and sunflower protein extracts are named E-BR, E-HE, E-PE, E-SO and E-SU, respectively.

2.4. Emulsion characterizations 2.4.1. Emulsion stability Immediately after the emulsion preparation, 25 mL of emulsion was poured in a graduated cylinder for 24 h at 20 ± 5 °C to measure the emulsion stability index (ESI) [START_REF] Sarkar | Esterification of guar gum hydrolysate and gum Arabic with noctenyl succinic anhydride and oleic acid and its evaluation as wall material in microencapsulation[END_REF]:

ESI % 1 - V !"#$%$&"' #($!" V &)&$* "+,*!-). 100
where Vseparated phase represents the volume of the separated phase and Vtotal emulsion represents the total volume of poured emulsion (25 mL). 2.4.2. Droplet size distribution The emulsion droplet size distribution was evaluated using a Malvern MasterSizer 2000 (Malvern Instruments Ltd, Malvern, Worcestershire, UK). Droplet size measurements are reported as the mean diameters for each population.

2.4.3. Morphology Emulsion morphology was observed by optical microscopy. A drop of obtained emulsion was dissolved in distilled water and placed on a glass plate. The sample was covered with a glass slide and observed in Leica DM2700M optical microscope (Leica Microsystems, Wetzlar, Germany).

2.4.4. Viscosity Emulsion viscosity was measured using a Physica MCR301 Rheometer (AntonPaar, Graz, Austria) at imposed shear rates of 0.1-100 s -1 . Measurements were made using stainless steel plate-plate geometry with a diameter of 50 mm and a gap of 1 mm. The apparent viscosity of emulsions was obtained at 100 s -1 shear rate.

2.5. Dry microparticles preparation The freshly homogenized emulsions were spray dried using a lab scale spray dryer Büchi, B-290 (Büchi Labortechnik, Flawil, Switzerland). The emulsions were fed into the main chamber with a peristaltic pump and the feed flow rate was controlled by the pump rotation speed. The applied air inlet temperature was 160 °C and the outlet temperature was measured at 90 ± 2 °C. The liquid flow rate was 9 mL/min. The aspirator rate was set at 100% for all drying processes. The powder samples were collected and weighed. The prepared microparticles were stored at 25 °C until further analysis.

Microparticles formulated with brown rice, hemp, pea, soybean and sunflower protein extracts are named M-BR, M-HE, M-PE, M-SO and M-SU, respectively.

2.6. Microparticles characterizations 2.6.1. Microparticles size distribution The size of the microparticles was measured at room temperature by laser diffraction using a Malvern Mastersizer 2000 equipment with Scirocco 2000 unit (Malvern Instruments Ltd, Malvern, Worcestershire, UK). Microparticles size measurements are reported as the mean diameter for each population.

2.6.2. Encapsulation efficiency The extraction of the microparticle surface oil was performed by following the method of [START_REF] Liu | Entrapment of Flaxseed Oil Within Gelatin-Gum Arabic Capsules[END_REF] with modifications. A Fisherbrand™ Porcelain Buchner Funnel was covered by a Whatman filter paper No. 1. Dry microparticles (1 ± 0.001 g) were weighted and placed on the filter paper. The microparticles were rinced three times with 6 mL of hexane. The organic phase was evaporated until constant weight to access complete solvent removal. Protein-to-oil ratio in microparticles was presumed to be equal to the protein-to-oil ratio in the emulsions. The encapsulation efficiency (EE) was calculated with the following equation:

EE % P &)&$* )-* -P !,%0$1" )-* P &)&$* )-* 100
where Psurface oil represents the percent ratio of oil content on the surface of the microparticles and Ptotal oil represents the percent ratio of oil content in dry matter of initial emulsion.

2.6.3. Moisture content and water activity The moisture content and the water activity of the microparticles were measured as described in section 2.2.2. Moisture content was calculated by the following equation:

Moisture % 2 W +-1%)#$%&-1*"! -W '%3 +-1%)#$%&-1*"! W +-1%)#$%&-1*"! 4 100
8 where Wmicroparticles represents the weight of original dry sample and Wdry microparticles represents the weight of sample after oven treatment.

2.6.4. Sunflower oil oxidative stability in accelerated storage test Accelerated oxidation tests were carried out on the crude oil and on the dry microparticles using a Rancimat apparatus (743 Rancimat METROHM, Switzerland). The samples were exposed to high temperature (100 °C) in Rancimat tubes. A stream of purified air at a flow rate of 10 L/h was injected inside the tubes to promote oxidation. The volatile oxidation products released in the atmosphere of tubes were carried by air to containers filled with water for conductivity measurements, thanks to electrodes connected to a measuring and recording device. The increase in conductivity was related to the oxidative stability of the oil. Based on this, the induction period (also called induction time) was calculated. It is defined as the time corresponding to the inflection point of the conductivity versus time curve (when the conductivity of water begins to increase rapidly). The higher the induction time, the more stable the oil.

2.6.5. Microparticles morphology The morphology of the spray dried microparticles was observed with an environmental scanning electron microscope (ESEM, Quanta 250 FEG, FEI Co., OR, USA). Samples were prepared by mounting the powders on an aluminum stub using an adhesive carbon tape. Samples were then coated with gold used as sputter coating and imaged at 20 kV accelerating voltage. Micrographs were taken at different magnifications in order to visualize the surface morphology of the microparticles.

2.7. Statistical analysis All of the characterization measurements of protein extracts, emulsions and microparticles were performed in triplicate. Results were expressed as the mean ± standard deviation and were statistically calculated using an analysis of variance (ANOVA). Comparison of means were performed by Tukey analyses at p<0.05.

Results and discussion

Solubility of protein powders

Thanks to their amphiphilic character, protein molecules have good emulsifying properties. When the protein extract has a good water solubility, a large amount of protein chains is able to diffuse to the oil/water interface and stabilize small droplets of emulsion. When this emulsion is transformed into an aerosol during the spray drying process, the small oil droplets are well distributed inside the spraying drops and are efficiently encapsulated inside the particles during the drying step. The powder is thus made of individual microparticles with a low surface oil content. On the contrary, if the protein extract has a poor water solubility, a small amount of protein chains is able to stabilize the emulsion and the oil droplets will be larger. During the spraying, the dispersed phase is not well distributed and a large amount of oil stays on the surface of the particles. This poor encapsulation will lead to the agglomeration of the microparticles (Fig. 1). The analysis of the solubility profiles of protein extracts is then mandatory to evaluate the effect of the wall material solubility on the microencapsulation process. In this study, the effect of the pH on protein solubility was studied (Fig. 2). The plant proteins had a U-shaped-like solubility profile, which is consistent with the literature [START_REF] Tang | Physicochemical and Functional Properties of Hemp ( Cannabis sativa L.) Protein Isolate[END_REF][START_REF] Tömösközi | Isolation and study of the functional properties of pea proteins[END_REF][START_REF] Withana-Gamage | Physicochemical, thermal and functional characterisation of protein isolates from Kabuli and Desi chickpea (Cicer arietinum L.): a comparative study with soy (Glycine max) and pea (Pisum sativum L.)[END_REF]. The lowest solubility observed at pH 4-5 corresponds to the isoelectric point of the proteins. When moving away from this point, the polar groups on protein chains are charged and the solubility increases. The differences between proteins solubilities can be explained by their composition in protein fractions. The protein compositions of the five extracts could then give an indication on the best candidates for the emulsification and the protection of sunflower oil by spray drying encapsulation. Plant seed proteins are divided into four main fractions, differentiated by their solubility properties: the albumin fraction, soluble in water; the globulin fraction, soluble in dilute saline solutions; the prolamin fraction, soluble in hydroalcoholic solvent (60-70% v/v) and the glutelin fraction, soluble in very alkaline water solutions (pH>10) [START_REF] Osborne | The vegetable proteins[END_REF]. In seeds, half or more of the total proteins are storage proteins. The major role of these proteins is to provide a store of nutrients for the plant growth [START_REF] Kawakatsu | Proteins[END_REF]. The composition of plant proteins used as wall materials in this study is presented in Table 1. Globulins are the major storage protein fractions in dicotyledonous plants seeds (pulses, legumes and oilseeds), whereas prolamins and glutelins are the main fractions (80-90%) in cereals [START_REF] Guéguen | Les protéines végétales : contexte et potentiels en alimentation humaine[END_REF]. This characteristic could explain the lower solubility of brown rice proteins compared to the other proteins (at pH 12, all protein extracts have a solubility percentage close to 100%, except brown rice proteins). On the other hand, the different solubility profiles observed among the other proteins can be explained by their proportion of soluble fractions. At pH 7.8, hemp, pea and soybean proteins appeared to be the most soluble macromolecules. According to the results obtained by [START_REF] Tang | Physicochemical and Functional Properties of Hemp ( Cannabis sativa L.) Protein Isolate[END_REF], hemp proteins contain almost exclusively soluble protein fractions (87% of globulins divided into 82% of edestin and 5% of vicilin, plus 13% of albumins). Our result showing that hemp protein extract is the most soluble at the pH of our study is therefore in accordance with this report. In addition to the content in soluble fraction, the protein fraction has an importance on the emulsifying properties of protein extracts. As it can be seen in Table 1, pea and soybean proteins contain a large amount of vicilins (also called 7S globulins) compared to the other plant proteins. According to [START_REF] Chen | Study on the emulsifying stability and interfacial adsorption of pea proteins[END_REF], pea proteins contain around 30% w/w of vicilins. Soy proteins contain more than 80% globulins with a ratio 7S/11S of 0.5-1.3 depending on the variety [START_REF] Nishinari | Soy proteins: A review on composition, aggregation and emulsification[END_REF]. It had been showed that 7S globulin has better emulsifying properties than 11S globulin, due to its lesser size and higher flexibility [START_REF] Chen | Study on the emulsifying stability and interfacial adsorption of pea proteins[END_REF][START_REF] Dagorn-Scaviner | Emulsifying Properties of Pea Globulins as Related to Their Adsorption Behaviors[END_REF]. Indeed, the solubility results and the protein fraction compositions suggest that soy and pea protein extracts could be good candidates for the microencapsulation of sunflower oil. 3.2. Properties of emulsions 3.2.1. Emulsion stability index Obtaining a stable liquid emulsion is a prerequisite for proper encapsulation with spray drying and the ESI mainly depends on the emulsifying properties of wall materials and homogenization technique [START_REF] Pinnamaneni | Comparison of oil-in-water emulsions manufactured by microfluidization and homogenization[END_REF]. The composition of a kinetically unstable oil/water emulsion is presented in Fig. 3. The top phase is the creaming phase. It appears when the dispersed phase of kinetically unstable emulsions migrates under the influence of buoyancy. The middle phase corresponds to the remaining emulsion. The bottom phase is composed of sedimented proteins. When the emulsification conditions do not allow to properly solubilize the proteins, they gradually migrate at the bottom of the cylinder under the influence of gravity. Since the effects of buoyancy and gravity also depend on the viscosity of the solution, the ESI values do not perfectly coincide with the water solubility but also with the viscosity of the protein extracts in solution. For example, the perfect stability of E-SO after 24 h of rest could be explained by the high viscosity of the emulsion limiting the migration of oil droplets (Fig. 4 and Table 2). Moreover, it could also be explained by the higher protein content of soy extract (90% w/w) compared to other protein extracts (54-78% w/w). The key point here is that after pH adjustment at 7.8 and HPH treatment, the emulsions were kinetically stable after 24 h of storage at room temperature, with unseparated phase fractions of 92-100%, showing the efficiency of the selected conditions to stabilize the emulsions.

Emulsion droplet size & viscosity

The droplet size is an important parameter for the stability of the emulsions. Small droplets are less affected by destabilization phenomena. The droplet size distributions before and after the optimization treatments are presented in Fig. 5. The emulsification by HPH significantly reduced the droplet size but also induced polydispersed distributions. The type of wall material had a significant effect on the droplet size since the mean diameter of the majority populations ranged from 0.2 to 22.6 µm after pH modification and HPH treatment (Table 2). The treatment was particularly effective on emulsions stabilized by soybean and pea protein extracts. Indeed, the mean diameter of the majority population decreased from 91.2 to 0.3 µm and from 66.0 to 0.2 µm, respectively. These results are in accordance with the solubility profiles of the protein extracts. Soybean and pea proteins were the most soluble macromolecules, with hemp proteins, at pH 7.8. When proteins are well solubilized in water, a large amount of proteins is able to diffuse to the oil/water interface. A large surface area can thus be stabilized and the emulsion droplets are smaller [START_REF] Hoffmann | Emulsions with unique properties from proteins as emulsifiers[END_REF]. This observation was also made by [START_REF] Wang | High pressure homogenization combined with pH shift treatment: A process to produce physically and oxidatively stable hemp milk[END_REF] in their study of HPH and pH shift treatments on hemp milk stability. In fact, when the pH is shifted from the isoelectric point of the protein, it increases the amount of electrostatic charges and thereby increasing the protein solubility. Moreover, the intense mechanical forces exerted on protein chains during HPH treatment may cause an increase of their flexibilities. When the protein macromolecules will reach the oil/water interface, they will be able to unfold and expose their hydrophobic regions at the interface and stabilize the emulsion by coating the interface [START_REF] Cabra | Emulsifying properties of proteins[END_REF]. Proteins can then provide physical stability to the emulsion [START_REF] Jiang | Interfacial Structural Role of pH-Shifting Processed Pea Protein in the Oxidative Stability of Oil/Water Emulsions[END_REF][START_REF] Nesterenko | Influence of soy protein's structural modifications on their microencapsulation properties: α-Tocopherol microparticle preparation[END_REF]. Concerning emulsions stabilized with hemp protein extracts, the relatively high mean diameter of the majority population (7.6 µm) after optimization treatments could be explained by the lower amount of protein in the extract compared to soybean and pea protein extracts (only 54% w/w of protein). Non-soluble residus of protein agglomerates contained in hemp protein emulsion (E-HE) are visible in Fig. 6 (yellow arrow), corresponding to the population 2 in Table 2. For E-BR and E-SU, the pH modification and the HPH treatment significantly decreased the droplet size but it also led to the formations of oil droplets aggregates (see red arrows in Fig. 6). Actually, the protein chains unfolding will expose non-polar regions and, thus, will increase the total surface hydrophobicity of the protein.

The following HPH treatment will then provide enough energy to the system to make hydrophobic groups interactions possible and lead to aggregation [START_REF] Lee | Effects of ultra-high pressure homogenization on the properties and structure of interfacial protein layer in whey proteinstabilized emulsion[END_REF]. Since brown rice and sunflower protein extracts are the less soluble materials (Fig. 1), it can be assumed that their chains contain more hydrophobic groups than other protein extracts, which could be related to the aggregations visible in Fig. 6. The viscosity is also an important parameter for solutions intended to be spray-dried. Low viscosities insure a proper formation of the aerosol and an efficient drying of the droplets. Usually, an adequate spraying is insured if the viscosity does not exceed 300 mPa.s [START_REF] Di Battista | The use of arabic gum, maltodextrin and surfactants in the microencapsulation of phytosterols by spray drying[END_REF]. The viscosities of the emulsions are shown in Table 2. Two behaviors are visible, depending on the nature of the plan protein. For non-hydrocolloid macromolecules (brown rice, hemp, pea and sunflower proteins), the HPH treatment induced an increase of the apparent viscosity. Indeed, the intense mechanical forces provided by the HPH lead to a structure modification of the proteins with the unfolding of proteinic chains, resulting in a better solubility and an increase of the viscosity. Owing to the structural properties of soybean proteins, the behavior of E-SO was different [START_REF] Hu | Gelation of soybean protein and polysaccharides delays digestion[END_REF]. Before optimization treatments, the viscosity of E-SO was significantly higher than the others. Indeed, soybean proteins have the ability to form gels with good holding capacity and they are often used by the food industry due to this functional property [START_REF] Utsumi | Heat-induced interactions between soybean proteins: preferential association of 11S basic subunits and .beta. subunits of 7S[END_REF]. The emulsification optimization treatment significantly decreased the viscosity of E-SO (from 210.6 to 3.3 mPa.s). [START_REF] Song | Effect of high-pressure homogenization on particle size and film properties of soy protein isolate[END_REF] also noticed that the viscosity of their non-homogenized soy protein isolate (SPI) suspension was much higher than homogenized samples. As said before, HPH provides extra energy in the system and this energy provokes the disruption of non-covalent interaction forces (electrostatic, Van der Waals and hydrophobic) possible between SPI molecules. This may lead to the degradation of the network structure of SPI gel and decrease the viscosity of the solution. Furthermore, the decrease of viscosity is an important parameter in emulsion characterization as it may affect the size of the spray-dried microparticles. The spray drying of high viscosity emulsions would form larger microparticles, due to the increase of solid content in each drop [START_REF] Nesterenko | Vegetable proteins in microencapsulation: A review of recent interventions and their effectiveness[END_REF][START_REF] Patel | Revealing facts behind spray dried solid dispersion technology used for solubility enhancement[END_REF][START_REF] Tonon | Influence of emulsion composition and inlet air temperature on the microencapsulation of flaxseed oil by spray drying[END_REF]. Table 2 also shows that, after optimization treatments, the viscosity was positively correlated to the droplet size of the emulsion. E-PE and E-SO showed the lowest viscosity, followed by E-HE, E-SU and E-BR (Table 2). This relation was also noticed by [START_REF] Turchiuli | Production of fine emulsions at pilot scale for oil compounds encapsulation[END_REF]. [START_REF] Tatar | Rheology of Emulsions[END_REF] described that, in concentrated solutions, smaller droplets have the ability to pack more efficiently than larger droplets, decreasing the emulsion viscosity. On the other side, the high viscosities of E-BR and E-SU could be explained by the formation of the aggregates after HPH treatment, corresponding to populations 2 in Table 2 and the red arrows visible in Fig. 6. Their presence tends to increase the resistance to flow, increasing the apparent viscosity of the emulsion. The different viscosities observed for E-BR and E-SU after optimization treatments could be explained by the natures of the extracts (constituents of the extracts and proteins composition) and by the distribution size of the emulsion droplets (higher volume proportion of population 2 of E-BR compared to the population 2 of E-SU).

3.3. Properties of microparticles 3.3.1. Moisture content and water activity Moisture content, which measures the total amount of water in a compound, is a critical parameter for formed microparticles. At high moisture contents, the properties of the wall materials change [START_REF] Velasco | Variables affecting lipid oxidation in dried microencapsulated oils[END_REF]. This change will induce stickiness of powder particles, resulting in the formation of inter-particles bridges that lead to caking, particle collapse and the release and oxidation of the core material during storage [START_REF] Beristain | Effect of Water Activity on the Stability to Oxidation of Spray-Dried Encapsulated Orange Peel Oil Using Mesquite Gum (Prosopis Juliflora) as Wall Material[END_REF][START_REF] Drusch | Physicochemical characterization and oxidative stability of fish oil encapsulated in an amorphous matrix containing trehalose[END_REF][START_REF] Harnkarnsujarit | Glass-Transition and Non-equilibrium States of Edible Films and Barriers[END_REF][START_REF] Partanen | Effect of Relative Humidity on Oxidation of Flaxseed Oil in Spray Dried Whey Protein Emulsions[END_REF]. The moisture contents of the plant proteins extracts and microparticles are presented in Table 3.

The residual water content ranged from 1.5 to 2.3% w/w for all five microparticles. Actually, there is no specification for the moisture content of dry food formulations made with plant proteins but these results comply with standard moisture content accepted for spray dried dairy powders (<4% w/w) [START_REF] Schuck | Relative humidity of outlet air: the key parameter to optimize moisture content and water activity of dairy powders[END_REF]. Since the process variables were held constant for all the experiments, the small difference in moisture content observed between microparticles is related to the affinity of the material for water and water diffusivity through polymer matrix. [START_REF] Bajaj | Pea Protein Isolates: Novel Wall Materials for Microencapsulating Flaxseed Oil[END_REF] obtained water contents ranging from 3.9 to 4.25% for their microparticles made of pea protein isolates and flaxseed oil and spray dried at 150 °C, which indicates the influence of the inlet temperature on the moisture content of the final product. Moisture content alone is not a sufficient indicator of food powder stability, since foods with the same water content do not necessarily have the same perishability [START_REF] Nielsen | Introduction to Food Analysis[END_REF]. The water activity (aw) also describes water content in powders but it provides information on how the water associates with other constituents. For example, when water is bounded to proteins, it is less available for chemical reactions and microbial growth. The aw of the protein extracts and microparticles formulated with the different wall materials are shown in Table 3. It was found that all microparticles had lower aw than the protein extracts used to produce them. These low aw (ranging from 0.118 to 0.269) could be attributed to the effective water evaporation during spray drying. Generally, foods with a aw<0.6 are considered microbiologically stable [START_REF] Quek | The physicochemical properties of spray-dried watermelon powders[END_REF].

3.3.2. Particle size and morphology Fig. 7 shows the particle size distribution of powders produced with the different wall materials. M-PE and M-SO showed monodispersed distributions, with d4,3 of 13.4 ± 3.4 and 16.3 ± 3.2 µm, respectively. These first peaks, around 10 µm, are mainly determined by the spray-dryer nozzle diameter used [START_REF] Di Giorgio | Encapsulation of fish oil in soybean protein particles by emulsification and spray drying[END_REF]. The other populations visible for M-BR, M-HE and M-SU correspond to agglomerates. The poor solubility of brown rice and sunflower protein extracts and the low protein amount in hemp protein extracts did not allowed to efficiently encapsulate sunflower oil and, thus, lead to high oil amount on the surface of the microparticles, leading to their agglomeration. Moreover, the contact between particles induced by this high oil surface content can lead to the formation of inter-particles bridges with the production of large size agglomerates [START_REF] Tonon | Influence of emulsion composition and inlet air temperature on the microencapsulation of flaxseed oil by spray drying[END_REF]. ESEM images of the microparticles prepared with the different wall materials are presented in Fig. 8. They revealed important differences in microparticles shapes and surface regularities. The agglomeration of M-BR, M-HE and M-SU is also visible on these images. Clusters with rough surfaces are visible for M-BR. Individual particles cannot be clearly distinguished in this sample. M-SU showed a smoother surface than M-BR but the microparticles are also highly agglomerated. Pores are visible on the surface of M-HE. On the contrary, M-PE and M-SO exhibited individual particles with smooth surfaces and no apparent cracks or fissures, which is important to ensure a lower gas permeability and a better oil protection. The shrinkage observed for these samples is typical from microparticles produced with proteins as wall materials [START_REF] Tonon | Influence of emulsion composition and inlet air temperature on the microencapsulation of flaxseed oil by spray drying[END_REF]. [START_REF] Gong | Emulsifying properties and structure changes of spray and freeze-dried peanut protein isolate[END_REF] suggested that this phenomenon may be due to the rapid formation of a dried crust layer on the surface of the microparticles followed by the high flux of moisture leaving the particle during drying. [START_REF] Tang | Microencapsulation properties of soy protein isolate and storage stability of the correspondingly spray-dried emulsions[END_REF] also noticed the shriveling for their microparticles made of soy protein isolate and soy oil. [START_REF] Xu | Effects of Emulsification of Fat on the Surface Tension of Protein Solutions and Surface Properties of the Resultant Spray-Dried Particles[END_REF] observed the similar characteristic for their microparticles containing sunflower oil protected by whey protein isolate and maltodextrin. The internal morphology of obtained M-PE microparticles showed the "sponge-like" structure of the protein matrix inside which oil droplets are located (red arrow in Fig. 8).

3.3.3. Encapsulation efficiency and oxidative stability In order to evaluate the protective efficiency of the wall materials on the oxidative stability of sunflower oil, accelerated oxidation tests were conducted with a Rancimat apparatus by comparing induction periods (IP) of microparticles with pure sunflower oil (Table 4). Data obtained from these experiments showed that pure sunflower oil (control) had an IP of 9.50 h. IPs of M-SO and M-PE were significantly higher than that found for the non-encapsulated oil. They presented IPs values of 12.49 and 21.26 h, respectively. This result is in agreement with previous observations on the solubility of protein extracts, the droplet size distributions of the emulsions and the morphologies of the microparticles. In comparison, [START_REF] Ahn | Antioxidant effect of natural plant extracts on the microencapsulated high oleic sunflower oil[END_REF] obtained microparticles of sunflower oil supplemented in natural plant extracts as antioxidants and protected by dextrin-MPI wall materials with an IP value of 16.26 h, which is lower than the IP value obtained for M-PE without the addition of antioxidant. The fact that particles stabilized by hemp proteins did not lead to a significant improvement of sunflower oil oxidative stability (IP value of 9.72 h) could be explained by the agglomeration and the porous nature of the microparticles, which facilitates the permeation of gas, moisture and oil release. Microparticles formulated with sunflower and brown rice protein extracts showed significantly lower IP compared to non-encapsulated sunflower oil. This could be directly related to the previous results, particularly on the observation of the external structure of the microparticles. These results showed a direct correlation with the values of the EE (Table 4). M-SU and M-BR had significantly lower EE values (79 and 69 %, respectively) compared to M-PE, M-SO and M-HE (93, 91 and 89 %, respectively). The proper retention and more efficient encapsulation of sunflower oil inside M-PE, M-SO and M-HE could then explain their higher IPs values.

Conclusion

This study intended to provide a comparative analysis of five plant protein extracts for the microencapsulation of sunflower oil with the aim of improving its oxidative stability. The microparticles prepared by the spray drying technique exhibited low moisture contents and low water activities. Our results demonstrate that the nature of plant protein extracts used for the microencapsulation of sunflower oil strongly affects the oxidative stability efficiency. A summary of the efficiency of the protein extracts on different parameters characterized during this study is presented in Table 5. Soybean and pea protein extracts are suitable wall materials for the encapsulation and the protection of sunflower oil. The microencapsulation remarkably improves the sunflower oil oxidative stability and allowed to multiply by 2.2 times the IP of the control for the microencapsulation by pea protein extracts. The improvement seems to be related to the solubility of the protein extracts and their structural properties. These findings are of importance for providing a solution to develop PUFA-enriched formulations for food and feed industries.

Process steps

Table 2

 2 Volume distribution size and viscosities of sunflower oil/water emulsions stabilized with brown rice (E-BR), hemp (E-HE), pea (E-PE), soybean (E-SO) and sunflower (E-SU) protein extracts. -d means in each column followed by different letters were significantly different (p < 0.05)

	Samples Droplet distribution size (µm)			Apparent viscosity at 100 s -
						1 (mPa.s)
		UT *		UT + HPH **		UT	UT + HPH
		Population 1 Population 2 Population 1 Population 2
	E-BR	47.3 ± 0.1 b	-	0.5 ± 0.1 b	22.6 ± 0.6 c	10.2 ± 0.1 d	25.0 ± 1.4 d
	E-HE	23.4 ± 0.2 a	126.6 ± 0.8	0.5 ± 0.1 b	7.6 ± 1.1 b	2.0 ± 0.1 b	6.1 ± 0.3 b
	E-PE	66.0 ± 1.0 c	-	0.2 ± 0.1 a	3.2 ± 0.2 a	1.3 ± 0.1 a	3.3 ± 0.2 a
	E-SO	91.2 ± 13.5 d -	0.3 ± 0.1 ab	-	210.6 ± 24.0 e	3.3 ± 0.2 a
	E-SU	67.9 ± 1.0 c	-	7.9 ± 0.1 c	67.6 ± 0.6 d	8.6 ± 0.8 c	12.3 ± 0.4 c

a* UT: emulsion had been pre-homogenized with an Ultra-Turrax ** UT+HPH: emulsion had been homogenized with an Ultra-Turrax followed by high pressure homogenization treatment

Table 3

 3 Moisture contents and water activities of the protein extracts and respective microparticles formulated with brown rice (BR and M-BR), hemp (HE and M-HE), pea (PE and M-PE), soybean (SO and M-SO) and sunflower (SU and M-SU).

	Samples	Moisture (%)	Water activity
	BR	5.2 ± 0.8 d	0.240 ± 0.001 d
	M-BR	2.3 ± 0.1 c	0.142 ± 0.011 ab
	HE	8.2 ± 0.8 e	0.440 ± 0.076 g
	M-HE	2.1 ± 0.3 bc	0.202 ± 0.022 c
	PE	5.8 ± 0.7 d	0.316 ± 0.003 f
	M-PE	1.4 ± 0.1 a	0.118 ± 0.019 a
	SO	7.2 ± 0.7 e	0.272 ± 0.005 e
	M-SO	1.6 ± 0.1 a	0.269 ± 0.008 e
	SU	7.5 ± 0.5 e	0.358 ± 0.004 g
	M-SU	2.0 ± 0.1 b	0.175 ± 0.028 bc

a-g means in each column followed by different letters were significantly different (p < 0.05)

Table 4

 4 Induction period and encapsulation efficiency (EE) for microparticles formulated with pea (M-PE), soybean (M-SO), hemp (M-HE), sunflower (B-SU), brown rice (M-BR) protein extracts and sunflower oil as control. -d means in each column followed by different letters were significantly different (p < 0.05)

	Samples	Induction period (h) EE (%)
	Pure sunflower oil (control)	9.50 ± 0.10 b	-
	M-PE	21.26 ± 0.44 d	88 ± 2 b
	M-SO	12.49 ± 0.38 c	91 ± 1 b
	M-HE	9.72 ± 0.13 b	89 ± 3 b
	M-SU	7.20 ± 0.28 a	79 ± 4 a
	M-BR	6.97 ± 0.07 a	69 ± 7 a

a

Acknowledgements

This work has been performed, in partnership with the SAS PIVERT, within the frame of the French Institute for the Energy Transition (Institut pour la Transition Energétique (ITE)) P.I.V.E.R.T. (www.institut-pivert.com) selected as an Investments for the Future (Investissements d'Avenir). This work was supported, as part of the Investments for the Future, by the French Government under the reference ANR-001-01.

Table 1

Fraction composition of proteins extracted from brown rice, soybean, pea, sunflower and hemp seeds, adapted from the works of [START_REF] Kawakatsu | Proteins[END_REF].

Species