Chady El Hachem

Kamilia Abahri

Sébastien Leclerc

Rachid Bennacer
email: rachid.bennacer@ens-paris-saclay.fr

Sebastien Leclerc

NMR and XRD quantification of bound and free water interaction of spruce wood fiber

Keywords: Bound water, free water, fiber saturation point, swelling, spruce wood, Nuclear Magnetic Resonance, X-Ray Diffraction

Introduction

Because of its thermal insulation properties and energetic performances, wood is being widely used in the construction field [START_REF] Dodoo | Effect of thermal mass on life cycle primary energy balances of a concrete-and a wood-frame building[END_REF]. Spruce wood, the studied material in this paper, is particularly known as well for its mechanical, thermal and hygric properties. This material presents an anisotropic microstructure, characterized by three principle directions: longitudinal, radial and tangential. However, this anisotropy is accompanied by a multi-scale behavior. In fact, at the microscopic scale, spruce wood is constituted of fibers (tracheids and parenchyma). Each fiber is formed of cell walls, containing porosity holes [START_REF] Mäkinen | Predicting wood and tracheid properties of Norway spruce[END_REF]Rafsanjani et al., 2012). This multi-scale behavior is as well observed at the cell walls' chemical constituents of wood, and biomaterials in general [START_REF] Salvati | Multiscale analysis of bamboo deformation mechanisms following NaOH treatment using Xray and correlative microscopy[END_REF]. In fact, the cell walls are constituted of many layers, in which mainly cellulose (crystalline and amorphous), hemicellulose and lignin coexist. Each of these constituents is characterized by its own thermal, hygric and mechanical properties, which strongly affect the macroscopic behavior of spruce wood [START_REF] Derome | The role of water in the behavior of wood[END_REF][START_REF] Almeida | Free shrinkage of wood determined at the cellular level using an environmental scanning electron microscope[END_REF][START_REF] Cueff | Experimental and numerical study of the thermomechanical behaviour of wood-based panels exposed to fire[END_REF].

On the other hand, wood is a hygroscopic bio-based material, often in interaction with humidity [START_REF] Gezici-Koç¸ | Bound and free water distribution in wood during water uptake and drying as measured by 1D magnetic resonance imaging[END_REF][START_REF] Zlahtic-Zupanc | Changes in moisture performance of wood after weathering[END_REF][START_REF] Zhou | Wetting enhanced by water adsorption in hygroscopic plantlike materials[END_REF]. Therefore, it tends to swell or shrink, whether its water content increases or decreases respectively [START_REF] Bonnet | NMR determination of sorption isotherms in earlywood and latewood of Douglas fir. Identification of bound water components related to their local environment[END_REF][START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF][START_REF] El Hachem | Original experimental and numerical approach for prediction of the microscopic hygro-mechanical behavior of spruce wood[END_REF]. In fact, the cell walls dimensional variations take place following the absorption (or desorption) of bound water, whereas free water does not affect the wood structure, because it occupies the space within the cavities (porosity). Along the sorption cycle, the material absorbs first the bound water until the cell walls get fully saturated [START_REF] Gezici-Koç¸ | Bound and free water distribution in wood during water uptake and drying as measured by 1D magnetic resonance imaging[END_REF]. This state of the material corresponds to the Fiber Saturation Point (FSP) [START_REF] Eitelberger | A comprehensive model for transient moisture transport in wood below the fiber saturation point: physical background, implementation and experimental validation[END_REF]. Once the FSP is achieved, there is no more swelling within the material, and every additional increase in the water content of the material is assimilated to free water within the porosity [START_REF] Gezici-Koç¸ | Bound and free water distribution in wood during water uptake and drying as measured by 1D magnetic resonance imaging[END_REF].

The microscopic dimensional variations caused by the relative humidity changes have been previously revealed [START_REF] Derome | The role of water in the behavior of wood[END_REF][START_REF] Almeida | Free shrinkage of wood determined at the cellular level using an environmental scanning electron microscope[END_REF][START_REF] Guo | Changes of wood cell walls in response to hygro-mechanical steam treatment[END_REF][START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF][START_REF] El Hachem | Hygromorphic characterization of softwood under high resolution X-ray tomography for hygrothermal simulation[END_REF]. In fact, at the cell walls' scale, there is an instability related to the swelling and shrinkage caused by the relative humidity variations, and particularly to the various chemical components properties [START_REF] Wang | Predictions of effective physical properties of complex multiphase materials[END_REF][START_REF] Neagu | Modelling of effects of ultrastructural morphology on the hygroelastic properties of wood fibres[END_REF][START_REF] Guo | Changes of wood cell walls in response to hygro-mechanical steam treatment[END_REF]. All these factors strongly influence the cell walls' response when subjected to various environmental conditions of relative humidity and temperature.

There is still lack of knowledge in literature concerning the interaction of the wood cell with relative humidity. There is no real distinction between the free and bound water with precise estimation when the FSP occurs.

In a previous work [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF], we investigated by X-ray tomography the pore size and the fibers dimensions' evolution function of the water content of wood. However, the X-ray tomography technique is limited to distinguish the water interaction within the wood microstructure.

Recently, Nuclear Magnetic Resonance (NMR) technique was used for wooden material [START_REF] Almeida | NMR imaging of fluid pathways during drainage of softwood in a pressure membrane chamber[END_REF][START_REF] Hernández | Magnetic resonance microimaging of liquid water distribution in sugar maple wood below fiber saturation point[END_REF]Fourmentin et al., 2016;[START_REF] Gezici-Koç¸ | Understanding the influence of wood as a substrate on the permeability of coatings by NMR imaging and wet-cup[END_REF][START_REF] Zhou | Magnetic resonance imaging evidences of the impact of water sorption on hardwood capillary imbibition dynamics[END_REF][START_REF] Rostom | Investigation of the effect of aging on wood hygroscopicity by 2D 1 H NMR relaxometry[END_REF]. For example, [START_REF] Almeida | NMR imaging of fluid pathways during drainage of softwood in a pressure membrane chamber[END_REF] studied the water restitution of initially saturated Douglas wood when subjected to dried conditions.

On the other hand, Fourmentin et al. ( 2016) studied the water adsorption/uptake of initially dried specimens of hemp when subjected to in-situ saturated conditions. The principal results of such a technique is the water content quantification.

Therefore, the present work aims to study the coupling between the hygric response of bound and free water at the microscopic scale and the local swelling related to the chemical components of the material, with respect to the FSP. In addition, the material in question is spruce wood, being a sapwood having a density of 450 Kg/m 3 , which has not been studied yet which such techniques (NMR and XRD).

The first part of the paper concerns the bound and free water quantification along the sorption cycle, using NMR. A special post-processing procedure and a specific device have been developed for this purpose. The relative humidity interactions with the material, at the scale of the cell walls' chemical constituents will be afterwards treated by X-Ray Diffraction (XRD) methods. Once more, a specific procedure has been set up in order to do the experiment considering each humidity state of the material. Finally, a global interpretation of the results will take place, in order to highlight the physical interactions at different scales.

Material and methods

An accurate study of bound and free water interactions with the solid matrix of wood requires very well organized and refined experiments' protocols. Both adopted techniques, Nuclear Magnetic Resonance (NMR) and X-ray Diffraction (XRD), as well as the corresponding experimental protocols and the specific devices used to maintain a constant relative humidity in each experiment will be first discussed.

Nuclear Magnetic Resonance (NMR)

NMR is a non-destructive experimental technique used in the characterization of materials. It consists in modifying the magnetization of nuclear spins of certain atoms, due to the effect of two magnetic fields B 0 and B 1 . The magnetic field B 0 presents a high and fixed intensity, whereas the magnetic field B 1 is rotational [START_REF] Keeler | Understanding NMR spectroscopy[END_REF]. The hydrogen atom is the most sensible nucleus [START_REF] Arends | Moisture transport in pine wood during one-sided heating studied by NMR[END_REF][START_REF] Wang | Penetration and distribution of paraffin wax in wood of loblolly pine and Scots pine studied by time domain NMR spectroscopy[END_REF] due to its ½ spin and high intrinsic gyromagnetic moment of 42.6 MHz/T. In addition, this atom is very abundant, which leads to a well-defined resonance phenomenon.

This technique is principally used in the medical field, and for organic molecules characterization purposes. Relaxometry methods were developed more recently for civil engineering applications. For example, this technique can be used to characterize the porous structure of materials, as well as the quantification of free water in the porosity and other water states that interact with the solid matrix of the material. In this work, only free water in the porosity, and bound water absorbed by the cell walls are identified. However, the signal coming from the chemical constituents of the wood is not detected with the experimental conditions used in this study.

An NMR experiment firstly consists in placing the sample studied in the homogeneous magnetic field B 0 . This is the basis of a macroscopic longitudinal magnetization M parallel to the intense field B 0 . Then, a radio frequency of short duration (a few microseconds), characterized by a rotating magnetic component B 1 , perpendicular to B 0 , causes the tilting of this magnetization on another plane. The longitudinal component of the induced magnetization will decrease, and the spins will be phased during the excitation of the atoms. This causes the transverse magnetization, which is the transverse component of the induced magnetization. At the end of the pulse, the nuclei being in an unfavorable energy level, magnetization will tend to return to its initial position. This is known as the relaxation phenomenon. An example showing the relaxation effects following a 90 ° pulse is shown in Figure 1.

There are two main types of relaxation, the transversal and the longitudinal relaxation.

Evolution of the magnetization M with relaxation and magnetic field B is given by the Bloch equations 1 and 2:

(1)

(2)

With T 1 and T 2 the time constants corresponding respectively to the longitudinal and transverse relaxation. [START_REF] Brownstein | Importance of classical diffusion in NMR studies of water in biological cells[END_REF] have proposed a theory that relates relaxation times to pore diameters. Therefore, each relaxation time measured in NMR corresponds to a given pore size. As a result, free water, located in the porosity, is easily discerned from the bound water that integrates into the cell walls.

X-ray Diffraction (XRD)

XRD is a technique that allows characterization of the crystalline structures of materials. A crystalline body is considered as an assembly of atomic planes more or less dense, called reticular planes. The characteristic distances separating these different planes are called interreticular distances. The existence of reticular planes and inter-reticular distances is at the basis of XRD.

When irradiating parallel beams of monochromatic X-rays, the atomic planes that form a crystalline material behave as a three-dimensional network. The X-ray beam is thus diffracted according to specific angles. This can provide information about the texture and degree of crystallinity of the material [START_REF] Liang | Method to evaluate the reactivity of microcrystalline quartz in carbonate rock[END_REF].

The theory of XRD was developed by W. L. Bragg. Its principle as well as its use are detailed in Bragg's relationship. Indeed, for an incident ray of monochromatic X-rays of wavelength λ, a ray is diffracted on Miller's index planes (h, k, l) if it satisfies Bragg's law [START_REF] Langford | Powder diffraction[END_REF]. The wavelength λ is between 0.1 Å for hard X-rays, and 50 Å for soft X-rays.

