Kamardine Selemani

Elodie Richalot

Olivier Legrand

Odile Picon

F Mortessagne

Energy localization effects within a

Keywords: Chaotic cavity, field homogeneity, field localization, reverberation chamber, statistical distribution

I. INTRODUCTION

I N electromagnetic compatibility (EMC) tests using a re- verberation chamber (RC), the device under test (DUT) is placed within the so-called working volume [START_REF]Electromagnetic Compatibility (EMC)-Part 4-21: Testing and Measurement Techniques-Reverberation chamber test methods, International Electrotechnical Commission[END_REF], situated in the cavity central volume, in which a uniform and isotropic field distribution is required. Since only the field in the working volume is useful to illuminate the DUT, it is important to optimize the power transmission towards this area, in particular when high field levels are required to test the device immunity. The understanding of wave function structure in complex geometries where the ray dynamics is nonintegrable has been a major interest of Wave Chaos research. Very early, the existence of large fluctuations in wave function intensity, with a large excess of intensity near the short unstable periodic orbits, was observed [START_REF] Mcdonald | Wave chaos in the stadium: statistical properties of short-wave solutions of the Helmholtz equation[END_REF] [3] [START_REF] Bogomolny | Smoothed wave functions of chaotic quantum systems[END_REF] and coined scarring. In chaotic cavities, this phenomenon was clearly identified as the main nonrandom features of chaotic wave functions. In systems where the dynamics of rays are not fully chaotic, wavefunctions may display localization effects related to trapping on regular islands [START_REF] Kim | Quantum-classical correspondence in localization of eigenstates for a system having mixedtype classical dynamics[END_REF].

Due to the strongly scattering nature of a classical mode stirrer and the multiple scattering environment created by the resonant cavity, such a phenomenon is likely to appear in an RC. The large dimensions of the stirrer in terms of the wavelengths under consideration precludes the appearance of strong localization effects, but a weak field enhancement can be envisaged. In this case, the field would no longer be homogeneous over the whole cavity, even if the homogeneity could still be observed in a restricted zone. For this reason, we investigate in this paper the field localization within an RC equipped with a mode stirrer, and propose a new cavity shape designed to avoid such a phenomenon, where the key ingredient is to replace the classical stirrer by a moving metallic hemisphere. The local field enhancement phenomenon is brought to light in two ways. Our first approach consists in calculating an indicator of the existence of an electric or magnetic energy enhancement area in the total volume of both cavities; it is called the Inverse Partition Ratio and has been widely used to quantify wave function localization either in chaotic or disordered wave systems [START_REF] Kaplan | Short-time effects on eigenstate structure in Sinai billiards and related systems[END_REF] [7] [START_REF] Pradhan | Correlations due to localization in quantum eigenfunctions of disordered microwave cavities[END_REF]. In our second approach, dedicated to the research of the field enhancement areas, we divide each cavity into sub-volumes, and compare the mean energies in the different sub-volumes. We will see that the localization effect is principally due to the presence of the mode stirrer and restricted to its vicinities, and is drastically reduced in our proposed alternative RC shape. Both approaches are applied at a fixed position of the stirrer and hemipshere as well as for different positions of them over a whole rotation.

It will be shown how such a field localization phenomenon affects the field statistical uniformity and isotropy. As a consequence, an improvement of the field properties is obtained when this phenomenon is avoided.

The improvement of the RC behavior with the suggested cavity shape modifications is finally shown with regard to standardized criteria from the norm governing EMC test in reverberation chambers [START_REF]Electromagnetic Compatibility (EMC)-Part 4-21: Testing and Measurement Techniques-Reverberation chamber test methods, International Electrotechnical Commission[END_REF].

II. STUDIED MICROWAVE CAVITIES

In order to illustrate simply the appearance of a localization phenomenon of electric or magnetic energy in a reverberant room, we consider in the first instance a 2D rectangular cavity provided with a metallic mode stirrer of complex shape (Fig. 1). The cavity is of length 756 mm and width 456 mm, walls and stirrer are of perfect metal. Using HFSS software for electromagnetic simulations, the eigenmodes of the cavity have been determined. Fig. 1 shows the distribution of electric and magnetic field amplitudes for the 53 rd eigenmode at 2.3873 GHz. We observe higher energy densities in the central area of the cavity (the useful zone), wherein the DUT could be placed. Thus this mode will contribute to achieve the required energy level in the working area. Indeed, the mean electric and magnetic energies calculated within the working area are higher than the ones calculated in the rest of the cavity volume: the ratio of the mean energies calculated in the working volume and outside is of 3.42 for the electric field and 3.03 for the magnetic one. A different behavior is observed with the 44 th eigenmode at 2.2087GHz depicted in Fig. 2. In this case, the electric and magnetic energy densities are higher in the area situated between the mode stirrer and one of the closest walls, so that the effect of this mode in the working area is reduced. In this case, the ratio of the mean energies calculated in the working volume and outside is of 0.655 for the electric field and 0.608 for the magnetic one.

This simple case highlights an energy localization phenomenon in a reverberation cavity disturbed by a mode stirrer. To examine its effect on the functioning of a reverberation chamber, this phenomenon will now be examined in 3D cavities and the effect of the mode stirrer displacement will be taken into acount.

Two parallelepipedic cavities of dimensions W = 0.785 m along (Ox), L = 0.985 m along (Oy) and H = 0.995 m along (Oz) are considered (Fig. 3). The classical RC (C1, Fig. 3a) consists of a parallelepipedic cavity equipped with a mode stirrer. The shape and location of the latter conform to an industrial RC, except for a global scaling factor. The metallic paddles are 3 mm thick whereas the vertical circular rod at the panels rotation axis has a diameter of 8 mm. All surfaces are considered as perfectly conducting. The previous 2D example suggests that the presence of a mode stirrer within the cavity volume could induce a field localization in its vicinity. Therefore we propose an alternative RC geometry with a hemispheric mode stirrer located on a cavity wall (C2, Fig. 3b). Apart from avoiding localization effects, the chosen geometry also has to lead to field uniformity and isotropy. As in chaotic cavities most of the resonant modes present uniform and isotropic field distributions [START_REF] Gros | Universal behaviour of a wave chaos based electromagnetic reverberation chamber[END_REF] [10] [START_REF] Stöckmann | Chaos: an introduction[END_REF], we draw inspiration from an already studied 2D chaotic cavity [START_REF] Méndez-Sànchez | Distribution of Reflection Coefficients in Absorbing Chaotic Microwave Cavities[END_REF] to propose a new 3D RC geometry. It consists of a parallelepipedic cavity provided with two metallic spherical caps and one metallic spherical hemisphere on its walls and in electrical contact with them, a spherical cap being a portion of a sphere cut off by a plane and whose height is, in our case, less than the sphere radius for a cap or equal for a hemisphere. In previous works we already studied this resonant cavity and showed its chaotic behaviour (with spatial and spectral statistics in agreement with those predicted by Random Matrix Theory [START_REF] Mehta | Random Matrices[END_REF]) as well as the improvement of the field properties with regard to the EMC standards in comparison to a classical RC [START_REF] Gros | Universal behaviour of a wave chaos based electromagnetic reverberation chamber[END_REF] [14] [START_REF] Richalot | Increasing contribution of stochastic methods in the characterization of complex systems: reverberation chamber and urban environment[END_REF]. The hemisphere has a radius of 15 cm and caps of 45 cm and 50 cm. The highest penetration depth of the caps within the cavity is of 15 cm. Whereas the two caps are fixed, the stirring process is ensured by the translation of the hemisphere. Its center rotates on the x=W plane along a circle of radius 27.85 cm and centered at (W; L/2.1; H/2). The initial position of the hemisphere center, used in the following while studying the cavity in a fixed configuration, is (W; 25.5 cm; 32 cm).

The two cavities are simulated, without considering any excitation antenna, by using the eigenmode solver of HFSS software. Simulations are performed without considering any losses in the cavity, so that the electric field components are real. For both cavities, simulations are first of all performed for a fixed position of the stirrer/hemisphere, then for a whole stirrer/hemisphere rotation.

III. EIGENMODES PROPERTIES A. For a fixed stirrer position

Originally proposed in quantum mechanics to distinguish the energy regions of localized states of a particle [START_REF] Wegner | Inverse Participation Ratio in 2+ε Dimensions[END_REF] [17], the Inverse Partition Ratio (IPR) has also been used in microwave cavities to evaluate the degree of localization of the electric field [START_REF] Pradhan | From chaos to disorder: Statistics of the eigenfunctions of microwave cavities[END_REF]. This indicator is calculated as :

L A k = V |A k | 4 dv V |A k | 2 dv 2 • V (1)
with k = x, y or z, A is the complex amplitude of electric or magnetic field. For a real field component following a normal distribution this indicator is equal to 3 (4.77 dB), whereas it increases with a localization phenomenon. A global indicator of the field localization is calculated as the quadratic mean of the three IPR related to field components :

L (E,H) tot = (L (E,H) x ) 2 + (L (E,H) y ) 2 + (L (E,H) z ) 2 3
(2)

Figure 4 shows the IPR of three rectangular electric field components as well as the global ones for the 330 first eigenmodes respectively in the chaotic cavity and in the classical cavity. We observe much higher IPR values attained in the classical RC compared to the chaotic one, the IPR values remaining close to the ideal value in the latter. It indicates a strong localization of the electric energy for many eigenmodes in C1. IPR values are particularly high at low frequencies, with numerous peaks above 20 dB below 0.564 GHz (corresponding to the 42 nd mode), but the localization effect remains at high frequencies with peaks above 15dB even after the 300 th mode (1.08 GHz). In the chaotic cavity (Fig. 4b), the maximal total IPR value is of 10 dB (only attained once) and IPR values remain below 6 dB after the 32 nd mode (0.541 GHz). As a comparison, the global IPR index exceeds 10 dB for 86 modes in the classical cavity. Thus it clearly appears that numerous modes of the classical RC undergo a strong electric energy localization; these eigenmodes present a heterogeneous distribution. In order to better quantify this phenomenon, Table I presents the means and standard deviations over the whole studied frequency band of the four examined indicators calculated from electric field distribution in each cavity. Due to the improvement of the energy uniformity with the frequency, these quantities are given while considering the 330 first modes then after excluding the first 30 and 60 ones. The resonant frequencies of the 30 th and 60 th modes are of 0.507 GHz and 0.632 GHz in C1, and 0.527 GHz and 0.666 GHz in C2. Table I confirms the previous observations. The mean IPR values close to the ideal one combined to their low frequential variation in the case of the chaotic cavity indicate a spatial and spectral uniformity of the electric field distribution. The much higher values of the means and standard deviations in the classical RC reflect the strong localization phenomenon suffered by numerous modes. By suppressing the first 30 then 60 modes, the standard deviations decrease in the chaotic cavity, whereas the means and standard deviations improve in the classical cavity. The performances of the chaotic cavity still remain far better than these of the classical RC.