The equation governing the Bragg relation is given in equation 3. This equation links the Bragg angles θ, the wavelength λ, the inter-reticular distance of the family of planes {h, k, l} which diffract the incident rays d hkl and the diffraction order n.

(3)

The principle of Bragg's law is illustrated in Figure 2.

The characterization of crystalline materials' structures is made from the diffraction laws. In fact, two diffraction modes exist: when the source is monochromatic, the incident X-ray beam has X-rays of the same wavelength, whereas if the source is polychromatic, the incident X-ray beam has X-rays of different wavelengths.

The X-rays diffracted during the angular sweep of the irradiated surface of the sample are collected by an X-ray detector. The recorded signal provides the diffraction diagram, and represents the intensity measured by the detector as a function of the value of 2θ. From the angles that correspond to the observed diffraction peaks, the inter-reticular distances of the families of crystalline planes can be calculated, using Bragg's law. This is the way the crystallographic structure of the material is identified.

At the cell walls' scale of wood, there are various chemical constituents. Hemicellulose and lignin are present only in the amorphous phase, whereas cellulose exists in the wood in crystalline and amorphous form. As a result, the crystalline peaks that will be detected in XRD should correspond to the crystalline cellulose. The chemical formula of the crystalline cellulose is (C 6 H 10 O 5 ) n , and the corresponding chemical structure is shown in Figure 3.

The crystalline cellulose present in the cell walls, is of type I β and of monoclinic form, with the mesh characteristic parameters represented in table 1.

The XRD spectrum of crystalline cellulose, as well as the different crystalline planes are shown in Figure 4-a

Experimental protocol for bound and free water quantification using NMR

The originality of this study is to quantify the bound and free water within the wooden material, according to the FSP. To do so, a dry state (conditioning with silica gel), and relative humidity states of 78% and 95% were chosen for this study. This choice was made since the FSP is reached for a relative humidity greater than 78%, and smaller than 95% [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF]. In addition, these two levels are identical to those adopted in the X-ray micro-tomography experiment, elaborated in a previous paper [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF]. That way, the results' interpretations will be more complementary.

Six spruce wood samples of dimensions 2.5x2.5x20 mm 3 were prepared. Each sample was big enough to contain a complete reproduction cycle of wood, where early wood and latewood coexist. The 20 mm side of each sample was along the longitudinal direction of wood. All samples were first dried in dessicators using silica gel. Actually, the signal measured in NMR is proportional to the amount of material present in the tube. Therefore, each sample was placed in an NMR spectroscopy tube. In this way, as each sample has six sides exposed to the same environmental conditions, the mass balance of all samples will be achieved much faster than when only one equivalent large volume sample is used.

The NMR device used is a Spectroscopy-Bruker Minispec MQ20 ND-Series, with a 1 H proton resonance frequency of 20MHz, at CRM2 (Laboratoire de Cristallographie, Résonance Magnétique et Modélisations) in Nancy, France. It is shown in Figure 5-a. The diameter of the NMR tube is 10 mm. At each hygric state, a saline solution that guarantees a well-controlled relative humidity was placed in a small tube of 6 mm diameter, placed on a teflon support, all inside the NMR tube. This support has several holes in order to let moisture into the tube, and thus provide an environment with homogeneous humidity. The NMR tube was then closed and covered with a parrafin film to ensure that no moisture exchange will take place with the external environment. The assembly formed of the teflon support, the NMR tube, the tube containing the saline solution, and the sample is shown in As soon as the tube is put in the NMR spectroscopy device, the transverse relaxation time T 2 is measured. Then the longitudinal relaxation time T 1 is evaluated every hour. Finally, when the mass balance is reached, the transverse relaxation time T 2 is measured once again.

T 1 relaxation times were calculated using an inversion recovery procedure. 16 scans were performed per measurement point, with a 2 s recycle time, a gain of 90 dB, and 50 points logarithmically spaced by an echo time between 0.1 ms and 1 s. T 2 relaxation times were measured using a Carr Purcell Meiboom Gill (CPMG) procedure. 512 scans were performed, with a recycle time of 2 s, a gain of 90 dB, an echo time of 80 μs, and 500 consecutive echoes acquired.

At each relative humidity level, once the mass balance is reached, the small tube containing the saline solution is extracted from the NMR tube, and is replaced by another tube of the same size containing another saline solution. Whenever the salt solution is changed, the procedure described for calculating the longitudinal and transverse relaxation times between two specific relative humidity conditions is completed.

The post-processing of the longitudinal and transverse relaxation times is done afterwards by Laplace inversion, using UPENWIN software [START_REF] Borgia | Uniform-penalty inversion of multiexponential decay data: II. Data spacing, T2 data, systematic data errors, and diagnostics[END_REF]. The corresponding diagrams represent the amount of signal measured for each relaxation time. On the other hand, treatment with inverse Laplace transformation leads to fairly wide peaks. Each peak is representative of a pore size in which the hydrogen protons are located. As a result, the bound water located in the cell walls, and the free water located in the porosity can be identified, for each relative humidity level. On the other hand, for the spectroscopy, it is necessary to take into account the area under the peaks, and not the exact positioning of the peaks. This is shown in Figure 6-a. In addition, from the approximate positioning of the peaks, the nature of the water present in the material can be identified. For example, for the case of spruce wood, concerning the longitudinal relaxation time T 1 , the bound water appears around 1 ms, whereas the free water appears around 50 ms. These relaxation times are not the same for all materials, since they depend greatly on the size of the pores. Moreover, the signal integral under each peak (the area under the peak), leads to the quantification of the concerned water type (bound or free water).

In fact, the solid protons (that correspond to the solid matrix) are identified by very short T 2 transverse relaxation times, which are not detected with the experimental protocol used.

Nevertheless, these protons can have high T 1 relaxation times. As a result, it is necessary to verify that the T 1 longitudinal relaxation times of the peaks correspond to the T 2 crosssectional relaxation times of the bound water and the free water. Figure 6 Since the T 1 NMR measurements were made every hour until the sample reaches equilibrium, a comparative study can be made. For each measurement, the areas under the peaks that correspond to the bound and free water were calculated, in order to quantify both water types function of the time. The results are finally interpreted, making the link with the swelling of the wooden cell walls, which is a function of the bound water within the material [START_REF] El Hachem | Original experimental and numerical approach for prediction of the microscopic hygro-mechanical behavior of spruce wood[END_REF].

Experimental protocol for chemical constituents' evaluation using XRD

The interest of this experiment is to identify the crystalline phase in the material, and to understand the various interactions between the chemical constituents of spruce wood and humidity. Therefore, three cylindrical samples of 20 mm diameter and 2 mm thickness were prepared (the circular plane corresponding to the tangential plane of wood), and each one of them was big enough to contain a complete reproduction cycle of wood, where early wood and latewood coexist. The samples were conditioned in a climatic chamber with controlled temperature and relative humidity. The temperature was always 25 °C, and the relative humidity states tested along the sorption cycle were 10%, 60%, 72%, 78% and 87%. These levels were chosen in order to have refine relative humidity states around the FSP, as observed by [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF].

The used XRD device is Bruker D2 Phaser Bragg-Brentano (θ / 2θ) diffractometer at SATIE (Systèmes et Applications des Technologies de l'Information et de l'Energie) in ENS-Paris-Saclay, France. The device in shown in Figure 7-a. It is equiped with an X-ray tube, provided with a cobalt anticathode (λ = 1.78897Å) and a LYNXEYE linear detector with a maximum opening of 10°, allowing a significant gain in intensity and a fast acquisition of diffraction diagrams. Slots of 0.6 mm and 8 mm were respectively adopted for the source and the detector.

A sample holder with an X-ray transparent dome was used. The dome is tightened to the support through a waterproof rubber wall. An anti-diffusion knife was also used.

At each relative humidity, when a sample reached mass equilibrium, it was covered by the dome inside the climatic chamber, and then used to perform the corresponding XRD test. The used dome is shown in Figure 7-b. After the XRD test, the sample was placed in a climatic chamber programmed to a different relative humidity level.

The diagrams presented in this study were completed over an angular interval between 5° and 70°. An acquisition time of 1s per point and a pitch of 0.02 ° were chosen. A high voltage of 30 kV and a current of 10 mA were used.

For this study, the post-treatment was completed using Diffrac.EVA software. The purpose is to follow the evolution of the diffraction diagrams from one water state to another.

Therefore, crystalline and amorphous phases as well as instrumental background noise were identified. In the quantification method of the amorphous phase, the instrumental background noise was subtracted. Then the crystalline phase was recognized based on the PDF-2 sheets. Indeed, in the diffraction diagram containing several phases, the peak intensity ratio of the various constituents depends on the ratio of their proportion. This is the basic principle of quantitative XRD methods. Since the notion of intensity does not mean the use of a diffraction peak height, as this can lead to errors that can reach 20% [START_REF] Cyr | Détermination, par diffraction des rayons X, de la teneur en phase amorphe de certains matériaux minéraux[END_REF], the diffraction peak intensity means the area between the curve and the continuous bottom line. This net surface is known by integral intensity. Therefore, qualitatively, in the case of a phase mixture, the greater the area of the peaks of a phase, more this phase is present in large proportion. It is subsequently possible to identify the composition of the material in weight proportions. In this work, since cellulose is the only chemical constituent present in the crystalline phase, any detected crystalline peak will be assimilated to it. Concerning the amorphous phase, it can be quantified by assuming that the intensity of the halo due to the scattered X-ray disorder is proportional to the surface under the proper hump [START_REF] Cyr | Détermination, par diffraction des rayons X, de la teneur en phase amorphe de certains matériaux minéraux[END_REF].

The surface method consists of calculating the amorphous phase ratio as the ratio of the halo surface to that of the entire diffraction diagram. The different surfaces are calculated using Diffrac.EVA software, after correcting the diffraction diagrams.

The procedure consists therefore to first substract background noise from the diffraction diagrams. Then, the area of the amorphous phase is deduced by subtracting the area of the background noise diagram from the area under the curve of the halo. Finally, the crystalline phase is defined by subtracting the areas of background noise and amorphous phase from the total surface of the diffraction diagram.

Hence, the amorphous phase rate and the crystalline phase rate can be calculated. Each phase is then represented on a diffractogram. The crystalline phases are compared with those tabulated in PDF-2 sheets. The evolution of the crystalline and amorphous phases as a function of the relative humidity are studied. This is governed by the monitoring of the area under the curve which represents the amorphous phase, by the monitoring of the half-height width of the crystalline phase, and finally by the follow-up of the inter-reticular distances of the crystalline phase.

Results and discussion