I L E k L E x L E y L E z L E tot µ σ µ σ µ σ µ σ C1 9.
So far, only the electric energy distribution has been considered. To verify that the observation of the magnetic energy distribution would lead to the same conclusions, the same indicators have been calculated using the IPR definition based on the magnetic field component amplitudes (Eq.1). A similar variation of both indicators, based on electric or magnetic fields, is observed. This similar behavior will allow us later on to restrict our study to the electric field distribution. We have highlighted a localization phenomenon in the classical RC, but IPR index doesn't indicate the region of higher energy density. To answer this question, we divide the studied cavities into several non-overlapping sub-volumes and compare their energy densities. In the classical RC (Fig. 5a), a volume is defined around the stirrer (of size 27.2 cm×32.2 cm along x-and y-axes) and another in the central area (at 50 mm from the cavity walls). The chaotic cavity is subdivided into four peripheral domains and a central one (at 50 mm from the cavity walls and caps) (Fig. 5b). We then evaluate the ratio between the mean electric or magnetic energy densities within the volume around the stirrer V b and the central volume V c in the classical cavity and the ratios between the mean energy densities within each peripheral volume V i and the central volume V c in the chaotic cavity.

R A i = Vi |A| 2 dv V i • V c Vc |A| 2 dv (3)
with V i the peripheral volume (V b in the conventional RC), V c the central volume, and A the vectorial complex amplitude of electric or magnetic field. A uniform field distribution corresponds to a unit energy ratio, whereas a deviation from unity indicates a higher energy concentration in one of the defined volumes. Fig. 6a representing the ratio of the mean energy around the stirrer and the one in the central volume of the conventional cavity shows a wide discrepancy with the ideal value, with numerous high ratio values (57 peaks of R E b among 475 over 2). In contrast, in Fig. 6b, showing the ratios of mean electric energies in the four peripheral volumes of the chaotic cavity and the one in the central volume, all the curves remain close to unity on the whole frequency band, with only three peaks over 2; it indicates a uniform energy distribution between the five defined sub-domains. As already observed with the IPR index, the comparison of the ratios calculated using the electric and magnetic energies (Fig. 6a) shows similar behaviors, so that the study can be restricted to the observation of only one of them. We have chosen to focus on the electric energy distribution. Table II summarizes the properties of the mean electric energy ratios in both cavities, either giving their mean and standard deviations on the whole frequency band or by excluding 30 and 60 modes. Whereas the means of the ratios are close to unity in the chaotic cavity, they are higher in the classical RC indicating an energy localization around the stirrer at the expense of the energy level required for EMC tests in the working area. In the chaotic cavity the ratios display low variation versus frequency, thus most of the modes ensure a uniform energy distribution between the five sub-volumes. On the contrary, the higher values of the ratio standard deviations in the conventional RC denote the spectral variability of the field spatial distribution properties.

As expected, all the mean values come closer to one and all the standard deviations decrease when the first 30 then 60 modes are excluded, but the performances of the chaotic RC remain better than those of the conventional RC.

To illustrate this phenomenon of energy localization around the mode stirrer, Fig. 7 depicts electric field amplitudes of a localized mode, namely the 205 th one (956 MHz,

L E tot = 16.238, L H tot = 22.211, R E b = 2.847, R H b = 3.119
). The localization of the electric energies around the mode stirrer is obvious. The field amplitudes are low within the cavity working volume, so that this mode contributes very little to the illumination of the DUT. It is also to notice that IPR and energy ratios are in agreement with field distributions, as they all indicate that energies are localized.

B. Effect of the stirrer rotation

In this section, we evaluate the previously introduced indicators of localization in the two cavities for different positions of the stirrer and the hemisphere. This study aims to verify that the localization phenomenon observed for fixed cavity geometries is not bounded to a specific configuration but occurs at several stirrer or hemisphere positions. This study has been carried out on three frequency ranges above the 60 th mode while considering 15 modes above 700 MHz, 10 modes above 800 MHz and 6 modes above 1.2 GHz. As we already mentioned the similar behaviors of the indicators based on electric and magnetic fields, only the electric energy distribution will be considered here for 360 positions of the stirrer and of the hemisphere on a whole rotation (for a rotation angular step of 1 • ).

For each considered mode of the classical RC, Fig. 8 indicates the variation of the related global IPR indicator L E tot versus the stirrer (Fig. 8a-c) or hemisphere position (Fig. 8df). Its mean and standard deviation over all the stirrer\ hemisphere angular positions are indicated in Fig. 9.

Fig. 9 shows a decrease with the frequency range of L E tot mean values over the stirrer or hemisphere rotation, corresponding to a weakening of the localization phenomenon. For the chaotic cavity, mean values close to the ideal value (3 in linear or 4.77 dB) are obtained on the two highest frequency bands. The decrease of L E tot standard deviations with the frequency is also obvious; thus fewer modes suffer from high localization effect. However, even if the same trends are observed for both cavities, on a given frequency band the performances of the chaotic cavity are always better than those of the classical RC, with mean values of L E tot closer to 3 and lower standard deviations. Thus, in the chaotic cavity, the eigenfield distributions are more uniform and less sensitive to the stirring object position than in the classical RC. From L E tot mean values over the stirrer rotation calculated for the different modes of each frequency band, a mean value over the modes is calculated. In order to quantify the behavior differences between the modes of a given frequency band, the standard deviation of the same ensemble L E tot mean values is also determined over the modes of each studied frequency band. The obtained means and standard deviations are given in Tab.III. These results highlight the better field uniformity in the chaotic cavity with mean values on the three frequency bands lower than the one attained on the highest frequency band in the conventional RC. Moreover, the much lower standard deviations indicate a better spectral uniformity obtained in this alternative RC. The improvement of the cavity properties while increasing the frequency is also noticeable with these results. In particular, the mean value attains the ideal one (namely 3) in the highest frequency band for the chaotic cavity.

IV. TOTAL FIELD PROPERTIES

A. Field reconstruction

Uniformity and anisotropy coefficients have been defined in the norm [START_REF]Electromagnetic Compatibility (EMC)-Part 4-21: Testing and Measurement Techniques-Reverberation chamber test methods, International Electrotechnical Commission[END_REF] to be applied to measurement results. They are calculated from the electric field within the excited RC and not the previously examined eigenmodes. Before applying them in order to compare the performances of both cavities regarding standardized criteria, the excited field has to be rebuilt from the numerically determined eigenmodes. As the eigenmodes have been numerically determined without considering any losses, this approach is approximate but can be considered as valid for low losses.

For a cavity excited at the angular frequency ω, the electric field at the same frequency can be expanded on the cavity eigenmodes as [START_REF] Cerri | Source stirring mode for reverberation chambers[END_REF]:

- → E (ω) = - ωµ V ∞ m=1 1 k 2 m -k 2 V J( r 0 ) • e m ( r 0 ) dv V | e m | 2 dv • e m (4) where k 2 = k 2 0 [1 -(-1 + ) • ω m / (Q m ω)],
with Q m the quality factor of the m th eigenmode. e m is the m th eigenmode field whose mean of the square electric field amplitude on the whole cavity volume V is normalized to one. k m is related to the m th resonant angular frequency ω m by k m = ω m /c, with c the speed of light in vacuum. The current J represents the source and r 0 its location in the cavity.

The quality factors Q m are calculated by considering the effects of the Joule losses on the cavity walls and of the cavity loading by an antenna, according to

(Q m ) -1 = (Q Joule m ) -1 + (Q ant m ) -1 . The Q Joule m
and Q ant m are calculated in the high frequency approximation [START_REF] Dunn | Local, High-Frequency Analysis of the Fields in a Mode-Stirred Chamber[END_REF] [21] (Eq. 5).

Q Joule m = 3 V 2 S δ and Q ant m = 16 π 2 V f c 3 (5)
with S the total surface area of the cavity walls and δ the wall skin depth.

In order to ensure the obtained results are not specific to a given excitation location, eight pointlike current sources are considered within both cavity volumes (Fig. 3), these unit current sources having an amplitude of 1/ √ 3 along each axis of the Cartesian coordinate system in order to excite similarly the three field components.

B. Impact of the localization to the eigenmodes contributions

The previous study of the eigenmodes properties showed that in the classical RC numerous modes have their energy highly concentrated around the stirrer to the detriment of the energy in the working area. According to Eq.4, the electric field within the working area is the sum of the contribution of each eigenmode; if an eigenfield is mostly concentrated around the stirrer, its contribution to the total field within the working volume will be reduced. To show the impact of the eigenfield localization outside the useful area on its contribution to the total field in the useful volume, we define the modal contribution E m of the m th mode as the excited electric field obtained while considering only the m th mode, without the contribution of the other modes. Derived from Eq.4, the modal contribution E m of the m th mode can be written :

- → E m(ωm) = - ωµ V (k 2 m -k 2 ) V J( r 0 ) • e m ( r 0 ) dv V | e m | 2 dv • e m (6)
and its amplitude is maximal at the m th mode resonant angular frequency ω m . From the extraction of the three electric field components of the m th eigenmode at the center of the working volume, we calculate the contribution at this central point of the first 10 eigenmodes above 700 MHz (Fig. 10.a for the conventional RC and Fig. 10.c for the chaotic one) while considering the S 1 excitation location in both cavities. The total electric field due to the contribution of these 10 modes and calculated using Eq. 4 is presented in Fig. 10.b for the conventional RC and Fig. 10.d for the chaotic one.

In Fig. 10.a, representing the evolution of the modal frequency weight |E m (ω)| of 10 modes of the classical RC, one distinguishes only 6 peaks, 4 invisible peaks corresponding to maximum amplitudes |E m (ω m )| below 0.4V/m or 2.6% of the largest peak (of 15.2V/m); these low weight modes weakly contribute to the total field and modal overlap (Fig. 10.b). For the chaotic cavity in contrast, the lowest peak amplitude is of 1.34V/m or 8.1% of the largest peak (of 16.5V/m). Thus all peaks are visible in Fig. 10.c and the contribution of the 10 modes is found on the total field (Fig. 10.d). The greatest disparity between the peak amplitudes for the different modes in the conventional cavity results in a normalized standard deviation between the amplitudes of 1.25 against 0.69 in the chaotic cavity. Let us now relate these modal contributions and localization indicators previously calculated. For the 10 modes considered, the mean values and standard deviations of L E tot are 14.6 and 17.12 in the conventional cavity against 3.95 and 0.25 in the chaotic cavity. Indeed, no mode presents a high energy localization in the chaotic cavity, contrary to the conventional cavity. In the latter, the most important localization effect is obtained for the 7 th mode (R b = 3.05, L E tot = 61.31) associated with an invisible amplitude peak (0.147 V/m). In contrast, the highest amplitude peak of the modal contributions (15.2 V/m for mode 3) is associated with parameters (R b = 0.43, L E tot = 4.34) indicating no energy localization around the stirrer. This shows the negative impact of energy localization patterns on spectral uniformity of the total field and on the modal overlap, the latter playing an important role in obtaining a uniform and homogeneous field in the RC [START_REF] Gradoni | Generalized extreme-value distributions of power near a boundary inside electromagnetic reverberation chambers[END_REF].