Results of the detailed above experimental techniques and procedures will show first the follow up and quantification of water sorption, then the evolution of the cell walls chemical constituents along the sorption cycle, and will be followed by a global interpretation of those results.

Follow up and quantification of water sorption by NMR

The follow up of different water states (free, bound and constitutive) in the material in the sorption cycle was studied by NMR, and was achieved function of time, in order to study the evolution of the state of water. In fact, the dimensional variations of the cell walls caused by the environmental hygric variations are function of bound water presence within the material.

This part of the study will lead to a better understanding of the water state, and a better evaluation of the dimensional variations function of the relative humidity.

The NMR diagrams corresponding to the longitudinal relaxation time T 1 , obtained from a dry state up to a relative humidity of 78% are shown in Figure 8-a. Acquisitions were taken with a gap of one hour, except that in this figure an interval of 12 hours is retained for clarity purposes. For each acquisition, two peaks are discerned: the peak around 1 ms corresponds to the bound water, and the peak between 50 and 100 ms corresponds to the free water.

From a general point of view, the peak representing the bound water evolves with time. As already mentioned, it is the area under the peaks that must be considered. The areas under each peak, calculated at each acquisition, are shown in Figure 8-b. This figure shows that the amount of bound water, initially zero in the material, increases rapidly over time, until it reaches a constant value. On the other hand, the free water remains more or less constant. This is justified by the fact that the saturation point of the fibers corresponds to a relative humidity greater than 78% [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF].

In a next conditioning step, the samples being conditioned at 78% relative humidity, a small tube containing a salt solution that guarantees a relative humidity of 95% was placed in the NMR tube. The acquisitions were triggered immediately. Figure 9-a shows the NMR diagrams corresponding to the longitudinal relaxation time T 1 , obtained from a relative humidity of 78% to 95%. As before, acquisitions are presented with an interval of 12 hours.

The two peaks representing bound and free water are always present. However, in this case, both peaks evolve with time. The areas under each peak, calculated at each acquisition, are shown in Figure 9-b. This figure shows that the areas below the two peaks increase as a function of time. On the other hand, especially the first hours of the experience, the bound water kinetics is greater than that of free water. In fact, the saturation point of the fibers corresponds to a relative humidity between 78% and 95%, due to which, up to the relative humidity of 95%, the quantities of bound water and free water evolve at the same time within the material. In addition, since the saturation point of the fibers has not been reached at 78% relative humidity, the bound water is first absorbed faster than the free water.

Evolution of the chemical constituents on the sorption cycle using XRD

A study at the level of the chemical constituents of the cell walls was achieved. In fact, the different properties of these constituents are at the base of the physical phenomena that take place on the scale of the cell walls. Firstly, the complex structure of the cell walls consists of several layers and sub-layers, each having mainly specific proportions of cellulose, hemicellulose and lignin. The water absorption associated with this scale is therefore highly dependent on the physical parameters of the chemical constituents. As a result, the use of XRD technique helps to quantify and evaluate the evolution of the different phases proportions as a function of the relative humidity, which will also provide information about the reactions between chemical molecules and relative humidity. On the other hand, the distances between the crystalline molecules can be calculated according to the hydrous state of the material, giving more information about the local compressive or tensile stresses that may exist due to local dimensional variations of the cell walls [START_REF] Yamamoto | Origin of the characteristic hygro-mechanical properties of the gelatinous layer in tension wood from Kunugi oak (Quercus acutissima)[END_REF].

An example of a diffraction diagram obtained is shown in Figure 10. A first estimate of the background noise and the amorphous phase was first performed, in order to identify the crystalline phase. Then, the half-height width of the crystalline peak of the [002] plane and the corresponding inter-reticular distance were calculated. This approach was applied to each hygric state.

The curves representing the background noise and the amorphous phase for each treated water condition are shown in Figure 11-a. These results showed that up to a relative humidity of 78%, the amount of amorphous phase decreases as the relative humidity increases. On the other hand, between relative humidity of 78% and 87%, the amount of amorphous phase increases.

Several elements should be considered for the results' analysis. Actually, the bound water is absorbed by the cell walls, and particularly by the amorphous phase of the walls. This takes place through the breaking of hydrogen bridges in hemicellulose and amorphous cellulose.

As relative humidity increases, there is a possibility that hemicellulose and amorphous cellulose in the cell walls dissolve [START_REF] Kulasinski | Water adsorption in wood microfibril-hemicellulose system: role of the crystalline-amorphous interface[END_REF], causing a decrease in the amorphous phase (Figure 11-a).

As the FSP corresponds to a relative humidity between 78% and 95%, free water takes place within the material, as demonstrated in section 3.1, but appears in XRD diagrams as amorphous phase [START_REF] Brownstein | Importance of classical diffusion in NMR studies of water in biological cells[END_REF][START_REF] Borgia | Uniform-penalty inversion of multiexponential decay data: II. Data spacing, T2 data, systematic data errors, and diagnostics[END_REF][START_REF] Keeler | Understanding NMR spectroscopy[END_REF]. The increase of the amorphous phase between 78% and 87% relative humidity is therefore due to the presence of free water within the material.

The curves representing the crystalline phase below the peak [002] for each treated relative humidity condition are shown in Figure 11-b. The changes in inter-reticular distances function of the relative humidity are shown in Figure 12. In addition, the areas under the curves and the widths at half height of the peaks show no significant variation from a relative humidity to another. Figure 12 further shows that the inter-reticular distances decrease when the relative humidity increases. Indeed, the crystalline phase, representing only the crystalline cellulose, is hydrophobic [START_REF] Kulasinski | Water adsorption in wood microfibril-hemicellulose system: role of the crystalline-amorphous interface[END_REF][START_REF] Park | Cellulose crystallinity index: measurement techniques and their impact on interpreting cellulase performance[END_REF]. It is well oriented and organized to prevent interactions with moisture. Its amount remains therefore more or less constant. This is proved by verifying the width at half height of the crystalline peak [002] and by the area under the curve. On the other hand, the reduction of inter-reticular distances from one hygric state to another can be caused by mechanical stresses of compression in the secondary layer of the cell walls due to swelling. In fact, swelling dimensional changes of hemicellulose and amorphous cellulose generate compressive stresses at the level of the crystalline cellulose in the secondary layer. This can be the cause of the inter-reticular distances decrease in the cellulose mesh.

Global interpretation

It has been shown, by the use of the NMR technique, that below a relative humidity of 78%, there was only bound water within the material and the FSP was not yet reached. On the other hand, between relative humidity of 78% and 95%, the bound and the free water coexisted within the material. The FSP therefore corresponds to a water content that has been reached between these two levels of humidity. Hence, it is consistent that above 78% relative humidity, the X-ray micro-tomography results obtained in the previously published paper [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF] were influenced by the presence of free water in the porosity.

On the other hand, the XRD reported constructive explanations regarding the shrinkage of cell walls observed below the FSP. Indeed, when absorbing moisture, a certain proportion of the hemicelluloses may dissolve due to hydrogen bridges rupture, or may swell and tend to generate compressive stresses at the level of the secondary wall, and more particularly on the crystalline cellulose. These stresses cause a decrease in the inter-reticular distances of the cellulose molecules. These phenomena can generate the local swelling along the longitudinal axis of the wood cell walls and the local shrinkage in the transversal plane of the cell walls.

It is after all important to evaluate the coupling of the previous results, concerning the bound and free water uptake within the material, the chemical components evolution due to the increase of the relative humidity, and the previous results obtained by [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF] about the microscopic morphological changes along the sorption cycle, many interesting. The results elaborated in this paper are totally consistent with previously described results based on X-ray micro-tomography experimental study [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF]. Therefore, it is to be noted that below the FSP, when the relative humidity increases, spruce wood is subjected to a macroscopic swelling, whereas, at the microscopic scale, the cell walls absorb bound water through the hemicellulose, which leads to local swelling of the cell walls along the longitudinal axis and to a stretching of the cell walls in the transversal plane. This leads to a shrinkage of the walls in that plane (the thickness of the cell walls decreases), and to an increase of the pore diameter. On the other hand, after reaching the FSP, there is no more obvious macroscopic swelling, whereas, at the microscopic scale, the cell walls are already saturated in bound water, and the free water takes place in the porosity holes.

Conclusion

Two experimental campaigns have been conducted in this study for a better characterization of the hygro-mechanical interactions within the wood anatomical structure. At a first stage, NMR was investigated to evaluate and quantify the increase in bound and free water within the material. This study was done along the sorption cycle and for different relative humidity levels. XRD was then performed for a better understanding of the interactions of moisture with the chemical constituents of wood cell walls.

The coupling of these two experimental campaigns, in addition to the works elaborated by [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF] has led to interesting results. The results showed an increase of bound water within the material between dry state and the FSP. Once the FSP is reached, the cell walls are already saturated in bound water, and the increase in relative humidity induces an increase in free water within the material. On the other hand, below the FSP, hemicellulose absorbs moisture and swells. In addition, a certain proportion of hemicelluloses may dissolve following the rupture of the hydrogen bridges. As for the amount of crystalline cellulose, it did not vary according to the relative humidity. A small narrowing of the crystalline cellulose molecules was observed, due to compressive stresses exerted by the swelling of the hemicellulose. These results are in agreement and even justify the results obtained with [START_REF] El Hachem | Microscopic swelling analysis of spruce wood in sorption cycle[END_REF].

These results also provide insight into the interaction of relative humidity and wood with more precision on the nature of water in the wood cells and its dimensional variations during swelling. These interactions should be considered seriously when modeling the multiphysical transfers phenomena inside this material.

  , based on the PDF-2 sheets (Powder Diffraction File 2), provided by the ICDD (International Centre for Diffraction Data). The diffraction graph shows two obvious peaks for 2θ values of 17.2 ° and 26.4 °. For the interpretations in sections 3.2 and 3.3, we will refer to the 2θ main peak of 26.4°, corresponding to plane [002]. Furthermore, some studies have been completed in the literature concerning the crystalline form of cellulose. For example, some values of inter-reticular distances and cellulose crystallinity, were calculated by Ju et al. (Ju et al., 2015): -d 002 = 0.389 -d 10ĩ = 0.531 -d 101 = 0.593 The geometric shape of the cellulose, as well as the main plane [002], and the inter-reticular distances d002 are shown in Figure 4-b-c.

Figure

  Figure 5-b.

  -b shows the NMR diagram of T 2 transverse relaxation times, in which two peaks are distinguished. The existence of these two peaks in this diagram validates the presence of bound and free water identified in Figure 6-a.


Chrystel Dezayes

Catherine Lerouge

Philippe Lach

Characterization of the basement / cover unconformity based on Saint Pierre Bois quarry observations (Rhine Graben, France)

est destinée au dépôt

Introduction

According to the literature, deep geothermal energy reserves are associated with naturally fractured/faulted reservoirs in various geological terrains. Fractures and faults provide pathways for fluid flow, which are crucial when primary porosity is low, such as in granite. Faults may be more permeable than fractures, but the fracture network may provide a more porous system than the former. Even though fractured reservoirs have been investigated for more than fifty years, it is difficult to model reservoir behavior because we lack knowledge of the fracture network efficient for fluid flow generally due to complex structural/tectonic history and episodic water/rock interactions associated with this history.

In this abstract, we present a combined works of structural and mineralogical analyses of the transition zone between the crystalline basement and its sedimentary cover. This zone is around 3km within the Rhine Graben and constitute a target for deep geothermal energy with a reached temperature between 150°C and 180°C, such as at Landau and Insheim in Germany, and at Soultz and Rittershoffen in France.

This characterization of the basement/sedimentary transition zone was performed at the Saint Pierre Bois quarry (Vosges, France), located on the western shoulder of the Rhine graben, in Central Vosges. The aim of this work is to define the major parameters that control the efficiency of a fracture network for fluid flow, based on:

-Description of the fracture network including fracture orientation measurements, fracture density, fracture distribution, microstructures, and deformation style, fracture filling, -Chronology of cementation, -Study of porosity and water/rock interactions at the fracture/fault scale.

The Saint Pierre Bois quarry is of particular interest for its contrasts in weathering, deformation, and rheology, and for three benches created by aggregate exploitation that provide a 3D picture of the fracture network in both granite and the Stephanian sedimentary cover, and of connectivity between them. The Hercynian granite belongs to the generation of two mica leucogranite emplaced along major crustal accidents during the late Carboniferous in the Central Vosges, 322-330 Ma ago. By measuring fractures and characterizing wall rock alteration and mineral filling, we are able to propose a conceptual model of the quarry-scale fracture network that integrates fluid circulation with geological time.

Methods

To characterize a fracture network, several key parameters are needed: fracture orientation (strike and dip), fracture distribution, aperture (open, filling), mineral filling, length (extent, height) and connectivity (abutment, cross-cut, termination). As is commonly done, we used the scanline method to describe the fracture network in the quarry and performed a statistical analysis of special distribution. A total of 229 fracture orientations were measured on all the walls of the quarry, mainly in the granite, which is more accessible. Infills of fractures were sampled in relation to their orientation and a petrological characterization was performed by optical microscopy and complementary scanning electron microscope.

Porosity data was acquired for granite and sandstone matrix in the Saint Pierre Bois by [START_REF] Kushnir | Characterizing the physical properties of rocks from the Paleozoic to Permo-Triassic transition in the Upper Rhine Graben[END_REF] and was completed by measurement of the infill fractures. Four techniques were applied to different sample types: triple-weighing method, mercury porosimetry, helium pycnometry and digital image processing applied to back-scattered electron.

Results

For all scanline directions measured in granite, the fracture patterns are equivalent in terms of orientation, whereas the density of each set varies. This suggests that our sampling is appropriate and avoids a high sampling bias. Two main/primary sets have been identified: F1, N90°E±20° nearly vertical and F2, N30°E±20° dipping 55° to SE. Fractures that do not belong to these sets are grouped in the F3 set, N150°E±40°.

Fracture frequency is approximately 0.7±0.07 fr/m and the most common spacing more than 1 m, except for UW250 scanline, which shows a higher fracture frequency, approximately 1 fr/m and a primary spacing less than 1 m. Whereas the coefficient of variation indicates a random fracture distribution, fracture clusters are identified (Figure 1).

Granite samples collected along major E-W tectonic corridors consist of cohesive cataclasite and breccia. Microquartz is the major cement in cohesive cataclastic zones in association with micron sized iron hydroxides and minor clay minerals. Clay minerals consist of illite, essentially, whose origin remains unknown. The same cements are present in the breccia; however, additional medium size subeuhedral to euhedral quartz grains also fill residual porosity. Rare macroscopic barite crystals (0.1 to 1.5 cm) are associated with euhedral quartz-filled fractures. At the scale of the quarry, reddish or white/greenish coatings are present on numerous fault/fracture surfaces where striations are occasional present. These coatings are a mm-to cm-thick cohesive cataclasite of granite cemented by silica/microquartz. White/greenish/reddish colors are mainly due to the absence/presence of iron hydroxides.

Throughout the quarry, granite wall rock alteration due to fluid circulation through fractures/faults rapidly decreases within several centimeters. In the highly weathered-fractured zone, the high fracture density induces a large volume of altered granite matrix. Close to the fractures and faults and in the fault corridor, biotite and in lesser amounts, muscovite are altered to clay minerals (illite, illite-smectite mixed layers) with titanium oxides and iron hydroxides, K-feldspar is partially dissolved, and primary plagioclase is replaced by illite with µm-sized aluminum-phosphate-sulfate (APS) minerals.

Porosities of breccias and cataclasites sampled in the structures range from 1.65 to 7.45 %. The lowest porosity values were measured in cataclasites. Porosity values of the bulk granite matrix in the core of the breccia corridors are rather low, whereas samples collected on the margins of a breccia corridor have bulk water porosities of 4.50 and 7.45 %. In the cohesive cataclasite, reddish zone cemented by iron hydroxides has He porosity of 6.74 %, whereas highly cohesive orange zone cemented by silica/microquartz has He porosity of 2.13 %.

The well in-place sedimentary cover overlying granite is accessible only in the northern part of the quarry and even there it is not easily accessible. Therefore, only 9 fracture orientations were measured, which does not constitute an adequate sample for statistics but does represent the overall orientation. The primary set is N10°E±10° and the secondary set is N100°E±10°, all are nearly vertical.

Sedimentary cover consists of heterogranular medium to coarse-grained arkoses. Detrital minerals consist mainly of quartz, plagioclase, and K-feldspar grains, with minor micas, clay matrix, iron hydroxides, and heavy minerals. Feldspars, and more specifically plagioclase, exhibit a variable degree of dissolution. Authigenic minerals are primarily quartz with minor K-feldspar, clay minerals including illite, aluminum-phosphate-sulfate (APS) minerals, galena and pyrite. APS and illite essentially fill porosity of plagioclase dissolution. A generation of micron-to medium-grained quartz crystallized at sandstone grain contacts indicates a protocataclastic deformation in the weathered-fractured zone in the eastern part of the quarry.

1 st Geoscience & Engineering in Energy Transition Conference 16 -18 November 2020, Strasbourg, France Bulk porosity and dissolution porosity of both types of feldspars measured in arkoses range between ~5 and 14 %. Total porosity is relatively well correlated with grain size. In all cases, plagioclase dissolution porosity is significantly higher than K feldspar dissolution porosity.

Discussion

The study area is located in a highly fractured region affected by Hercynian tectonics and the Tertiary opening of the Rhine Graben. Whereas the main directions are equivalent in the Saint Pierre Bois quarry than in the all scales in the Rhine graben and its surrounding, the distributions of each set are different in the Vosges Massif and in the Rhine Graben. The Hercynian directions, mainly NE-SW and NW-SE, are more common than the Rhenish direction in the Vosges Massif. However, the E-W direction, which is mainly present in the Saint Pierre Bois quarry and its surroundings, is less frequent at the Rhine graben scale. This fracture set could be related to Hercynian tectonics and linked to the direction of the Villé basin, which opened during the Permian at the end of the Hercynian orogen.

Inside the quarry, two areas have been differentiated: the fractured-weathered zone to the east and a more structured zone in the western and northern parts (Figure 1). These zones are part of a larger fault zone, and are defined as the core zone and damage zone, respectively [START_REF] Bauer | Architecture, fracture system, mechanical properties and permeability structure of a fault zone in Lower Triassic sandstone, Upper Rhine Graben[END_REF]. The damage zone includes three dominant fracture directions, sometimes grouped in clusters. This fracture sets can be linked to the Riedel model. Therefore, we consider the F1 fracture set (N90°E±20°) to correspond to shear fracture in the R-direction. Likewise, the F2 fracture set (N30°E±20°) corresponds to P' fractures and the F3 fracture set (N150°E±40°) to P fractures. This implies that the large fault zone should have a N80°E direction with dextral motion during its formation in an ESE-WNW σ 1 stress field (Figure 1). The three alteration stages are present in the core zone and damage zone and in granite and arkose, indicating that fluid circulation has been active in all zones since the deposition of arkose, due to pre-Triassic granite disaggregation and weathering. However, the intensity of each alteration/weathering stage varies in the damage zone and core zone, and with rock type (granite/arkose) in relation to fluid pathways (Figure 1).

Conclusion

Based on an integrated geological, structural and petrological study at the quarry scale, we propose a conceptual model of the fracture network at the sedimentary cover/basement transition that integrates fluid circulation with geological time (Figure 1).

This model could help the geothermal exploration in the deeper part of the Rhine Graben by identified the main circulation/plugging zones at the transition zone of the crystalline basement and the sedimentary cover.

  1 st Geoscience & Engineering in Energy Transition Conference 16 -18 November 2020, Strasbourg, France

Figure 1

 1 Figure 1 Main fracture direction in the Saint Pierre Bois quarry. The color rectangles on the wall represent the fracture clusters identified by the scanline analysis. Red: F1 set; green: F2 set; blue: F3 set. The large hatched sector represents the altered-fractured part of the quarry.

  & Engineering in Energy Transition Conference 16 -18 November 2020, Strasbourg, France Petrology of fracture filling minerals in the granite and diagenetic minerals cementing the overlying arkoses provide evidence of common cementation stages, strongly suggesting fluid circulation between the two lithological compartments. A first common stage is microquartz and illite, followed by coarser quartz (euhedral in granite, xenomorphic in arkose) ± barite. Illite formation indicates temperatures of at least ~130-150°C[START_REF] Liewig | K-Ar dating of varied microtextural illite in Permian gas reservoirs, northern Germany[END_REF]. Therefore, this cementation stage results from deep circulation of hot fluids associated with high fluid/rock ratios at the cover/granite interface. A second common stage consists of hydroxides with minor clays and microquartz that fill residual porosity in fractures, faults, and breccia corridors in granite and in the arkose matrix. The presence of iron hydroxides trapped in cements filling porosity of plagioclase dissolution in granite testifies to pre-Triassic granite weathering conducting to the arkose deposition. All these observations highlight three major fluid circulation stages and fluid/rock interactions: (1) pre-Triassic weathering, (2) deep hot fluid circulation at depth (hydrothermal stage), and (3) recent weathering.

Acknowledgements

This study was co-funded by ANR (French National Research Agency) in the framework of the ANR project CANTARE-Alsace (ANR-15-CE06-0014). The authors would like to thank Jérémie De Bonneval, Director of Leonhart for permitting the access of the quarry.


Keywords: System identification, Nonlinear systems, Volterra series, Least-squares identification