C. Field homogeneity and isotropy

According to Eq.4, the spectral bandwidth of each resonance frequency is non zero. As a consequence, the eigenmodes are not only excited at their resonance frequency, and even if their weight rapidly decreases outside the resonance width, modal overlap occurs. It implies that the properties of the excited field are not strictly equivalent to those of the isolated eigenmodes, and that the rebuilt field of Eq.4 has to be used to apply uniformity and anisotropy standardized criteria. In the following, the rebuilt field is calculated for frequency steps chosen equal to the third of the smallest frequential mode bandwidth in each frequency band, respectively of 14.67 kHz, 12.22 kHz and 5.91 kHz for the chaotic cavity and 13.74 kHz, 11.01 kHz and 5.45 kHz for the classical RC.

To characterize the field uniformity, four standard deviations are used. Three of them, associated to the orthogonal field components E i (i = x, y or z in our case), are calculated as [START_REF]Electromagnetic Compatibility (EMC)-Part 4-21: Testing and Measurement Techniques-Reverberation chamber test methods, International Electrotechnical Commission[END_REF]

σ i,dB = 20 • log 1 + σ i E i max (7)
where the standard deviation σ i and mean value E i max are calculated from the maximal field amplitudes attained over a stirrer or hemisphere rotation at each of the eight vertices of the working volume, the latter one being chosen to present the same volume in both cavities. A global uniformity indicator σ tot is also calculated using the maximal field amplitudes obtained by considering the three field components. Only independent stirrer/hemisphere positions are considered, i. e. 20, 22 and 32 positions on the three frequency bands for the chaotic cavity, and 18, 20 and 25 positions in the classical one, as indicated by correlation coefficients between normalized eigenfields. Fig. 11 shows the four standard deviations obtained in three bandwidths for both cavities while considering the excitation S 1 (indicated in Fig. 3). The standard deviations are lower in the chaotic cavity and seldom reach 3 dB. According to the norm, the total standard deviation has to be below 3 dB with a tolerance of a maximum of three peaks above 3 dB per octave. In the chaotic cavity, the standard deviations stay below the 3dB limit in the two higher frequency bands, whereas it is only respected in the highest band in the classical RC.

The frequential averages of the four curves of Fig. 11 are given in Tab.IV. They highlight the lower standard deviations in the chaotic cavity, as all the average values are lower in this cavity; it indicates that a better field uniformity is obtained. We also notice in this cavity lower differences between the mean standard deviations related to the three field components; it shows a better field isotropy. These mean values also confirm the decrease of the standard deviations with the frequency, but the properties of C2 always stay better than those of C1. Using the same procedure for the field reconstruction at the eight vertices of the working volume and the stirring process, we also determined the field anisotropy coefficients (A αβ ) given by [START_REF]Electromagnetic Compatibility (EMC)-Part 4-21: Testing and Measurement Techniques-Reverberation chamber test methods, International Electrotechnical Commission[END_REF]:

A αβ = |E α | 2 -|E β | 2 |E α | 2 + |E β | 2 (8)
where α and β indicate directions of mutually orthogonal Cartesian components, and the average is performed firstly on the hemipshere or stirrer positions then on the eight field extraction points. From these three partial anisotropy coefficients, a global anisotropy indicator A tot can also be calculated using :

A tot = A 2 xy + A 2 yz + A 2 zx /3 (9)
Fig. 12 shows the amplitude of these four electric field anisotropy coefficients on the three frequency bands in both cavities.

First we note that the anisotropy coefficients decrease with frequency in both cavities. In all frequency bands, the anisotropy coefficients are lower in the chaotic cavity, indicating a better electric field isotropy. According to the standard [START_REF]Electromagnetic Compatibility (EMC)-Part 4-21: Testing and Measurement Techniques-Reverberation chamber test methods, International Electrotechnical Commission[END_REF], with a number of independent positions between 10 and 30 as in our case, a good stirring quality is characterized by a total anisotropy coefficient less than or equal to -5 dB, whereas a medium stirring quality is associated to the limit -2.5 dB. In the chaotic cavity, the total anisotropy coefficient remains lower than -2.5 dB in the three frequency bands and decreases with frequency, whereas in the cavity with a stirrer this criterion is only respected in the band above 1.2 GHz. However, the standard is respected below the -2.5 dB limit, and this is achieved on the highest frequency band for both cavities.

Based on these values and according to the norm [START_REF]Electromagnetic Compatibility (EMC)-Part 4-21: Testing and Measurement Techniques-Reverberation chamber test methods, International Electrotechnical Commission[END_REF], the chaotic cavity permits to obtain a good stirring quality (total anisotropy coefficient below -5dB) whereas the conventionnal RC can only provide a medium stirring quality.

To confirm these remarks, the amplitudes of the four anisotropy coefficients have been averaged on each frequency band (Tab.V). We notice the decrease of these values with the frequency, and the best performances of the chaotic cavity, with lower mean values. We have shown that, for a point like source located at S 1 , the electric field over a stirrer or hemisphere rotation is more uniform and isotropic in the chaotic cavity. We will now verify that this conclusion is independent of the excitation location.

Tab.VI shows the frequency means on the three frequency bands of the total standard deviation σ tot and the total anisotropy coefficient A tot for the seven other source locations S i presented in Fig. 3. Similar results are obtained at all the source locations, and the conclusions drawn at the first source location remain valid: the electric field is always more uniform and isotropic in the chaotic cavity and its properties improve in both cavities with frequency increasing. The latter results allow us to confirm the better properties of the chaotic cavity regardless of the source position.

V. CONCLUSION

It has been demonstrated that an energy localization phenomenon appears around the mode stirrer of a classical RC even at high frequencies. As the resonant modes undergoing this localization weakly contribute to the total field within the working volume, it degrades the spatial and spectral properties of the field in this useful area. To overcome this problem, a new RC shape has been proposed, that has been obtained through simple modifications of the conventional RC. The localization phenomenon is avoided in this case, so that the energy distribution is more uniform within the whole cavity and the modes contribution to the field within the working volume is stronger and less frequency dependent. The resulting improvements of the field properties are clear by considering standard uniformity and isotropy criteria, as the well operating indicators move below the required threshold at lower frequency and, on a given frequency band, always indicate a better operating of the chaotic cavity. The proposed simple and low-cost modifications of classical RC could thus permit to improve the field properties within a classical RC and also to decrease its Lowest Useable Frequency.

  (a) Electric field amplitude (b) Magnetic field amplitude Fig. 1. Normalized amplitudes of electric and magnetic fields at f = 2.3873 GHz for a 2D cavity with a metallic stirrer (in red). The red line indicates the working area situated at a quarter-wavelength from the walls and stirrer.

  (a) Electric field amplitude (b) Magnetic field amplitude Fig. 2. Normalized electric and magnetic field amplitudes at f = 2.2087 GHz for a 2D cavity with a metallic stirrer (in red).

  (a) Classical RC (C1) (b) Chaotic cavity (C2) Fig. 3. Chaotic and classical cavities with 8 points sources S i .

Fig. 4 .

 4 Fig. 4. Electric field localization index (dB), (a) in the classical RC, (b) in the chaotic cavity (blue for x-, red for y-and green for z-component and black for L E tot ).

Fig. 5 .

 5 Fig. 5. Integration volumes (a) within the classical RC : volumes V b and Vz, and (b) within the chaotic cavity: peripheral volumes (V 1 , V 2 , V 3 , V 4 ) and central volume Vc.

Fig. 6 .

 6 Fig. 6. Mean electric energy ratios for R E i in (a) C1, (b) C2.

Fig. 7 .

 7 Fig. 7. Electric field amplitude of the 205 th mode (f=908 MHz) represented by iso-amplitude surfaces. Only the amplitudes above an arbitrarily fixed limit are represented.

8 .

 8 Global L E tot (in dB) versus stirrer/hemisphere position (whole rotation with 1 • step) in conventional (a-c) and chaotic (d-f) cavities, (a)&(d) 15 modes above 700 MHz, (b)&(e) 10 modes above 800 MHz and (c)&(f) 6 modes above 1.2 GHz.

Fig. 9 .

 9 Fig. 9. Means (blue) and standard deviations (green) over stirrer/hemisphere rotation of the global L E tot (linear scale), in conventional (a-c) and chaotic (d-f) cavities, (a)&(d) 15 modes above 700 MHz, (b)&(e) 10 modes above 800 MHz and (c)&(f) 6 modes above 1.2 GHz. The values for each mode are indicated at its resonant frequency at the initial stirrer/hemisphere position.

Fig. 10 .

 10 Fig. 10. Amplitude of the modal contribution |E m(ω) | at the center of working volume for the eigenmodes m i ((a) for classical cavity and (c) for chaotic cavity); and amplitude of the total electric field induced by their contribution ((b) for classical cavity and (d) for chaotic cavity).

Fig. 11 .

 11 Fig. 11. Standard deviations of the three Cartesian field components and total one (blue for σx, red for σy, green for σz, and black for σtot) for (a 1 , a 2 , a 3 ) cavity C 1 and (b 1 , b 2 , b 3 ) cavity C 2 in the three frequency bands for a source located in (S 1 ).

TABLE

  Frequency means on the three frequency bands of the four standard deviations and maximal difference between the three ones related to field components for both cavities, source at S 1 .

			TABLE IV
			σx	σy	σz	σtot	δσ k
	700M Hz	C1 C2	2.99 1.98	2.57 2.05	2.90 2.15	3.61 2.37	0.42 0.17
	800M Hz	C1 C2	2.18 1.76	2.99 1.69	2.64 1.89	3.32 1.90	0.81 0.20
	1200M Hz	C1 C2	1.90 1.48	1.63 1.57	1.58 1.54	1.94 1.63	0.32 0.09

  Fig.[START_REF] Méndez-Sànchez | Distribution of Reflection Coefficients in Absorbing Chaotic Microwave Cavities[END_REF]. Amplitude of the field anisotropy coefficients (in dB) versus frequency (blue for Axy, red for Ayz, green for Azx, and black for Atot) for (α 1 , α 2 , α 3 ) cavity C 1 and (β 1 , β 2 , β 3 ) cavity C 2 , source location (S 1 ) in the three frequency bands.

	-1.8				(β 1 )
	-2.5
	-3
		(α 1 )
	-4	0.71 0.72 0.73 0.74	0.71 0.72 0.73 0.74
	-2.5
	-3				(β 2 )
	-4	(α 2 )
	-5
		0.81	0.82	0.83	0.805 0.81 0.815 0.82
	-5 -4.5				(β 3 )
	-5.5
		(α 3 )
	-6.5
		1.202 1.204 1.206	1.202 1.204 1.206 1.208
		f (GHz)		f (GHz)

  Frequency means of anisotropy coefficient amplitudes (in dB), source at S 1 .

			TABLE V
			Axy	Ayz	Azx	Atot
	700M Hz	C1 C2	-2.44 -3.23	-2.40 -3.44	-2.54 -3.26	-2.67 -3.31
	800M Hz	C1 C2	-3.05 -4.28	-2.96 -4.29	-2.90 -4.24	-2.97 -4.27
	1200M Hz	C1 C2	-4.88 -6.03	-4.94 -6.13	-4.81 -6.24	-4.88 -6.13

  Frequency means on the three frequency bands of the total standard deviation σtot and the total anisotropy coefficient Atot for 7 different source locations S i .

				TABLE VI
			700MHz	800MHz	1200MHz
			C1	C2	C1	C2	C1	C2
	S 2	σtot Atot	3.693 -2.323	2.209 -3.261	3.344 -3.065	1.837 -4.280	2.448 -4.693	1.651 -6.051
	S 3	σtot Atot	3.315 -2.525	2.251 -3.286	3.246 -3.047	1.983 -4.253	2.107 -4.828	1.662 -6.089
	S 4	σtot Atot	3.744 -2.459	2.255 -3.262	3.158 -3.071	1.881 -4.300	2.074 -4.872	1.675 -6.068
	S 5	σtot Atot	3.386 -2.601	2.218 -3.281	3.150 -3.168	1.901 -4.302	2.215 -4.852	1.562 -6.033
	S 6	σtot Atot	3.501 -2.496	2.344 -3.275	3.221 -3.102	1.994 -4.250	2.044 -4.793	1.611 -6.033
	S 7	σtot Atot	3.451 -2.410	2.269 -3.262	3.265 -3.019	1.889 -4.291	1.929 -4.834	1.566 -6.072
	S 8	σtot Atot	3.608 -2.462	2.396 -3.269	3.134 -3.067	1.986 -4.251	2.125 -4.804	1.604 -6.052

This work was supported by the French National Research Agency (ANR) under the project CAOREV. K. Selemani is with CMN Cherbourg, 51 rue de la Bretonnière, 50105 Cherbourg Cedex, France (kselemani@cmncherbourg.com). E. Richalot and O. Picon are with the Université Paris-Est, ESYCOM (EA 2552), UPEMLV, ESIEE-Paris, CNAM, F-77454 Marnela-Vallée, France (elodie.richalot, odile.picon)@univ-mlv.fr. O. Legrand and F. Mortessagne are with the LPMC, CNRS UMR 7336, Université de Nice-Sophia Antipolis, F-06108 Nice Cedex 2, France (olivier.legrand, fabrice.mortessagne)@unice.fr.

Fabrice MORTESSAGNE is currently a Physicist, a Professor at the University Nice Sophia Antipolis (UNS), and the Director of Laboratoire de Physique de la Matière Condenseé (UNS and CNRS). His research activities include field of wave propagation in complex media and mesoscopic physics. His current works include chaotic scattering experiments in optics and microwaves, and microwave realizations of artificial graphene and other relativistic meta-materials.


G De Filippis

J P No Ël

G Kerschen
email: g.kerschen@ulg.ac.be

L Soria

C Stephan

Experimental Nonlinear Identification of an Aircraft with Bolted Connections

Keywords: Nonlinear system identification, frequency-domain subspace method, full-scale aircraft, cubic splines, spurious poles

Aircraft structures are known to be prone to nonlinear phenomena, especially as they constantly become lighter and hence more flexible. One specific challenge that is regularly encountered is the modeling of the mounting interfaces between aircraft subcomponents. Indeed, for large amplitudes of vibration, such interfaces may loosen and, in turn, trigger complex mechanisms such as friction and clearances. In this context, the present work intends to investigate the nonlinear dynamics of the Morane-Saulnier Paris aircraft, accessible at ONERA. This aircraft possesses multiple bolted connections between two external fuel tanks and wing tips. The objective of the paper is specifically to carry out an adequate identification of the numerous nonlinearities affecting the dynamics of this full-scale structure. Nonlinearity detection and the subsequent subspace-based parameter estimation have been performed on experimental data, collected during an on-ground test campaign. Nonlinearity detection is first achieved by the comparison of frequency response functions estimated at low excitation level, with those obtained at high amplitude level, yielding insight towards accurately characterizing the behavior of the bolted connections. Then, a nonlinear subspace identification method is applied to measured data to estimate the linear and nonlinear parameters of the structure and novel strategies and tools that overcome specific arisen problems are developed.

INTRODUCTION