This article addresses identification of nonlinear systems represented by Volterra series. To improve the robustness of some existing methods, we propose a pre-processing stage that separates nonlinear homogeneous order contributions from which Volterra kernels can be identified independently. Unlike existing separation methods that use amplitude relations between test signals, we propose another order separation method based on phase. This method gives access to a new type of complex-valued output signals, which can be used to improve kernel identification. First, the underlying ideas are introduced via the presentation of a theoretical method using complex-valued test signals. Second the proposed order separation method using real-valued signals is described. Third, a new identification process is given, combining existing leastsquares identification method with the previous results. Finally, numerical experiments are used to compare the proposed order separation method with state-of-the art, as well as to evaluate the new Volterra series identification process.

Introduction

Nonlinear system identification is a very challenging topic. One approach, using Volterra series [START_REF] Volterra | Sopra le funzioni che dipendono da altre funzioni[END_REF], has begun in the 40's with the seminal work of [START_REF] Wiener | Response of a nonlinear device to noise[END_REF]. One characteristic that led to the widespread use of Volterra series is the large class of models that they can represent [START_REF] Boyd | Fading memory and the problem of approximating nonlinear operators with Volterra series[END_REF]. But this advantage comes with important difficulties concerning identification, such as:

• the high number of kernels coefficients to estimate, which increases drastically with respect to the truncation order and the kernels memory; this leads to high variance on the estimated parameters. • the simultaneous estimation of all kernels from the direct output (as shown in Fig. 1a); indeed, for a system that allows a Volterra series representation, amplitudes generally decreases as the order n increases, thus making estimation of the higher-order kernels less reliable.

The first problem can be alleviated in different ways. Tensor decomposition [START_REF] Batselier | Tensor Network alternating linear scheme for MIMO Volterra system identification[END_REF][START_REF] Favier | Nonlinear system modeling and identification using Volterra-PARAFAC models[END_REF], polynomial kernel regression [START_REF] Franz | A unifying view of Wiener and Volterra theory and polynomial kernel regression[END_REF] or projections on orthogonal basis, such as Laguerre [START_REF] Marmarelis | Identification of nonlinear biological systems using Laguerre expansions of kernels[END_REF], Kautz (da Rosa, Campello, & Amaral, 2007) or Generalized Orthogonal Basis [START_REF] Kibangou | Selection of generalized orthonormal bases for second-order Volterra filters[END_REF], greatly reduces the number of parameters to estimate, and therefore the overall complexity. The system can also be approximated using a block structure [START_REF] Schoukens | Identification of block-oriented nonlinear systems starting from linear approximations: A survey[END_REF] , e.g. Hammerstein, Wiener or Wiener-Hammerstein structure, where the system is decomposed with only linear filters and static nonlinearities; for specific structures, estimation methods with good results are available [START_REF] Rébillat | Identification of cascade of Hammerstein models for the description of nonlinearities in vibrating devices[END_REF][START_REF] Tiels | Wiener system identification with generalized orthonormal basis functions[END_REF]. Alternatively, regularization methods can be used to improve kernels estimation [START_REF] Birpoutsoukis | Efficient multidimensional regularization for Volterra series estimation[END_REF][START_REF] Birpoutsoukis | Regularized nonparametric Volterra kernel estimation[END_REF].

The second problem is more difficult to circumvent. One approach is to divide the identification into 2 steps (as shown in Fig. 1b):

1. Order separation: nonlinear homogeneous orders y n are separated from a set of output signals; 2. Kernel identification: Volterra kernels h n are identified separately from each signal y n .

This two-step approach decreases computational cost and increases robustness of the overall kernels estimation by splitting the identification process into N smaller and independent System u Identification problems. Furthermore, this procedure is modular, meaning that any existing identification methods can be used in step 2. But it has to be noted that, due to the fact that order separation methods require specific user-defined input, this approach can not be used with any set of measurements (e.g. pre-recorded datasets).