The increasing levels of performance demanded to aircraft structures can require the implementation of design procedures and options that may, in turn, lead to the activation of some nonlinear phenomena. Numerous works reporting on the evidence of existing nonlinearities in aerospace applications can be found in the technical literature, as, to give only some examples, Ref. [START_REF] Ahlquist | Assessment of nonlinear structural response in A400M GVT[END_REF][START_REF] Carrella | Identifying and quantifying structural nonlinearities in engineering applications from measured frequency response functions[END_REF][START_REF] No Ël | Nonlinear dynamic analysis of an F-16 aircraft using GVT data[END_REF] . Classical linear subspace-based methods [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation and Applications[END_REF] have recently been generalized to handle these occurrences [START_REF] Marchesiello | A time domain approach for identifying nonlinear vibrating structures by subspace methods[END_REF][START_REF] No Ël | Frequency-domain subspace identification for nonlinear mechanical systems[END_REF] . Specifically, the frequency-domain nonlinear subspace identification (FNSI) method has been proved to successfully deal with grey-box models, affected by many nonlinearities [START_REF] No Ël | Grey-box identification of a non-linear solar array structure using cubic splines[END_REF] .

In this paper, the FNSI methodology is applied for the first time to a full-scale structure, the Morane-Saulnier Paris aircraft. This aircraft presents multiple bolted connections between the wings and two external fuel tanks, located at the wing tips. The experimental identification requires that a considerable number of output locations are utilized on the aircraft main body and at the two the ends of the mounting interfaces, where nonlinear behaviors are expected to arise. Due to the presence of many nonlinearities, a relevant number of basis functions are included in the FNSI routine. The development of novel tools and strategies is carried out to overcome the problems related to the high dimensionality of the consequent inverse problem. Moreover, classical subspace methods do not provide the user with the possibility of selecting among the poles computed at a given model order. This means that, in the usual procedure, all the calculated poles have to be included, the physical as the spurious ones. To avoid this drawback, an appropriate strategy is employed, that allows for removing the undesired spurious poles.

Linear identification is first performed at low amplitude level, with the aim of collecting the modal parameters. Then, the comparison between frequency response functions (FRFs) estimated at low and at high amplitude level is employed for detecting the overall nonlinear behavior of the structure. Afterward, the FNSI method is applied to estimate the underlying linear properties and the nonlinear restoring forces. Finally, a novel strategy for spurious poles' removal is applied to the identified model.

DESCRIPTION OF THE AIRCRAFT AND OF THE EXPERIMENTAL SETUP

The structure is the airframe of the Morane-Saulnier Paris aircraft, represented in Figure 1. This French jet aircraft was built during the 1950s and was used as a trainer and liaison aircraft. The structural configuration under consideration corresponds to the aircraft without its jet engines and standing on the ground through its three landing gears with deflated tires. A specimen of this airplane is present in ONERA's laboratory. It is used for training engineers and technicians for ground vibration tests (GVT) and for research purposes. GVT exhibit nonlinear behavior in the connection between the wings and external fuel tanks located at the wing tip. A total of 38 accelerometers is employed to instrument the structure, of which 8 tri-axial sensors are positioned at the two ends of each of the 4 wing-tank bolted attachments (Figure 2), where, thus, nonlinearities are mainly expected to be located. Two small-size shakers, one visible in Figure 1, are used to apply pure broad-band random excitation inputs to the two wings, along the vertical direction. The adopted sampling frequency is equal to 400 Hz. Acceleration and force signals are both recorded at the driving points, located below the fuel tanks (Figure 3), by means of impedance heads. In this paper, all the output channels related to bolted attachments are denoted by a label of the form abc:D, where a is Ri or Le, for indicating the right or the left side of the aircraft, respectively, b is Fr or Re, when referring to the front or the rear bolted connection, and c is Wi or Ta, if the considered node belongs to the wing or to the tank, respectively. More over, D represents the channel direction (X, Y, or Z), evaluated with respect to a fixed reference frame, having the X-axis coinciding with the longitudinal axis of the aircraft, but directed from the nose to the tail, the vertical, Z-axis, pointing upward, and the Y-axis chosen to form a right-handed cartesian coordinate system. Both the ends of the bolted connections are included when the character c is imposed equal to BC. To give an example of this nomenclature, the label LeReWi:Y denotes the output channel referred to the left rear bolted connection, node on the wing, along the Y direction. In the case of input channels, instead, we need only to distinguish between the left (LeIn:Z) and the right (RiIn:Z) excitation.

NONLINEAR SUBSPACE IDENTIFICATION IN THE FREQUENCY DOMAIN

The FNSI method derives models of mechanical systems possessing localised nonlinearities directly from measured data, and without resorting to a preexisting numerical model, e.g. a finite element model [START_REF] No Ël | Frequency-domain subspace identification for nonlinear mechanical systems[END_REF] . It is applicable to multi-input, multi-output structures with high nonproportional damping and high modal density, and makes no assumption as to the importance of nonlinearity in the measured dynamics [START_REF] No Ël | Subspace-based identification of a nonlinear spacecraft in the time and frequency domains[END_REF] . The vibrations of such nonlinear systems are governed by Newton's second law M q(t) + C q(t) + K q(t) + g(q(t), q(t)) = p(t)

where M, C, K ∈ R n×n are the mass, linear viscous damping and linear stiffness matrices, respectively; q(t) and p(t) ∈ R n are the generalized displacement and external force vectors, respectively; g(t) ∈ R n is the essentially nonlinear, i.e. nonlinearizable, restoring force vector encompassing elastic and dissipative contributions, and n is the number of degrees of freedom (DOFs) of the structure obtained after spatial discretization. The amplitude, direction, location and frequency content of the excitation p(t) determine in which regime, linear or nonlinear, the structure behaves. The effects of the r lumped nonlinear components in the system are represented using a linear-in-the-parameters model of the form

g(q(t), q(t)) = r ∑ a=1 sa ∑ b=1 c a,b h a,b (q(t), q(t)). (2)
In this double sum, sa is the number of nonlinear basis functions h a,b (q(t), q(t)) selected to describe the a-th nonlinearity, and c a,b are the associated coefficients. The total number of nonlinear basis functions introduced in the model is equal to s = ∑ r a=1 sa. Linearity in the parameters avoids an iterative optimization process, and issues related to initialization and convergence thereof.

Given measurements of p(t) and q(t) or its derivatives, and an appropriate selection of the functionals h a,b (t), the FNSI algorithm aims at computing estimates of (i) the FRF matrix

H(ω) = ( -ω 2 M + j ω C + K ) -1 ( 3
)
describing the underlying linear properties of the system in Eq. 1, where ω is the frequency and j the imaginary unit, and (ii) the nonlinear coefficients c a,b in Eq. 2.

The identification methodology essentially builds on a block-oriented interpretation of nonlinear structural dynamics, which sees nonlinearities as a feedback into the linear system in the forward loop [START_REF] Adams | A frequency domain method for estimating the parameters of a non-linear structural dynamic model through feedback[END_REF] , as illustrated in Fig. 4. This interpretation boils down to moving the nonlinear internal forces in Eq. 1 to the right-hand side, and viewing them as additional external forces applied to the underlying linear structure, that is

M q(t) + C q(t) + K q(t) = p(t) - r ∑ a=1 sa ∑ b=1 c a,b h a,b (q(t), q(t)). (4)
Based on Equation 4, the parameter estimation problem is solved in the frequency domain by processing a user-selected number Underlying linear system: M, C, K

Nonlinear feedback: c a,b , h a,b (q(t), q(t)) + p(t) q(t), q(t), q(t)

Figure 4: Feedback interpretation of nonlinear structural dynamics [START_REF] Adams | A frequency domain method for estimating the parameters of a non-linear structural dynamic model through feedback[END_REF] . of frequency lines in the measured band. A state-space representation of Eq. 4 is first extracted directly from frequency-domain data using robust tools from numerical analysis, including QR and singular value decompositions. The calculated state-space parameters are secondly converted into the modal space to estimate the underlying linear FRF matrix H(ω), and to the physical space to estimate the nonlinear coefficients c a,b . Fig. 5 presents an overview of the identification methodology. The interested reader is referred to Ref. [START_REF] No Ël | Frequency-domain subspace identification for nonlinear mechanical systems[END_REF] for a detailed introduction to the theoretical and practical aspects of the FNSI method.

LINEAR ANALYSIS AT LOW LEVEL AND NONLINEARITY DETECTION

When the system is excited with low amplitude level inputs, the nonlinearities affecting the bolted connections are not activated. This means that data acquired at low level can be utilized to perform a linear identification of the structure. To this aim, a band-1. Select appropriate basis functions h a,b (q(t), q(t)) to represent the nonlinearities.

2. Choose the number of processed frequency lines in the measured band.

3. Determine the order of the state-space representation of Eq. 4 and calculate its parameters from frequency-domain data.

4. Convert the state-space model to modal space to estimate the underlying linear FRF matrix H(ω).

5. Convert the state-space model to physical space to estimate the nonlinear coefficients c a,b .

Figure 5: Overview of the FNSI methodology. limited (7 -45 Hz) white-noise signal is applied to the aircraft, having a root-mean-square (RMS) amplitude of 4.4 N. In this conditions, the FNSI method can be employed for obtaining the linear modal parameters of the structure. Since we, here, focus only on the identification of the first 5 modes, the frequency band of analysis ranging from 7 to 21.5 Hz is considered. A detailed motivation of this choice is reported in subsection 5.2. Figure 6 shows the linear properties extracted by the FNSI method for model orders up to 50. This stabilization diagram leads to the identification of 5 mode shapes, whose natural frequencies and damping ratios are listed in Table 1. The obtained results are therein compared with those estimated by using another state-of-the-art algorithm, i.e. the poly-reference Least-Squares Complex Frequency Domain (p-LSCFD [START_REF] Peeters | The PolyMAX frequency-domain method: A new standard for modal parameter estimation?[END_REF] ). Furthermore, in Figure 7 the comparison of results is performed by means of the modal assurance criterion (MAC [START_REF] Heylen | Modal Analysis Theory and Testing[END_REF] ). By considering what is shown we may conclude that the two methods give basically the same results, both in terms of natural frequencies and mode shapes, while a bias error inevitably affects the estimate of the damping ratios. However, a good agreement is found when comparing FRFs synthesized by the FNSI method and those obtained with the H1 estimator (Figure 8).

A second random data set measured at high level (44 N RMS) is then analyzed to detect and locate in the frequency domain the modes triggering nonlinear mechanisms. Since all the considered mode shapes involve to a certain extent the wings' deformation, all the frequencies seem to be affected by nonlinear distortions. To give an example, the deformed shapes of modes 2 and 4 are depicted in Figures 9a and9b, respectively. In Figure 10, the comparison between the low-and high-level FRFs is presented. A softening behavior is clearly visible for all the first 3 modes, included in the range from 9 to 13 Hz. Modes 4 and 5 represent the most critical configurations in terms of deformed shapes. During the resonance, in fact, the bolted attachments are activated through a combination of wing bending and torsion. As a consequence, the related modal parameters are strongly affected by nonlinear distortions. Since the peak located at 18.7 Hz is more depressed at high level, the interpretation of the nonlinear behavior becomes arduous. However, by looking at the phase plot, one may conclude that the noticeable curve translation is still induced by softening effect. We instead demonstrate in what follows, that this phenomenon can likely to be attributed to loosening in the bolted connections. In fact, even if at high level an overall softening effect has been detected, we later show that some connections can exhibit a combined softening-hardening behavior.

p-LSCFD Discrepancy Mode fn (Hz) ζ (%) fn (Hz) ζ (%) δ f (%) δ ζ (%) 1

NONLINEARITY IDENTIFICATION AT HIGH LEVEL

To obtain a nonlinear subspace model of the aircraft at the high 44 N RMS level, a five-step procedure is followed (see Figure 5). The first step is an appropriate selection of the nonlinear basis functions h a,b (q (t) , q (t)) defined in Eq. 2. The analysis reported in Section 4 has revealed that the physical phenomena activating the nonlinear behavior of the structure are complex and numerous. Since the aircraft is equipped with 4 bolted attachments and the acceleration is measured at the two ends along each direction, 24 output channels are candidate to form the right basis functions. In the following section, two useful tools are introduced: i) cubic splines, that avoid the selection of physics-based functional forms; ii) the comparison approach, that allows for evaluating the influence of the selected basis functions on the retrieved underlying linear properties.

Selection of appropriate basis functions

For the sake of simplicity and without loss of generality, the scope of this section is restricted to an elastic nonlinear restoring force g(q(t)) in the single-DOF case. Let q be divided into L segments of arbitrary length and defined by their abscissas, denoted by q k for k = 1, . . . , L + 1. Each abscissa is associated with an ordinate g k , together defining a knot (q k ,g k ) of the spline. Thus, if q is a displacement value in between knots k and k + 1, the corresponding point of a cubic-spline-based approximation of g(q(t)) is given by [START_REF] De Boor | A Practical Guide to Splines[END_REF] g(q(t)) = (2t 3 -3t 2 + 1)

g k + (-2t 3 + 3t 2 ) g k+1 + (t 3 -2t 2 + t)(q k+1 -q k ) g ′ k + (t 3 -t 2 )(q k+1 -q k ) g ′ k+1 ( 5
)
where t is the normalised displacement q-q k q k+1 -q k . The computation of the first derivatives g ′ k = ∂g k /∂q k can be achieved by forcing the cubic spline and its first two derivatives to be continuous across each of the interior knots. This results in L -1 linear constraint equations

g ′ k-1 q k -q k-1 + 2 ( 1 q k -q k-1 + 1 q k+1 -q k ) g ′ k + g ′ k+1 q k+1 -q k = 3 ( g k -g k-1 (q k -q k-1 ) 2 + g k+1 -g k (q k+1 -q k ) 2 ) . ( 6
)
Since the essentially nonlinear restoring force g(q(t)) is zero and has zero slope at equilibrium, one should also enforce, in the segment containing the abscissa of the equilibrium point, that

(t 3 0 -2t 2 0 + t0)(q k+1 -q k ) g ′ k + (t 3 0 -t 2 0 )(q k+1 -q k ) g ′ k+1 = -(2t 3 0 -3t 2 0 + 1) g k -(-2t 3 0 + 3t 2 0 ) g k+1 (7)
and

(3t 2 0 -4t0 + 1)(q k+1 -q k ) g ′ k + (3t 2 0 -2t0)(q k+1 -q k ) g ′ k+1 = 6(t0 -t 2 0 )(g k -g k+1 ) (8)
where t0 = -q k q k+1 -q k . Eqs. 6, 7 and 8 constitute L + 1 relations that uniquely define the L + 1 parameters g ′ 1 , . . . , g ′ L+1 as functions of the ordinates of the knots g1, . . . , gL+1. The first derivatives can thus be substituted in Equation 5to compute the basis functions associated with the ordinates. These basis functions correspond to the terms h a,b (t) introduced as additional external forces in the FNSI algorithm, as described in Section 3.

The optimal number of knots in the definition of the splines should rigorously be sought by minimizing the difference, in some metric, between the predictions of the nonlinear model and measured data. In practice, an acceptable number of knots can be found by a trial and error approach, by maximizing the accuracy of the underlying linear properties of the system, estimated from nonlinear data. However, when a large number of outputs is employed to form the basis functions, the inverse problem becomes difficult to be solved. This is mainly due to the high dimensionality of data matrices involved in the computation. This means that the maximum allowed number of basis functions is inevitably influenced by the size of the data set utilized in the identification. In this framework, the comparison approach represents a novel strategy, developed to solve the above described problems. The main idea is to single out the minimum number of outputs and the simplest functional forms, which are able to retrieve at high level enough accurate underlying linear properties. The first step is to collect natural frequencies, damping ratios and real mode shapes estimated by the FNSI method at low level. Tolerance thresholds are imposed by the analyst to these three parameters. In particular, the tolerance on natural frequencies is utilized to form the acceptance regions. Of course, the number of these regions has to be equal to the number of the collected natural frequencies. Then, at high level, the outputs that are candidates to form the basis functions have to be combined in a series of scheduled experiments. The underlying linear properties are extracted for all the combinations and for each model order. The obtained solutions are stacked in the acceptance regions and the discarded linear properties are immediately classified as non comparable in natural frequency. At each model order, at least one admitted solution must be found in all the acceptance regions. An ultimate step is needed to reveal the so called comparable model orders, in which the underlying linear properties are retrieved within certain tolerance limits. A graphical tool, that we name comparison diagram, is developed to make simple and effective the visualization of the comparable model orders and of equivalent solutions. On the x-axis, the natural frequencies are represented, while the model orders are given on the y-axis; horizontal green lines indicate the comparable model orders. As in the case of stabilization diagram, several situations are labeled with different symbols: The dot denotes no equivalence in natural frequency, the cross represents equivalence in natural frequency, the square and the circle stand for extra equivalence in damping ratio and in MAC value, the triangle refers to fully equivalent solutions.

For the investigated aircraft, the modal parameters estimated at low level are collected in Table 1. Comparison thresholds for natural frequency, damping ratio and MAC value are set to 1.25 %, 10 % and 0.95, respectively. The motivation of these choices is later discussed in Subsection 5.4. We choose to include in each experiment, only one direction for each bolted connection. Thus, employing a full factorial Design of Experiments with 4 factors (the 4 attachments) and 3 levels (the 3 directions), a table of 81 possible combinations is obtained. With this strategy, only 8 of the 24 output channels related to the bolted connections have to be included in the FNSI routine. By the comparison approach, it is found that splines with 5 knots on relative displacements are the simplest functional forms able to suitably capture the physical mechanisms activating the nonlinear phenomena. It is worth to here stress the reduced dimensionality of the consequent inverse problem, that encompasses 20 basis functions, 2 external forces, 32 outputs for the aircraft main body and 8 outputs for the bolted connections.

Choice of the processed bandwidth

As anticipated in Section 4, the bandwidth processed for the identification ranges from 7 to 21.5 Hz. This choice is principally made for reducing the computational burden. In fact, the dimensionality reduction operated through the comparison approach may become not sufficient when the FNSI method has to deal with a considerable number of basis functions. The restriction of the bandwidth is helpful for avoiding memory issues, even if does not allow for preserving the nonlinear distortions usually appearing outside the whole band of excitation (from 7 to 45 Hz). Owing to this dimensionality, time-domain identification techniques, as, in particular, the subspace method proposed in Ref. [START_REF] Marchesiello | A time domain approach for identifying nonlinear vibrating structures by subspace methods[END_REF] , are likely to suffer from computational memory issues, because of their inability to reject less informative measured samples. Moreover, few methods in the technical literature are capable of reliably handling a large amount of nonlinear coefficients. For instance,the conditioned reverse path method Ref. [START_REF] Richards | Identification of multi-degree-of-freedom non-linear systems under random excitations by the reverse path spectral method[END_REF] estimates the coefficients sequentially, i.e. errors are accumulated throughout the identification process, and would probably lead to very large errors in the present case study. Another example is the nonlinear identification through feedback of the outputs method Ref. [START_REF] Adams | A frequency domain method for estimating the parameters of a non-linear structural dynamic model through feedback[END_REF] , whose estimation accuracy is known to decrease rapidly as the number of parameters increases.

Determination of the model order

The third step in Figure 5 is the determination of the model order. One can notice that, here, the information about the poles' stabilization may become ambiguous and difficult to be interpreted (Fig. 11). In modal analysis, the stabilization diagram is often used to assist the analyst in separating the physical system poles from the mathematical ones. Since the linear identification has been accomplished at low level, the correct physical poles have been already collected. This means that in the comparison diagram the spurious poles are a priori discarded, since they are classified as non comparable in natural frequency at the beginning of the analysis. On the contrary, the stabilization diagram may lead to the selection of stable, but not comparable underlying linear properties. These considerations remark the importance of the proposed strategy, in which the comparison approach represents a useful methodology that allows for: i) refining the localization of nonlinearities, and ii) revealing the presence of comparable model orders.

Estimation of the underlying linear properties

In Figure 12 we show the comparison diagram corresponding to the experiment that involves the outputs denoted as RiReBC:Z, RiFrBC:Z, LeReBC:Y and LeFrBC:X. Three consecutive model orders, i.e. the 42, 44 and 46, are found to be comparable. An extra equivalence in the MAC value is evident only for the first mode. Depending on the chosen model order, the damping ratio may result comparable on the third or on the second mode. Of course, these outcomes are absolutely dependent on the imposed tolerance limits. We stress that the comparison approach is really sensitive to tolerance fixed on the natural frequencies, mainly owing to the construction of the acceptance regions. Obviously, the larger are these regions, the higher is the chance of finding comparable model orders. On the opposite, using strict tolerance values leads to the estimate of more accurate underlying linear properties. Moving from these considerations, the analyst is suggested to set the comparison thresholds according to the application. In Figure 10 a small shift can be noticed for the first 3 peaks of the FRF estimated at high level. This means that, for the present investigation, the adoption of severe tolerance limits is strongly recommended.

The FRF LeReWi:Y/RiIn:Z synthesized by the FNSI at model order 42 is depicted in Figure 13. When compared to the synthesis obtained at low level, many additional peaks are now visible, which corresponds to the numerous spurious poles present at this model order. Due to the extra equivalence on the damping ratio, a good agreement is observable only at the third resonant frequency (12.75 Hz). A useful way to improve this result is presented in Section 6.

Estimation of the nonlinear coefficients

With the last step in Figure 5, the state-space model is converted to the physical space, in order to visualize the restoring forces, fitted by using cubic splines. To this end, the conversion scheme detailed in Ref. [START_REF] No Ël | Frequency-domain subspace identification for nonlinear mechanical systems[END_REF] is utilized. As a result, the nonlinear coefficients identified from the extended FRF matrix He (ω) are spectral quantities, i.e. are complex-valued and frequencydependent. A reliable identification, together with an appropriate selection of the basis functions, should make the imaginary parts much smaller than the corresponding real parts. The coefficients' frequency dependence should also remain small. These indications are usually utilized as criteria for the quality assessment of the estimated nonlinear parameters. In Figure 14, we illustrate the real and the imaginary parts of the complex-valued coefficients related to the output LeFrBC:X, whose estimation results to be that more affected by errors. The variations along the frequency axis are located at the same positions of the selected modal parameters. As one can clearly observe, on both the curves the spurious poles give rise to sharpest peaks, in turn, heavily affecting the average values of the estimated coefficients. To mitigate this effect, the frequency range for the conversion is restricted to the flat zone ranging from 13 to 14 Hz (Fig. 14). In this range, the ratios between the averaged real and imaginary parts are computed, actually revealing the output LeFrBC:X as the worst estimated among the four of the experiment (see the fourth row of Table 2). This result is even improved in the next section (Sec. 6), by using an appropriate strategy specifically designed for removing the spurious poles.

PROPOSED STRATEGY FOR REMOVING SPURIOUS POLES AND IMPROVED RESULTS

The presence of the spurious poles can be attributed to many different causes, e.g. noise on the data, modeling errors, computational issues, etc. Classical subspace methods do not provide the user with the possibility of selecting poles computed at a given model order, not allowing for discarding the spurious ones. In the usual procedure, in fact, all the poles have to be included.

To avoid this drawback, an appropriate novel strategy is employed for removing the spurious poles. In the first part of the FNSI method, the underlying linear properties are directly computed from the continuous state and output matrices Ac ∈ R 2n×2n and C ∈ R l×2n . Then, the discrete extended input and direct feedthrough matrices B e d ∈ R 2n×(m+s) and D e ∈ R l×(m+s) are calculated by defining the transfer function matrix Gs (k) associated with the state-space model

Gs (k) = C ( z k I n×n -A d ) -1 B e d + D e (9)
and minimizing the weighted difference between the measured and modeled output spectra in a linear least-squares sense

B e d , D e = arg min F ∑ k=1 W 2 (k)|Y (k) -He (k) E (k)| 2 (10)
where z k = e j2πk/N is the Z-transform variable for discrete-time models, W(k) is a real-valued weighting function, Y(k) and E(k) represent the output spectra and the extended input vector, F denotes the number of frequency samples, the subscripts c and d refer to continuous and discrete time models, while l and m indicate the number of outputs and the number of inputs involved in the identification, respectively. The continuous-time state-space matrices can be converted in discrete-time form by means of the following relations, considering a sampling frequency fs

A d = expm (Ac/fs) B e d = ( expm (Ac/fs) -I 2n×2n ) A -1 c B e c ( 11
)
where expm indicates the matrix exponential. The employed procedure for the removal of spurious poles consists in the reduction of the dimensionality of the state-space matrices. From the eigenvalue decomposition of Ac the poles of the system and the eigenvectors are collected in Λc ∈ C 2n×2n and V ∈ C 2n×2n , respectively. Once np poles are selected by the analyst, the corresponding diagonal terms in Λc are extracted and a reduced matrix of system poles is obtained Λc ∈ C 2nr ×2nr (nr = 2 • np).

The same procedure of extraction is repeated on the columns of V, leading to V ∈ C 2n×2nr . The reduced output matrix C ∈ C l×2nr is, thus, calculated in modal coordinates

C = C V (12)
and the following Eq. 13 is utilized to replace Eq. 9 in the minimization problem, in order to derive the Be d and D e matrices

Gs (k) = C( z k I n×n -Λd ) -1 Be d + D e . (13)
The array Be d ∈ C 2nr ×(m+s) is the discrete extended input matrix represented in modal coordinates. We comment that for converting Λc and Be d to discrete and continuous time forms the same relations holding for Ac and B e d (Eq. 11) can be employed, respectively. Finally, the extended FRF matrix He (ω) is estimated and the typical scheme for conversion to the physical space can be applied to retrieve the nonlinear coefficients. In the comparison diagram of Figure 12 the blue marks indicate the selected system poles at model order 42. In Table 3, the linear properties identified at low level (see Section 4), together with the associated frequencies estimated at high level are reported. As one can see, the error values are respectful of the imposed comparison thresholds and this means that, since strict tolerances can be fixed on natural frequencies, a consequent matching is found between estimates at low and high levels.

FNSI at Low

Differently, the overall estimate of damping ratios is less satisfactory. The damping ratios that present higher discrepancies are those corresponding to modes 4 and 5. Probably, this fact may be related to the activation of nonlinear phenomena involving damping forces. This behavior, already singled out during the detection step, is further stressed and put in evidence by the comparison diagram. In fact, the last region ranging from 17 to 21 Hz is very confused and really affected by the presence of not comparable solutions. We even comment that the splines utilized in the present identification are not able to capture this kind of behavior, being solely defined on the relative displacement. Future investigations will be made in this direction. The results presented in Table 3 have a direct influence on the synthesis of FRFs. The same FRF described in the previous section (Sec. 5), i.e. the one of output LeReWi:Y/RiIn:Z, is depicted in Figure 15, where 5 distinct peaks are now clearly visible also at high level. Due to a more reliable estimate of the damping ratio, a good agreement is confirmed for the third mode, and a relevant improvement is also achieved in the cases of the first two frequencies, as evidenced by the phase plot. Again by inspecting the comparison between FRF syntheses, we may conclude that mode 4 and 5 are the most distorted and, as expected, a huge discrepancy is observed in the magnitude plot.