y {h 1 , h 2 , • • • , h N } (
In the literature, order separation is processed by using a collection of amplitude gains applied to a common input signal [START_REF] Halme | Polynomial operators in non-linear systems theory[END_REF][START_REF] Simpson | Correlation techniques for the identification of non-linear systems[END_REF]; it has been used for identification purposes in e.g. [START_REF] Bard | Modeling of Nonlinearities in Electrodynamic Loudspeakers[END_REF]; [START_REF] Boyd | Measuring Volterra kernels[END_REF]; [START_REF] Lee | Estimation of non-linear system parameters using higher-order frequency response functions[END_REF]; [START_REF] Zhang | Volterra series truncation and kernel estimation of nonlinear systems in the frequency domain[END_REF]. But this method relies on the inversion of a real-valued Vandermonde matrix, which is ill-conditioned when its size grows (Pan, 2016). In practice, this method requires numerous amplitude gains which have to be carefully chosen in order to provide good separation.

Recently, a separation method using phase deviations and amplitude gains has been proposed [START_REF] Bouvier | Nonlinear Homogeneous Order Separation for Volterra Series Identification[END_REF], for which previous remarks on conditioning of the Vandermonde matrices still apply. Those ideas have been further developed in [START_REF] Bouvier | Homophase signals separation for Volterra series identification[END_REF] leading to another kind of separation into homophase signals, which has been used to develop new identification processes. The use of phase deviations can be related to the results of [START_REF] Gardner | Exploitation of cyclostationarity for identifying the Volterra kernels of nonlinear systems[END_REF], where, in a probabilistic framework, cyclostationary signals are used to obtain orthogonality between orders1 (in a cross-correlation sense). But even with those new developments, robust order separation can not yet been achieved: the method of [START_REF] Bouvier | Nonlinear Homogeneous Order Separation for Volterra Series Identification[END_REF] still uses a set of amplitude gains leading to conditioning problems, and the method proposed in [START_REF] Bouvier | Homophase signals separation for Volterra series identification[END_REF] does not separate homogenous orders. This paper presents a new parameter-free order separation method exclusively based on phase deviation between excitation signals. This separation depends on the inversion of a 2D Discrete Fourier Transform and sparse matrices. Its contribution to Volterra series identification is illustrated with the introduction of a new Least-Squares identification method. All presented algorithms are available in an open-source Python toolbox 2 . This paper is organized as follows: in Section 2, we give some mathematical foundations for the Volterra series paradigm and recall the existing separation method. Section 3 recalls basics on the phase-based order separation for the abstract case of complex-valued inputoutput signals; then, its robustness properties to noise and truncation error are analysed. Section 4 presents the main result, that is, an adaptation of the separation method to real-valued input-output signals to be used on physical systems. From this result, Section 5 derives a new Volterra series identification method. Finally, evaluation and comparison between several methods are made in Section 6, before giving some conclusions and perspectives.

Volterra series and order separation

Notations and functional setting

Introduce the following notations and definitions (with n ∈ N * where N * = N \ {0}):

• L p (X, Y): standard Lebesgue-space of p-normed functions from vector spaces X to Y;

• T: set of times R, R + or [0, T ] with T > 0; • U = Y := L ∞ (T, R): sets of input and output signals; • VS: set of series {h n } n∈N * with h n ∈ L 1 (T n , R) such that • ρ h > 0 is a non zero convergence radius of • φ h : x → n∈N *
h n 1 x n , the gain bound function;

• V n := ML n U, Y : set of bounded (continuous) n-linear operators V n : U n → Y.

Volterra series

This section presents some recalls on Volterra series for bounded-input bounded-output (BIBO) systems. More details can be found in e.g. [START_REF] Boyd | Analytical foundations of Volterra series[END_REF]; [START_REF] Rugh | Nonlinear system theory[END_REF].

Definition 2.1 (Volterra series). A BIBO causal time-invariant system is described around its equilibrium y 0 ∈ R by a Volterra series {h n } n∈N * ∈ VS if, for all inputs u ∈ U such that u ∞ < ρ h , its output y is given by, for all t ∈ T,

y(t) = V [u] (t) := V 0 =y0 + ∞ n=1 V n [u, . . . , u] (t) =yn(t) , (1)
with V 0 the constant operator of value y 0 , and the Volterra operators of order n ∈ N * defined in V n by

V n [u 1 , . . . , u n ] (t) := T n h n (τ 1 , . . . , τ n ) n i=1 u i (t -τ i ) dτ i , (2)
where h n are called the Volterra kernels.

Note that equilibrium y 0 is usually chosen as zero, using the change of output y → y -y 0 . In the sequel, y 0 is kept non zero for sake of completeness, and all the y n 's are called nonlinear homogeneous order contributions (or orders for short). It has been shown in [START_REF] Boyd | Fading memory and the problem of approximating nonlinear operators with Volterra series[END_REF] that a Volterra series operator V can be used to approximate any fading-memory time-invariant continuous nonlinear system. In practice, bounds of convergence radius ρ h can be computed for large classes of well-posed systems [START_REF] Brockett | Convergence of Volterra series on infinite intervals and bilinear approximations[END_REF][START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linear-analytic systems and application[END_REF]. Volterra kernels of a system are not unique: permuting time variables in h n preserves y n in (1). For identification issues, it is convenient to use uniquely-defined kernels [START_REF] Rugh | Nonlinear system theory[END_REF]. Symmetric ones are those invariant with respect to permutations, i.e. such that for all permutations π of n variables,

h n τ π(1) , . . . , τ π(n) = h n τ 1 , . . . , τ n .
(3)

To improve identification of V from measurements, this paper exploits the following straightforward properties.

Property 1 (Operator V n ). Let V n ∈ V n as defined in (2). Then the following properties holds:

(i) Symmetry: if h n is symmetric (w.r.t. its time variables) then V n is symmetric (w.r.t. its input signals);

(ii) n-linearity: for all signals u 1 , . . . , u n , v and scalar α,

V n [u 1 , . . . , u k + α v, u k+1 , . . . , u n ] = V n [u 1 , . . . , u n ]+α V n [u 1 , . . . , v, u k+1 , . . . , u n ] , (4) with 1 ≤ k ≤ n;
(iii) Homogeneity: for all signals u 1 , . . . , u n and scalar α,

V n [αu 1 , . . . , αu n ] = α n V n u 1 , . . . , u n .
(5)

In the sequel, kernels h n and operators V n are considered symmetric and convergence conditions are assumed to be met.

Amplitude-based Separation (AS) Method

This section recalls the order separation method introduced in the 70's in [START_REF] Halme | Polynomial operators in non-linear systems theory[END_REF]; [START_REF] Simpson | Correlation techniques for the identification of non-linear systems[END_REF]. Denote V ∈ VS the operator of a truncated Volterra series with degree N . Consider a reference input signal u ∈ U and a set of gains α k ∈ R \ {0}, k = 0, . . . , N . Introduce the collection of input signals 1) and ( 5), it follows that

u k = α k u ∈ U output signals z k = V [u k ] ∈ Y. From (
     z 0 z 1 . . . z N      (t) =      1 α 0 . . . α N 0 1 α 1 . . . α N 1 . . . . . . . . . . . . 1 α N . . . α N N      •      y 0 y 1 . . . y N      (t) z(t) = A • y(t) , (6)
where the y n 's are the homogeneous orders of the truncated series in response to the input signal u.

Since A is a Vandermonde matrix, it is invertible if and only if the α k are all different. In this case, it is possible to recover orders y n . In the following, this method will be referred as the Amplitude-based Separation (AS) method.

In practice, inverting the Vandermonde matrix A poses numerical difficulties due to its bad conditioning (Pan, 2016). Those can be partially alleviated by using Newton or Lagrange Recursive methods (Golub & Van Loan, 2012, Algorithm 4.6.1 and 4.6.2), or considering more than N + 1 gains (leading to more test signals) and using a pseudo-inverse matrix. Also, making use of several amplitudes and a unique waveform poses several limitations: for high amplitudes, the system can overload or be excited in a regime outside the Volterra series convergence domain; for small ones, higher orders are more and more hidden in measurement noise. Thus a trade-off has to be made between the condition number quality of A and the amplitude range narrowing. To the knowledge of the authors, the truncation order N exceeds

Complex-valued Phase-based Separation (CPS)

This section first presents ( § 3.1) some recalls on the phase exploitation for order separation in the case of 'abstract' (non physical) case of complex-valued inputs and outputs:

U C = Y C := L ∞ (T, C).
Second, in § 3.2, the main properties of this method are analysed and stated.

Method presentation

Denote V ∈ VS the operator of a truncated Volterra series with degree N . Consider a reference input signal u ∈ U C and w = e j2π/(N +1) the first (N + 1) th root of unity. Replacing gains α k by w k , k = 0, . . . , N in the AS method leads to

     z 0 z 1 . . . z N      (t) =      1 1 . . . 1 1 w . . . w N . . . . . . . . . . . . 1 w N . . . w N 2      •      y 0 y 1 . . . y N      (t) z(t) = W • y(t) , (7)
where W is the Discrete Fourier Transform (DFT) matrix of order N + 1 with normalization factor 1. Its optimal conditioning makes the recovery of orders y(t) = W -1 z(t) robust, contrary to the AS method. Note that in ( 7), the DFT is applied to homogeneous orders y n (produced by the reference input signal u = u 0 ) with respect to order variable n (and not the time variable). This method is referred below as the Complex-valued Phase-based Separation (CPS) method.

From the formal algebraic point view, AS and CPS methods look very similar. But from the signal point of view, they operate in very different ways: AS exploits amplitudes (whatever the signal shape) as a contrasting parameter of homogeneous orders on a validity range to be maximal, whereas CPS exploits shape (whatever the amplitude) through phase deconstruction and reconstruction. Differences on excitations can be summarized by comparing factors α k , w k and their power in the complex plane (see Fig. 2): in CPS, all the factors and their power lie on the unit circle and are one of the N + 1 unit-roots; in AS, they all lie on the real axis and powers of α k converge towards 0 (if

|α k | < 1) or diverge (if |α k | > 1)
as the truncation order N increases.

Properties of CPS method

Due to the properties of the DFT, the CPS method comes with interesting properties.

Zero-padding and number of phases

As will be shown in § 3.2.2 and § 3.2.3, the number of phases P has a positive impact on the robustness of the CPS method. For a given maximal order N , it can be artificially increased using zero-padding. Let the number of phases P be greater than P min = N + 1. Introduce w = e j2π/P , and the collection of input signals

u k = w k u ∈ U C , k = 0, . . . , P -1 and output signals z k = Re Im α 1 α 2 α 3 α 4 w 0 = α 0 w 1 w 2 w 3 w 4 Figure 2.
Powers of typical factors used in AS and CPS methods, compared in the complex plane: (AS) red circles ( ) for α > 1; (CPS) blue squares ( ) for w = e j2π/(N +1) with N = 4. 7) remains valid, where W is now the DFT matrix of order P and the vector of orders is y = y 0 . . . y N 0 . . . 0 T , with P -N -1 zeros. Considering more than P min phases increases the number of tests signals K (equal to P for CPS method).

V [u k ] ∈ Y C . Then relation (

Nonlinear order aliasing

Now, consider that V ∈ VS is the operator of an infinite Volterra series. Consider a number of phases P ≥ P min = N + 1, where N is the chosen approximation order. Since w n+rP = w n for any P th root of unity w and for all integers n and r, property (5) yields

V [wu] (t) = P -1 n=0 w n ∞ r=0 y n+rP (t) , (8)
and (7) becomes

z(t) = W y(t) + ∞ r=1 y (r) (t) , (9)
with y = y 0 . . . y P -1 T and y (r) = y rP . . . y rP +P -1 T . This reveals that order y n is perturbed by a residual term ∞ r=1 y rP +n , which is structured as a nonlinear order aliasing. This effect permits to use the amplitude as a contrast factor to control the Signal-to-Residual Ratio (SRR). Replacing gains w k by numbers γ k = ρw k , with 0 < ρ < 1, leads to

z(t) = W • R • y(t) + ∞ r=1 ρ rP y (r) (t) , (10)
with R := diag ρ 0 , ρ 1 , . . . , ρ P -1 . Therefore, in estimated orders y(t) = R -1 W -1 z(t), there is a ratio 1/ρ P between the desired signal and the first residual y (1) (t). This ratio permits to improve the SRR by choosing a high P or a low contrast factor ρ. For example, a SRR gain of 20 dB is obtained by setting P = 100 and ρ ≈ 0.977.

Noise reduction

Suppose that each measured output z k is perturbed by an independent and identically distributed (i.i.d.) complex Gaussian noise, i.e.

z(t) = z(t) + ε r (t) + jε i (t) , (11)
where ε r and ε i are both multivariate i.i.d. Gaussian noise of null mean and covariance matrices Σ r/i = σ 2 r/i I, with I the identity matrix. The best estimator y of y in a least meansquares sense is

y(t) = R -1 W -1 z(t) = y(t) + ε r (t) + j ε i (t) , (12)
where ε r/i are multivariate i.i.d. Gaussian noise of null mean and covariance matrices

Σ r/i = σ 2 r/i P R -2 . ( 13
)
Equation ( 13) shows that using CPS method, Signal-to-Noise Ratio (SNR) between estimated order y n and error due to noise measurement is reduced by a factor ρ n √ P . For example, setting P and ρ as in the previous example ( § 3.2.2) and considering a maximum order of N = 9 leads to a SNR gain going from 20 dB (for the zeroth order) to 18.2 dB (for order 9). This illustrates how increasing the number of phases P (and therefore the number of test signals K) only leads to a more robust order separation method.

But the CPS method is abstract and useless in practice since only real-valued signals can be used on a physical system. Its recast from complex-valued to real-valued signals is not possible because equation ( 7) cannot be recomposed using an independent processing of real and imaginary parts.

Real-valued Phase-based Separation (RPS)

This section provides a method based on the phase exploitation, adapted to the use of realvalued signals for handling real-life cases:

U = Y := L ∞ (T, R).
The main theoretical result of this paper is stated in Theorem 4.4 ( § 4.1). It allows the derivation of the main practical result (RPS method, § 4.2), the properties of which are detailed in § 4.3.

Preliminary results

In this section, let u ∈ U C be a complex-valued signal, and operators V n ∈ V n be associated with symmetric real-valued Volterra kernels h n . Definition 4.1 (Interconjugate signal). The interconjugate signals y n,q ∈ Y C , q = 0, . . . , n, are defined by y n,q (t) := V n u, . . . , u n-q times , u, . . . , u q times (t) ,

where a denotes the complex conjugate of a. By extension, we also denote y 0,0 (t) := y 0 .

Remark 1. Signals y n,0 (respectively y n,n ) correspond to the output of order n of a series excited by signal u (resp. u). Signals y n,q with 0 = q = n do not have similar correspondence in terms of input/output relationships.

Property 2. By symmetry of V n , it is straightforward to show that y n,n-q (t) = y n,q (t), and that for even n, signal y n,n/2 is real.

Interconjugate signals allows the rewriting of the output of order n for a real-valued input signal, as the following lemma shows.

Lemma 4.2. Let r = Re [u] ∈ U. Then y n = V n [r, . . . , r] ∈ Y is given by y n (t) = 1 2 n n q=0 n q y n,q (t) , (15)
where n q = n! q!(n-q)! is the binomial coefficient.

Proof. The proof stems from the symmetry and n-linearity of operators V n fed by r(t) = (u(t) + u(t)) /2.

Example 4.3. Consider the cosinusoidal signal r(t) = cos(2πf t) and its analytic version u(t) = e j2πf t . Interconjugate signals y n,q are complex sinusoids of frequency (n -2q)f , with phase and amplitude depending on operator V n . Furthermore, if V n is such that V n [r, . . . , r] (t) = r(t) n , (15) leads back to the Chebyshev expansion of cos(2πf t) n .

Lemma 4.2 shows that separating orders y n could be achieved from the knowledge of all terms y n,q . To this end, a first solution has been proposed in [START_REF] Bouvier | Nonlinear Homogeneous Order Separation for Volterra Series Identification[END_REF]. It combines phase deviations (as in CPS) and amplitude gains (as in AS). In [START_REF] Bouvier | Homophase signals separation for Volterra series identification[END_REF], a notion of homophase signals is introduced and used to derive an iterative identification method. To reject numerical difficulties due to Vandermonde matrices and reserve the amplitude of inputs as a free parameter, a better alternative for order separation is proposed below. It involves linear combinations of pairs of signals with various phase deviations, according to the following theorem.

Theorem 4.4. Consider a Volterra series with a maximal order N ∈ N ∪ {+∞}. Denote θ i , i = 1, 2 two complex scalars on the unit circle, and let y ∈ Y be the output of the series excited by 2 Re [(θ 1 + θ 2 ) u] ∈ U. Then

y(t) = |p1|+|p2|≤N θ p1 1 θ p2 2 ψ (N ) p1,p2 (t) , (16)
with, ∀(p 1 , p 2 ) ∈ Z 2 , ψ p1,p2 (t) := N n=|p1|+|p2| n≡(p1+p2) [2] κ (n) |p1|,|p2| y n, n-p 1 -p 2 2 (t) , (17)
where where a ≡ b [k] stands for 'a and b are congruent modulo k' and where

κ (n) a1,a2 := n-a2 l=a1 l≡a1 [2] n l+a1 2 , n-l+a2 2 , l-a1 2 , n-l-a2 2 , ( 18
)
with n q1,q2,q3,q4 = n! q1!q2!q3!q4! the multinomial coefficient. Proof. See Appendix A.
Corollary 4.5 (Symmetries on ψ p1,p2 ). Using Property 2 and κ

(n) a1,a2 = κ (n)
a2,a1 , the following properties holds:

(i) equality: ψ p1,p2 = ψ p2,p1 (ii) conjugation: ψ -p1,-p2 = ψ p1,p2 and ψ p,-p is real.

Theorem 4.4 shows how the output of a Volterra series can be decomposed onto a 2dimensional phase plane (p 1 , p 2 ) (see Fig. 3). This makes possible to separate signals ψ p1,p2 using a 2-dimensional DFT. Furthermore, this repartition exhibits the tetrahedron structure of a Pascal's pyramid, which is due to the multilinearity of the V n operators (see Fig. 4). This property allows the extraction of signals y n,q from ψ p1,p2 . Those two steps constitute the main part of the proposed separation method, which is now detailed.

RPS Method

Denote V ∈ VS the operator of a truncated Volterra series with degree N , and u ∈ U C a complex-valued reference signal. Let u R = Re [u] ∈ U the associated real-valued excitation signal, and y n = V n [u R ] ∈ Y the homogeneous orders we want to separate.

(i) Creation of the signals collection: Let P = 2N + 1 and w = e j2π/P . Introduce the collection of input signals U k1,k2 = 2 Re w k1 + w k2 u ∈ U, k 1 , k 2 = 0, . . . , P -1, and output signals

Z k1,k2 = V[U k1,k2 ] ∈ Y.
(ii) Separation of signals ψ p1,p2 : Introduce the matrix Ψ such that ψ p1,p2 (t) = Ψ p1[P ],p2[P ] (t), ∀p 1 , p 2 = -N, . . . , N . Let F P,P denote the 2D-DFT operator 3 of order (P, P ). Then Theorem 4.4 gives

Z(t) = F P,P {Ψ} (t) . ( 19
)
3 As for the one dimensional DFT, the normalization factor is put in the inverse operator, meaning that

F M,N {x} k,l = M -1 m=0 N -1 n=0 xm,ne -j2π( km M + ln N ) ,
where x is a two-dimensional signal in matrix form. denoting the integer part. Also let

p 1 p 2 -3 - -3- -2 - -2- -1 - -1- 0 - 0- 1 - 1- 2 - 2- 3 - 3
v p = y Q, Q-p 2 y Q-2, Q-p 2 -1 . . . y |p|, |p|-p 2 T
regroup all interconjugate signals y n,q with q = n-p 2 . Then, using (17), the following relationship holds:

φ p (t) = Λ p v p (t) , (20)
where

Λ p =            κ (Q) r0,rQ 0 . . . 0 κ (Q) r1,rQ-1 κ (Q-2) r1,rQ-1 . . . 0 . . . . . . . . . . . . κ (Q) rQ/2,rQ/2 κ (Q-2) rQ/2,rQ/2 . . . κ (|p|) rQ/2,rQ/2 . . . . . . . . . . . . κ (Q) rQ,r0 0 . . . 0            (21) p 1 p 2 n 0 1 2 3 -3 -2 -1 0 1 2 3 0 1 2 3 y 0,0 y 1,1 y 1,0 y 1,0 y 2,2
y 2,0 2 y 2,0 y 2,0 2 y 2,1 4 y 2,1 y 3,3 y 3,0 3 y 3,0 3 y 3,0 y 3,0 3 y 3,1 3 y 3,2 9 y 3,2 9 y 3,1 9 y 3,1

Figure 4. Repartition of interconjugate signals yn,q up to order 3 in the 3D half-space (p 1 , p 2 , n), with p 2 ≥ 0; the gray cylinder of axis (p 1 = 1, p 2 = 0) regroups all terms constitutive of ψ 1,0 .

for p even, and

Λ p =            κ (Q) r0,rQ . . . 0 . . . . . . . . . κ (Q) r(Q-1)/2,r(Q+1)/2 . . . κ (|p|) r(Q-1)/2,r(Q+1)/2 κ (Q) r(Q+1)/2,r(Q-1)/2 . . . κ (|p|) r(Q+1)/2,r(Q-1)/2 . . . . . . . . . κ (Q) rQ,r0 . . . 0            (22)
for p odd. Inverting (20) for all p = -N, . . . , N gives all signals y n,q , n = 0, . . . , N, q = 0, . . . , n.

(iv) Reconstruction of orders y n : Once the interconjugate signals y n,q are known, homogeneous orders y n can be reconstructed using (15).

Those four steps constitute the proposed method, which will be referred as the Real-valued Phase-based Separation (RPS) method. In practice, one does not need to make P 2 measurements, due to the symmetry U k1,k2 = U k2,k1 . So the number of test signals is only K = P (P +1) 2 = (2N + 1)(N + 1) (number of terms in the triangular part of matrix Z), and the remaining signals of matrix Z are found by symmetry. Robustness of RPS method greatly depends on the numerical conditioning of matrices Λ p . This value is fixed for a given truncation order N and phase index p, and increases with N and |p|. Therefore, as for AS method, RPS robustness will decrease when N is large.

Properties of the RPS method

Zero-padding and number of signals

As in CPS method, it is possible to use a number of phases P greater than the required minimum P min = 2N + 1. As long as P ≥ P min and the truncation hypothesis is valid, retrieving signals ψ p1,p2 from Ψ is possible.

Nonlinear order aliasing

Now consider that V ∈ VS is the operator of an infinite Volterra series. As in CPS, a nonlinear order aliasing appears in step (ii). Here it is manifested by the fact that the desired signal ψ p1,p2 is mingled with all ψ i,j such that p 1 ≡ i [P ] and p 2 ≡ j [P ], deteriorating its estimation. The only way to minimize this aliasing is to choose P 2N + 1. It would also be interesting to diminish the effect of higher-order residuals y n,q , n > N in signals ψ p1,p2 . Unfortunately, the use of a contrast factor as described in § 3.2.2 is not feasible for RPS. Indeed, because ψ p1,p2 contains contributions from all orders n ≥ |p

1 | + |p 2 | such that p 1 + p 2 ≡ n [2],
the SRR gain would be 1/ρ 2 .

Number of phases and noise reduction

Suppose that each measured output Z k1,k2 is perturbed by an i.i.d. complex Gaussian noise. Identical developments as in § 3.2.3 show that the use of a 2D-DFT improves the SNR in estimation of signals ψ p1,p2 by a factor P . But the independence hypothesis between noise is false, since Z k1,k2 and Z k2,k1 comes from the same measurement (with the same noise). Similar developments as in § 3.2.3 shows that the SNR improvement has the same order of magnitude, but the error noise between signals ψ p1,p2 is not independent (i.e. the noise correlation matrix is not diagonal). This shows that, for a given maximal order N , a large value of P allows a better estimation of signals ψ p1,p2 in step (ii), and therefore of y n,q and y n in steps (iii-iv). However a large P leads to a large number of test signals K.

Amplitude of the input collection

All signals in the input collection U k1,k2 do not have the same amplitude. It greatly depends on the chosen signal u, and also, for a given u, on the factor a = |2(w k1 + w k2 )|. The maximum a max is reached for k 1 = k 2 and is equal to 4. The minimum a min is reached when the two values w k1 and w k2 are the furthest apart on the complex circle, i.e |k 1 -k 2 | = (P -1)/2 for an odd4 P . Basic calculus shows that a min = 4 cos π 2 P -1 P . Therefore the ratio between maximum and minimum amplitudes in the input collection will be closely related to cos π 2 P -1 P -1 . Therefore, increasing the number of phases P also increases the amplitude range needed for the excitation of the system5 . For example, setting P = 19 (i.e. minimum number of phases when considering a maximum order N = 9) leads to an amplitude range of ∼ 21.7 dB.

Computational complexity

Because RPS method relies only on linear algebra operations, it is easy to compute its computational complexity6 . In the following, we consider the case of a (discrete-time) signal of length L P ≥ 2N + 1 (i.e. the number of data points in the signal is much greater than the number of phases used, which is a reasonable assumption).

Step (ii) consists of the application of an inverse 2D-DFT at each time step. Using an implementation relying on one-dimensional Fast-Fourier Transforms, the 2D-DFT has a complexity of O P 2 log P 2 . Therefore the complexity of step (ii) is O 2LP 2 log P .

Step (iii) consists in solving (20) at each time step for all phases p = -N, . . . , N , i.e. 2N + 1 times. Matrices Λ p have different sizes, but Λ 0 is always the largest with size (N + 1) × ( N/2 + 1). Because L N , the pseudo-inversion of Λ p is negligible in terms of complexity compared to the products Λ † p φ p (t) for all time steps, where Λ † p denotes the pseudo-inverse of Λ p . Therefore the complexity of step (iii) is O ((2N + 1) (N + 1) ( N/2 + 1) L) = O LN 3 . The reconstruction of orders in step (iv) can be shown to have a complexity of O LN 2 , which is negligible compared to step (iii). The total complexity of RPS method is therefore O max 2LP 2 log P, LN 3 . If log P > N , this complexity is dominated by the first term and hence determined by the number of phases P ; otherwise it is cubic in the truncation order N . In both cases, complexity is linear with respect to the signal's length. In comparison, AS method has a complexity of O LN 2 . In practice, because N and P are not extremely large, computation times for RPS are very short, and close to those of AS.

Least-squares identification based on separation

This section presents, in the Least-Square (LS) formalism, how the homogeneous orders y n or the interconjugate signals y n,q can be used instead of the direct output y for discrete-time Volterra kernel estimation. The ideas behind processes (A), (B) and (C) presented in § 5.2 could easily be conveyed to other identification methods. Discrete-time signals and systems are now considered. All previous definitions and results are available by replacing T by T d := Z, N or a subset [0, L] N with L ∈ N * , and replacing integrals in (2) by discrete sums.

Least-squares formulation and combinatorial matrix

Suppose that the system-to-be-identified has a finite memory of M samples. Given L input-outputs measurements of the system to identify, relations (1) and ( 2) can be summarized in vectorial form by y = Cf , where y is the column vector of output data, f = f 0 . . . f N T a column vector regrouping all kernel coefficients, and C = C 0 . . . C N the input combinatorial matrix containing all the delayed-products of the input signal (see [START_REF] Rugh | Nonlinear system theory[END_REF] for more details about construction of this matrix). In order to reduce the number of coefficients to estimate, only non-redundant terms of the symmetric form7 of kernels h n are considered in vectors f n .

Example 5.1. For an order 2 kernel of memory length M = 3, C 2 and f 2 are given by8

C 2 =         u[0] 2 u[1] 2 u[2] 2 . . . 0 u[1]u[0] u[2]u[1] . . . 0 0 u[2]u[0] . . . 0 u[0] 2 u[1] 2 . . . 0 0 u[1]u[0] . . . 0 0 u[0] 2 . . .         T and f 2 =         h 2 [0, 0] h 2 [0, 1] h 2 [0, 2] h 2 [1, 1] h 2 [1, 2] h 2 [2, 2]        
.

Similarly, a vectorial form of ( 14) is given by y n,q = C n,q f n , where C n,q is the input combinatorial matrix corresponding to the column vector of interconjugate signal y n,q . Matrix C n,q is constructed in a similar way as C n , but conjugated terms have to be taken into account when doing symmetrization9 .

Example 5.2. For an order 2 kernel of memory length M = 3, combinatorial matrix C 2,1 is given by

C 2,1 =          |u[0]| 2 |u[1]| 2 |u[2]| 2 . . . 0 u[1]u[0]+u[1]u[0] 2 u[2]u[1]+u[2]u[1] 2 . . . 0 0 u[2]u[0]+u[0]u[2] 2 . . . 0 |u[0]| 2 |u[1]| 2 . . . 0 0 u[1]u[0]+u[1]u[0] 2 . . . 0 0 |u[0]| 2 . . .          T

Identification methods

(A) Direct estimation: Identifying kernels directly from the measured output comes down to solving y = Cf . This requires that the data length L is larger than the number of coefficients to estimate, i.e. L ≥ M +N N (see [START_REF] Mathews | How bad are Vandermonde matrices[END_REF]).

y n,n    =    C n,0 . . . C n,n    f n . (23)
Signal length must fulfill the condition L ≥ 1 n+1 M +n-1 n .

In those methods, computation time and memory usage directly depends on the size of the matrix to be inverted (and also varies with the algorithm used for solving the linear problem).

Size of matrix C grows linearly with the data length L but exponentially with the truncation order N and system memory M , which makes this method quite cumbersome (time and memory-wise) when high nonlinearities and/or long memory are required. Matrices C n have the same number of rows as C, but the number of columns, which grows with the current order n as well as memory M , is always smaller. Therefore, the inversion of C n has a computation time and memory usage reduced compared to the inversion of C. Furthermore, we have observed in practice that method (B) in total (i.e. inversions of all C n done sequentially) is computationally quite faster than (A). Each matrix C n,q has the same size as C n . But, for method (C), the matrix to invert for each order has a number of rows n + 1 times higher. Therefore the computation time and memory usage of method (C) is higher than method (B). Compared to method (A), the augmentation in the number of rows is accompanied by a reduction in the number of columns. In terms of computation time, the latter compensates the former when M is large enough (for a given N ). Furthermore, we have observed in practice that method (C) (i.e. resolution of ( 23) done sequentially for each order n) is computationally quite faster than (A).

Evaluation and comparison

This section evaluates the proposed methods using numerical simulation experiments. The code used for the separation and identification methods is available in an open-source Python toolbox10 .

Order separation evaluation

Experiment description The simulated system is a discrete-time Volterra series truncated to order 9, with a finite memory length of M = 5 samples. Kernels values are set according to the following formula:

h n [k 1 , . . . , k n ] = 10 1-n (-1) i ki e -max k 10 (24)
Those choices are arbitrarily made in order to have (a) exponentially decaying kernels, (b) amplitude of orders y n decreasing w.r.t. n and (c) a relatively low number of parameters to estimate (i.e. 2001). Real-life systems can be expected to meet criteria (a) and (b), but memory length is usually longer. AS and RPS methods are compared using the standard relative Root-Mean-Square (RMS) error value between true and estimated order, given for each order n by n = RMS (ỹ n -y n ) / RMS (y n ), where ỹn is the estimation of order y n . For RPS method, a complex white Gaussian noise of zero-mean and unit variance is used as the reference input signal u, with a length of L = 5000 samples. For AS, the signal 4 Re [u] is used as the input reference11 . A perturbation noise (also white and Gaussian) is added to the measured outputs, with a SNR of 60 dB w.r.t. to the maximum measured RMS level. For each method, 10 experiments with different realization of the signal u are done; presented results are the median over all experiments. In order to have a meaningful comparison between AS and RPS, the same number of test signals K = 190 is used for both methods. This corresponds to the required minimum number by the RPS method. Furthermore, for AS method, gains α k are chosen such that α k = (-1) k α k/2 . Value α is set so that the minimum amplitude gain corresponds to the minimum amplitude in the collection of RPS test signals (see § 4.3.4), i.e. α 0.974. This choice was empirically observed to optimize AS method.

Results Fig. 5 displays the separation error for AS and RPS methods. It shows that, for both methods, estimation of orders higher than 5 are unusable (due to a relative error exceeding 0 dB). Furthermore, error does not monotonically increases with n: even orders (6 and 8) perform better than previous odd orders (respectively 5 and 7). This is due to the choice of an odd truncation order, i.e. N = 9. Fig. 5 also shows that AS method has lower error than RPS, with a difference between ∼ 6 dB (for order 1) and ∼ 3 dB (for order 9). If a smaller number K of test signals is used for AS, or if gains α k are not properly tuned, separation error drastically increases.

Kernel identification evaluation

Experiment description The same system and experiment as before is used. The error for kernel estimation is given by ε n = RMS hn -h n / RMS (h n ), where hn is the estimation of order h n . Least-Squares identification methods are tested and compared, using either:

• method (A) on output y corresponding to the reference signal;

• method (B) on orders y n estimated via AS;

• method (B) on orders y n estimated via RPS;

• method (C) on signals y n,q estimated via RPS, i.e. RPS without step (iv).

The direct approach relies only on one measurement of L samples, and is therefore a far quicker experiment than using separation methods. For a meaningful comparison, a measurement of KL samples should be used; but this would drastically increase the size of the ), orders yn estimated via RPS method ( ) and terms yn,q estimated via RPS method ( ). Median value for all experiments are represented. combinatorial matrix C, and consequently the memory and CPU usage12 . For application on real-life systems with longer memory M , this approach would lead to extremely large matrix sizes and intractable computation.

Results Estimation errors are shown in Fig. 6. As for the separation, the same non-monotonic behaviour is observed, which is still due to the choice of the truncation. Direct identification (A) has an error similar to the SNR for order 1; when n increases, this error grows rapidly. Previous observations on AS and RPS separation error still apply for identification error of both (B)-AS and (B)-RPS. In comparison to (A), those methods have lower error for kernels of orders 1 to 4, but worse estimation for orders 7 to 9. Using interconjugate signals instead of homogeneous orders improves greatly the estimation of Volterra kernels. Method (C)-RPS is ∼ 20 dB lower than (A) for all orders. In comparison to (B)-RPS, this method ameliorates identification, especially for higher orders: this is the only method where identification results are lower than 0 dB for orders higher than 7. This is explained by the fact that, for the same number of coefficients to estimate, more data is available when using signals y n,q (i.e. (n + 1)L instead of L).

Conclusions and perspectives

In this paper, the exploitation of phase deviation in Volterra series output for separating homogeneous orders has been studied. These results have then been used to devise new kernel identification methods in the Least-Squares formalism. Numerical experiments have been carried out to illustrate the amelioration brought in terms of estimation error. The concept of interconjugate signals has been introduced. This notion permits, when working with a real-valued input, to consider the effect of the Volterra series on the complex-valued analytic signal. Using this concept, a new separation method has been proposed, relying on the 2-dimensional Discrete Fourier Transform. Furthermore, the use of interconjugate signals has been shown to improve kernel identification for higher orders. The new separation method presented in this paper has been successfully applied to a reallife system [START_REF] Bouvier | Identification de systèmes non linéaires représentés en séries de volterra: applications aux systèmes sonores[END_REF][START_REF] Lebrun | Estimation de paramètres d'un hautparleur électrodynamique non linéaire[END_REF]. Future work will focus on the use of this prior separation stage for identification of more structured systems, such as block-structured or polynomial nonlinear state-space systems.

Figure 1 .

 1 Figure 1. Identification process of kernels {hn} n=1,...,N without (a) and with (b) prior order separation

Figure 3 .

 3 Figure 3. Repartition of interconjugate signals yn,q up to order 3 in the 2D phase plane (p 1 , p 2 ); the dashed line ( ) corresponds to the symmetry axis, and the dotted line () to the hermitian symmetry axis; e.g. the gray circle corresponds to the location of ψ 1,0 .

(

  B) Estimation on homogeneous orders: If the homogeneous orders y n are available, kernel estimation is done by solving y n = C n f n for each order n. Signal length must fulfill the condition L ≥ M +n-1 n . (C) Estimation on interconjugate signals: If the interconjugate signals y n,q are available, kernel estimation is done by solving, for each order n,

Figure 5 .

 5 Figure 5. Comparison of separation error w.r.t order n for methods AS ( ) and RPS ( ). Median values over all experiments are represented.

Figure 6 .

 6 Figure 6. Comparison of identification error w.r.t order n for LS identification using output y ( ), orders yn estimated via AS method (), orders yn estimated via RPS method ( ) and terms yn,q estimated via RPS method ( ). Median value for all experiments are represented.

The links between the notion of order separation and order orthogonality (which is the key basis of the Wiener series framework) are not addressed here and deserve further study; for identification methods using orthogonalization approaches, see(Kibangou & Favier

, 2010;[START_REF] Tseng | Identification of cubic systems using higher order moments of iid signals[END_REF]. 2 Available at https://github.com/d-bouvier/pyvi

in only one work[START_REF] Zhang | Volterra series truncation and kernel estimation of nonlinear systems in the frequency domain[END_REF]. To design a robust identification, the main idea of this paper is to exploit phase deconstruction/reconstruction in a preliminary order separation step.

An even P leads to a min = 0 and must be avoided.

It is important to note that using RPS, the amplitude variations between all signals of the input collection are not wanted and chosen (as in AS method) but are a side effect of the method.

The memory usage of RPS method will not be addressed here; it is generally not an issue, and computation can be done separately for each time step.

Or equivalently the nonzero terms of the triangular form.

The lexicographical order is used here, but other choices can be made, without impact on the identification.

For example, combinatorial terms u[l]u[l -1] and u[l -1]u[l] are both linked to the same kernel coefficient h 2 [0, 1], and should thus contribute equally to the corresponding term in C 2,1 .

https://github.com/d-bouvier/pyvi

The factor 4 is here so that both collection of test signals have the same maximum amplitude.

For example, in the limited case of this experiment, increasing the number of samples from L to KL would make the memory usage to store the matrix C alone going from 76 MiB to 14.5 GiB.

Appendix A. Proof of Theorem 4.4

Consider the hypothesis of Theorem 4.4. Because 2 Re [(θ 1 + θ 2 ) u] = θ 1 u + θ 1 u + θ 2 u + θ 2 u, the n-linearity and symmetry of V n gives the following expression for the output of order n:

Consider the change of variables

Then it is straightforward to show that the new variables l i , p i , i = 1, 2 respect the following conditions

Applying the change of variables (A2) to (A1) gives

which can be simplified as

The second summation is equal to κ

|p1|,|p2| as defined in (18). This is easily seen by observing that the multinomial coefficient is invariant under any sign change of p 1 and/or p 2 . The total output is then

The two sums can be swapped as follows

which concludes the proof.