Location

Coeff Again, the last step is the conversion of the reduced state-space model to the physical space. In Figure 16 the complex-valued coefficients estimated by the FNSI method in the whole frequency bandwidth are represented. With respect to Figure 14 the diagram appears highly improved. The sharpness of the peaks is now mitigated both for the real and the imaginary parts.

Extremely interesting results are obtained by performing the estimate in the frequency range from 15 to 16.5 Hz (see Table 4), where the ratio between the averaged real and imaginary parts, in log scale, are now positive for all the coefficients, meaning that, finally, a good estimate has been accomplished. We comment that the trend of the estimated nonlinear restoring forces is consistent with the nature of the bolted connections (see Fig. 17). In fact, a combined softening-hardening behavior can be observed for RiReBC:Z, RiFrBC:Z and LeReBC:Y, while in the case of output LeFrBC:X, a simple softening effect is found out.

It is worth to stress that due the high stiffness, the nonlinear behavior of bolted connections is activated in the very small scale of relative displacements (10 -2 mm). Moreover, a high density of points is present inside the softening region (±3 • 10 -2 mm), that may probably explain the overall softening behavior observed in the detection step (Sec. 4), and that can now certainly considered imputable to loosening in the attachments.

CONCLUSIONS

In the presented research, new tools and strategies have been implemented in the FNSI methodology to tackle the experimental identification of a full-scale aircraft, affected by the presence of many nonlinearities. The strength of the proposed comparison approach consists in the possibility of reducing the dimensionality of the inverse problem, avoiding the memory issues that may arise when large data sets are employed. On the other hand, the nonlinear identification can be carried out without loosing relevant information that might be related to the not included outputs. In addition, this strategy can be utilized to evaluate the capability of the chosen basis functions of encompassing the activated nonlinear phenomena. In this context, the limits of the stabilization analysis have been exploited at high level. A new visualization tool, named comparison diagram, is used to reveal the existence of comparable model orders in which the underlying linear properties are recovered within certain tolerance limits. At high level, cubic splines with 5 knots on the relative displacement, related to the selected outputs, are included in the FNSI routine. An improved agreement is found when comparing FRFs synthesized at low and high levels, after a novel procedure for spurious poles' removal is applied. In the end, the estimated nonlinear restoring forces have proved to be consistent with the nature of the bolted connections. Despite the fact that an overall softening behavior has been observed in the detection step, the estimates have revealed an additional combined softening-hardening behavior, that is certainly imputable to loosening in the attachments. A less satisfying agreement is found out on estimates of damping ratios at high level, that may be probably related to the activation of nonlinear phenomena involving the damping forces. In the present study splines on velocity were not considered, but this will be the focus of our future investigations. More over, to refine the achieved results effort will be spent for improving the proposed novel tools and procedures.

Figure 1 :

 1 Figure 1: Morane-Saulnier Paris aircraft.

Figure 2 :

 2 Figure 2: Connection between external fuel tank and wing tip (top view). Close-up of (a) front and (b) rear bolted attachments.

Figure 3 :

 3 Figure 3: Layout of the experimental setup (top view). Black circles: output locations; red circles: driving points.

Figure 6 :

 6 Figure 6: Stabilization diagram computed by the FNSI method at low level (4.4 N RMS). Dot: no stabilization in natural frequency; cross: stabilization in natural frequency; square: extra stabilization in damping ratio; circle: extra stabilization in MAC value; triangle: full stabilization. Stabilization thresholds for natural frequency, damping ratio and MAC value are 2 %, 10 % and 0.95, respectively. The green line indicates the selected model order.

Figure 7 :

 7 Figure 7: Comparison between FNSI and p-LSCFD results. (a) AutoMAC for FNSI method; (b) MAC between FNSI and p-LSCFD mode sets.

Figure 8 :

 8 Figure 8: FRF LeReWi:Y/RiIn:Z. Comparison between H1 estimate and FNSI synthesis at low level (4.4 N RMS).

FNSI

Figure 9 :

 9 Figure 9: Deformed shapes corresponding to modes 2 and 4 identified at low level (4.4 N RMS). Geometrical dimensions are given in meters. (a) Mode 2 (10.64 Hz) and (b) mode 4 (18.71 Hz). Blue and red lines indicate the undeformed and deformed shapes, respectively.

Figure 10 :

 10 Figure 10: FRF LeReWi:Y/RiIn:Z. Comparison between H1 estimate at low (4.4 N RMS) and high (44 N RMS) levels.

Figure 11 :

 11 Figure 11: Stabilization diagram computed by the FNSI method at high level (44 N RMS). Dot, crosses, squares, circles and triangles are defined similarly to Fig. 6.

Figure 12 :

 12 Figure 12: Comparison diagram computed at high level (44 N RMS). Splines with 5 knots on relative displacements of the outputs denoted as RiReBC:Z, RiFrBC:Z, LeReBC:Y and LeFrBC:X are included in the FNSI routine. Dot: no equivalence in natural frequency; cross: equivalence in natural frequency; square: extra equivalence in damping ratio; circle: extra equivalence in MAC value; triangle: full equivalence. Comparison thresholds for natural frequency, damping ratio and MAC value are 1.25 %, 10 % and 0.95, respectively. The green lines are the comparable model orders and the blue marks indicate the selected system poles.

Figure 13 :

 13 Figure 13: FRF LeReWi:Y/RiIn:Z. Comparison between FNSI synthesis at low (4.4 N RMS) and high (44 N RMS) levels.

Figure 14 :

 14 Figure 14: Complex-valued estimates of the 5 nonlinear coefficients related to the output LeFrBC:X in the full bandwidth (7 -21.5 Hz).

Figure 15 :

 15 Figure 15: Improved FRF LeReWi:Y/RiIn:Z. Comparison between FNSI synthesis at low (4.4 N RMS) and high (44 N RMS) levels. The selected poles are depicted in Fig. 12.

Figure 16 :

 16 Figure16: Improved complex-valued estimates of the 5 nonlinear coefficients related to the output LeFrBC:X in the full bandwidth (7 -21.5 Hz). The selected poles are depicted in Fig.12..

Figure 17 :

 17 Figure 17: Improved nonlinear stiffness curves reconstructed at high level (44 N RMS) at different output locations. Red circles localise the 5 chosen knots of the splines. The selected poles are depicted in Fig. 12.

TABLE 1 :

 1 Comparison

		9.02	0.29	9.02	0.32	0.00	9.38
	2	10.64	0.78	10.64	0.78	0.00	0.00
	3	12.75	0.67	12.70	0.56	0.39	19.64
	4	18.71	1.34	18.64	2.10	0.38	36.19
	5	20.23	0.86	20.11	0.84	0.60	2.38

between FNSI and p-LSCFD results. Estimated natural frequencies (fn) and damping ratios (ζ) in the 7 -21.5 Hz band.

TABLE 2 :

 2 Ratio between averaged real and imaginary parts in logarithmic scaling. The estimate is performed in the 13 -14 Hz band.

TABLE 3 :

 3 Comparison between FNSI estimates at low and high levels. Estimated natural frequencies (fn) and damping ratios (ζ) in the 7 -21.5 Hz band. The fourth and the fifth columns refers to underlying linear properties selected in Fig.12.

			Level FNSI at High Level	Error
	Mode fn (Hz)	ζ (%)	fn (Hz)	ζ (%)	ϵ f (%) ϵ ζ (%)
	1	9.02	0.29	8.94	0.42	0.89	44.83
	2	10.64	0.78	10.54	0.50	0.94	35.90
	3	12.75	0.67	12.62	0.72	1.02	7.46
	4	18.71	1.34	18.78	0.16	0.37	88.06
	5	20.23	0.86	20.00	0.20	1.14	76.74

TABLE 4 :

 4 Improved ratio between averaged real and imaginary parts in logarithmic scaling. The estimate is performed in the 15 -16.5 Hz band. The selected poles are depicted in Fig.12.

		. 1 Coeff. 2 Coeff. 3 Coeff. 4 Coeff. 5
	RiReBC:Z	0.96	1.72	1.94	1.68	0.95
	RiFrBC:Z	0.98	1.68	1.25	1.78	1.05
	LeReBC:Y	0.52	0.85	0.57	1.06	0.08
	LeFrBC:X	0.36	1.00	0.40	1.64	0.97

ACKNOWLEDGMENTS

The author J.P. No ël is a Postdoctoral Researcher of the Fonds de la Recherche Scientifique -FNRS which is gratefully acknowledged.


Rabaa Youssef

Anis Kacem

Sylvie Sevestre-Ghalila

Christine Chappard

Graph structuring of skeleton object for its high-level exploitation

Skeletonization is a morphological operation that summarizes an object by its median lines while preserving the initial image topology. It provides features used in biometric for the matching process, as well as medical imaging for quantification of the bone microarchitecture. We develop a solution for the extraction of structural and morphometric features useful in biometric, character recognition and medical imaging. It aims at storing object descriptors in a re-usable and hierarchical format. We propose graph data structures to identify skeleton nodes and branches, link them and store their corresponding features. This graph structure allows us to generate CSV files for high level analysis and to propose a pruning method that removes spurious branches regarding their length and mean gray level. We illustrate manipulations of the skeleton graph structure on medical image dedicated to bone microarchitecture characterization.

Introduction

Skeletonization is used in various applications such as biometrics [START_REF] Zhu | Hand Dorsal Vein Recognition Based on Hierarchically Structured Texture and Geometry Features[END_REF][START_REF] Zhao | Preprocessing and postprocessing for skeleton-based fingerprint minutiae extraction[END_REF], medical imaging [START_REF] Chappard | Subchondral bone micro-architectural alterations in osteoarthritis: a synchrotron micro-computed tomography study[END_REF][START_REF] Sevestre-Ghalila | Texture Image Analysis for Osteoporosis Detection with Morphological Tools[END_REF][START_REF] Tschirren | Shock Segmentation, Skeletonization, and Branchpoint Matching A Fully Automated Quantitative Evaluation of Human Intrathoracic Airway Trees[END_REF] and character recognition [START_REF] Wshah | Segmentation of Arabic Handwriting Based on both Contour and Skeleton Segmentation[END_REF] since it provides features that enables user to access high-level analysis of the image objects. In fact, object matching methods based on skeleton features are used in biometric identification through minutiae comparison of hand vein [START_REF] Zhu | Hand Dorsal Vein Recognition Based on Hierarchically Structured Texture and Geometry Features[END_REF] or digital fingerprint [START_REF] Zhao | Preprocessing and postprocessing for skeleton-based fingerprint minutiae extraction[END_REF], in bronchial airway trees monitoring [START_REF] Tschirren | Shock Segmentation, Skeletonization, and Branchpoint Matching A Fully Automated Quantitative Evaluation of Human Intrathoracic Airway Trees[END_REF], in symbols identification [START_REF] Siddiqi | Shock graphs and shape matching[END_REF] and in character recognition [START_REF] Wshah | Segmentation of Arabic Handwriting Based on both Contour and Skeleton Segmentation[END_REF]. In addition, classification methods using morphometric features extracted from skeleton helps diagnose osteoporosis [START_REF] Sevestre-Ghalila | Texture Image Analysis for Osteoporosis Detection with Morphological Tools[END_REF] and osteoarthritis [START_REF] Chappard | Subchondral bone micro-architectural alterations in osteoarthritis: a synchrotron micro-computed tomography study[END_REF] diseases. Graph-based representation of the skeleton is widely investigated for matching issues since the correspondence between skeleton branches with graph edges and nodes with its vertices is natural and intuitive. The proposed solutions in the literature are based on shock graphs [START_REF] Siddiqi | Shock graphs and shape matching[END_REF] or attributed relational graphs [START_REF] Di Ruberto | Recognition of shapes by attributed skeletal graphs[END_REF]. These representations based on adjacency matrix can only consider acyclic graphs which restricts the application domains. Our contribution is to develop a solution for the exploitation of skeleton features in binary and grayscale domain and to save them in a re-usable format easily manipulated. We propose to construct a graph able to represent any skeleton and to store its structural and morphometric features. To this purpose, we identify skeleton with multigraph [START_REF] Gross | Graph Theory and Its Applications[END_REF] where multiple edges between nodes and emergence of mass of junctions are permitted. We store the skeleton structural and morphometric features in data structures that link graph vertices to their adjacent edges. This solution allows us to manipulate the skeleton information and to implement a pruning procedure that takes into account not only segments lengths but also their average gray values in order to decide which branches are spurious and thus, delete them. The use of this contribution is illustrated in the context of medical application that cover almost all usable features. This work is composed of five sections. The first one presents skeleton features for medical application which covers biometric and character recognition structural features. The second section describes the used topological definitions for the graph construction. The third section is devoted to the graph data structures proposed in this work. The fourth section is dedicated to the implemented pruning and the manipulation of the structured graph skeleton on image of subchondral bone in the tibial knee with the aim of characterizing its microarchitecture.

Skeleton features in medical imaging

Concerning the characterization of bone microarchitecture, studies has shown that bone changes from a healthy person to a person with osteoporosis [START_REF] Sevestre-Ghalila | Texture Image Analysis for Osteoporosis Detection with Morphological Tools[END_REF] or osteoarthritis [START_REF] Chappard | Subchondral bone micro-architectural alterations in osteoarthritis: a synchrotron micro-computed tomography study[END_REF]. These changes are quantified using morphometric parameters of bone microarchitecture such as the number of pixels of the skeleton, the halfwidth and the length of trabeculae, the number of trabeculae and the number of nodes and ends. These morphometric parameters covers also structural descriptors (detection of nodes and extremities) used in biometric and character recognition applications [START_REF] Zhu | Hand Dorsal Vein Recognition Based on Hierarchically Structured Texture and Geometry Features[END_REF] [2] [START_REF] Wshah | Segmentation of Arabic Handwriting Based on both Contour and Skeleton Segmentation[END_REF]. We notice that skeletonization methods used to extract features for applications cited above are binary or gray. Consequently, we choose to make this feature extraction tool applicable on both types of skeleton. Furthermore, we notice that skeleton gray levels are not actually exploited. We propose to consider in the extracted features the segment mean gray level which is useful for the pruning of insignificant skeleton branches.

Identification of topological configurations

First of all, we need to identify the topological nature of each skeleton point in order to position it in the graph. A skeleton is formed of ridge, junction and end pixels. A ridge is a skeleton pixel placed on a crest line that splits the local background in 2 4-connected components. A junction splits the local background to 3 or 4 4-connected components and an end is a skeleton pixel with a unique skeleton neighbor. Our graph proposal is designed to accept as input any type of skeleton. To do so, it is important to identify critical situations that a conventional graph cannot handle. Indeed, a skeleton could have a mass of junctions that cannot be thinned. In Figure 1, pixels noted x and identified as mass of junctions divide the background in two 4-connected components according to their 8-neighbors. However, these pixels belong to a mass of nodes that should be entirely considered as a unique junction. In the literature, an erosion at nodes is performed before the construction of adjacency matrix in order to eliminate such configurations and to avoid a cyclic graph. A decomposition of the possible mass junctions in 3 × 3 configurations (with rotations) is made for our proposal in order to identify these pixels when processing the skeleton for graph construction. Once the ambiguous configurations identified, we can proceed with the construction of graph data structures.

Graph construction

Since we intend to conceive a graph that takes into account any kind of skeleton, simple graphs constructed using an adjacency matrix are not appropriate. A graph with multiple edges (multigraph) that manages mass of junctions and stores other morphometric features of skeleton branches is chosen in this work. We detail in the next the data structures used to construct the proposed multigraph and the linking step that permits to connect graph vertices and edges.

Data structures for the graph construction

First, we define in Table 1 P oint data structure that differentiates end, junction and ridge pixels and establish a primary link between current pixel and its direct skeleton neighbors. Having defined the point structure, we can define the Segment structure that corresponds to a graph edge. Intuitively, a segment has two extremities. There are 3 types of segments. The first one has two free extremities. The second has its extremities identified as junctions and the third one has a free end and a junction at its respective extremities. Therefore, we need to identify the segment extremities in the data structure in order to point out the segment type. The data structure segment is defined by the attributes listed in Table 2. Tail and head attributes correspond to segment extremities. These attributes are useful to establish links between graph edges. Since our solution is intended to be general, covering all types of skeleton graph, we consider segments having multiple heads and tails as illustrated by Figure 2.(a). Hence, we choose for the head attribute (respectively tail) an array of points of size 3. Furthermore, and unlike graphs based on adjacency matrix that cannot handle loops, we get our solution to consider a general case that defines more than one segment between two nodes. Consequently, we add two other attributes to the data structure, which are "Next to head" and "Next to queue" in order to distinguish the multiple segments appearing between two junction pixels. Such a case is illustrated by Figure 2.(b).

Linking step for segment structure

This module enables the implementation to fill point data structures and then, create segment structures by browsing the skeleton. It proceeds by first searching for graph nodes (junctions or ends) and second, by monitoring the ridge pixels of the same segment structure. Using the "neighbor" attribute stored in the point structure, it is possible to follow the segment and fill its remaining attributes. A stop condition of filling a segment is the detection of an end or other junction. This means that the linking step reached the end of the segment, can finally calculate its features and move to the next segment via the "tail" attribute. Figure 3 Finally, the flow chart of Figure 4 resumes the complete processing stages for the extraction of skeleton features. The graph construction is applicable to gray or binary skeleton resulting from any method of the literature. Skeletonization methods, and especially homotopic thinning are often subject to the emergence of spurious branches. For better reliability of results and pertinence of analysis, a pruning step is performed to eliminate these branches. In some applications, it is the orientation of the segment that indicates its "significance". However, orientation cannot be pertinent for microarchitecture quantification since trabeculae moves in various directions. Through experiments, we notice that the average gray level of spurious branches is informative and generally smaller than the average gray level of significant ones. We propose to implement a pruning procedure based on our graph construction that takes into account not only the length but also the mean gray value of segments to suppress noisy branches. In fact, according to the resulting skeletons of Figure 5 and their respective features, we can establish an upper threshold for segment lengths and gray values at a fixed percentile (50% in this case). These thresholds adequately adjusted permit to detect spurious branches and delete them in both cases. Fig. 5. Results of gray skeletonization using methods published in [START_REF] Youssef | Statistical control of thinning algorithm with implementation based on hierarchical queues[END_REF] and in [START_REF] Couprie | Grayscale image processing using topological operators[END_REF].

Color intensity (white to blue) reflects segments lengths.

An application of such process is performed on high resolution peripheral computed tomography (HR-pQCT) slice image from Voxelo project 1 as illus-trated in Figure 6. We note respectively N b.S, N b.N , N b.E, N b.seg the number of skeleton pixels, nodes, ends and segments, and N b.pix seg, Length of seg, Gray level of seg the mean pixelwise length, euclidian length an gray level of segments. We observe in the update of skeleton features a gain in pertinence after pruning spurious branches while preserving significant trabeculae.

Conclusion

The proposed graph construction applies to binary and gray skeletons and facilitate a high level manipulation of the skeleton via the structural and morphometric features it stores. The analysis of the trabecular microarchitecture aiming to differentiate normal bone from osteoarthritis is an example of skeleton graph structure usefulness. The pruning procedure exploiting graph and skeleton features eliminates efficiently spurious branches according to thresholds of the mean gray level and length of segments. This high level skeleton manipulation can be used in biometric and character recognition applications in the matching process. Our perspective is to create a GUI for easy handling of the graph: tracking lines, measuring distance between nodes/edges, updating interactively the features. In addition, other features such as branches orientation could be easily added to the data structure segment in order to monitor skeleton lines and add another optional threshold to the pruning.

Fig. 1 .

 1 Fig. 1. Topological configurations of skeleton pixels.

Fig. 2 .

 2 Fig.2. Particular situations for segment structure. In blue: segment end. In red: segment head (or tail) labeled as junction. In yellow: segment ridge point.

Fig. 3 .

 3 Fig. 3. Segments data structures for image from Figure 2.(b).

Fig. 4 .

 4 Fig. 4. Flow chart of features extraction tool.

  surface N b.S N b.N N b.E N b.seg N b.pix seg Length of seg Gray level of seg

Fig. 6 .

 6 Fig. 6. Skeleton pruning of the trabecular microarchitecture in a slice image of suchondral bone in tibial knee: Gain in relevance of features.

Table 1 .

 1 Data structure: P oint

	Structure P oint		Description
	identifier	pixel index in the image (unique)
	Flag	(1): end	(2): ridge	(3): junction
	Neighbors	only one 8

-connected neighbor Table of size [2..4] Table of size [3..8]

Table 2 .

 2 Data structure: Segment

	Structure Segment	Description
	Label	segment identifier
	Head	Table of point indentifiers for the segment head (maximum size is 3).
	Next to head	Index of the point that follows the head in the segment
	Tail	Table of point indentifiers for the segment tail (maximum size is 3)
	Next to tail	Index of the point that follows the tail in the segment
	Stack	stack containing segment pixels (index of ridge points)
	Segment features	Number of pixels, Euclidean length, Mean gray level, Half width

ANR Voxelo TECS-0018

Acknowledgments

This work is financed by ANR Voxelo TECS-0018.