Ludovic Moreau
email: ludovic.moreau@univ-grenoble-alpes.fr

Jérôme Weiss

David Marsan

Accurate estimations of sea-ice thickness and elastic properties from seismic noise recorded with a minimal number of geophones: from thin landfast ice to thick pack ice

come     Accurate estimations of sea-ice thickness and elastic properties from seismic noise recorded with a minimal number of geophones: from thin landfast ice to thick pack ice

Introduction

The propagation of seismic waves guided in sea ice has been exploited for decades to develop methods for the monitoring of the ice properties [START_REF] Anderson | Preliminary results and review of sea ice elasticity and related studies[END_REF][START_REF] Marsan | Sea-ice thickness measurement based on the dispersion of ice swell[END_REF][START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF][START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF], which are ingredients for climate and sea ice models.

These methods exploit the dispersion characteristics of the guided modes that compose the wavefield. With appropriate forward modeling, an inverse problem can be defined to infer the ice thickness and elastic properties, based on a fit between the dispersion of the guided modes in the model and in the data. Such approaches are quite common, not only in geophysics, but also at the ultrasonic scale for nondestructive testing [START_REF] Mitra | Guided wave based structural health monitoring: A review Smart Mater[END_REF] or medical acoustics [START_REF] Bochud | Predicting bone strength with ultrasonic guided waves[END_REF] applications. As far as sea ice applications are concerned, the main challenge with such monitoring methods are the in situ logistics, which require the deployment of seismic antennae with many geophones, as well as the use of active sources.

Given the hostile conditions and the difficulty to access polar environments, these are considered to be the main limitations of such methods, despite their potential for accurate sea ice properties estimations. Therefore, the long-term monitoring of sea ice with seismic methods remains unlikely as long as autonomous systems with minimal deployment logistics can be used.

Thanks to rapid technological and methodological progress, these constraints are less and less limiting. For example, current technology allows miniature triaxial autonomous geophones to record continuously the seismic wavefield in polar conditions for more than 30 days, at a sampling frequency of 500 Hz. The next generations is expected to record for several months while continuously transmitting the data via satellite communication. Moreover, the need of active sources was recently made unnecessary thanks to analyses based on seismic noise interferometry [START_REF] Marsan | Sea-ice thickness measurement based on the dispersion of ice swell[END_REF][START_REF] Marsan | Characterizing horizontally-polarized shear and infragravity vibrational modes in the Arctic sea ice cover using correlation methods[END_REF][START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF].

In [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF], a passive seismic experiment was performed to prove the concept of accurately measuring sea ice thickness and elastic properties with seismic noise. An array of 247 geophones was deployed at the surface of the frozen Van Mijen fjord in Svalbard (Norway), to record the ambient seismic field between 1 and 26 March 2019. The frequencywavenumber dispersion curves of the guided modes propagating in the ice cover were extracted from both icequakes and the noise correlation function (NCF). These dispersion curves were inverted for sea ice thickness, Young's modulus and Poisson's ratio, yielding very accurate and robust estimates when compared with active seismic acquisitions. With this approach, the wavefield must be spatially sampled. This requires a seismic array with stations spacing that meets Nyquist's sampling criterion, with typically ∼ 50 stations. Although it is expected that about 20-30 stations can be sufficient to apply this technique [START_REF] Moreau | Monitoring ice thickness and elastic properties from the measurement of leaky guided waves: A laboratory experiment[END_REF], developing new methods able to achieve similar results with even fewer sensors is desirable to reduce in-situ logistics.

We here introduce a different approach where fewer than 5 stations are sufficient to obtain accurate estimates of sea ice thickness and elastic properties. The method combines noise interferometry for estimating the elastic properties, with a Bayesian inversion of the dispersion in the waveforms of icequakes for inferring ice thickness. It takes advantage of the oneto-one relationship that exists between the time-frequency spectrum of the waveforms on the one hand, and ice thickness and propagation distance on the other hand.

The method is first applied to thin landfast ice. Based on a very rich dataset recorded in 2019 in the Van Mijen fjord, in Svalbard [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF], we demonstrate its robustness, stability and accuracy by processing waveforms from many active sources and from icequakes (section 3.1). Then, the method is applied to icequakes recorded on thick pack ice in the Arctic ocean in 2007 [START_REF] Gascard | Exploring Arctic transpolar drift during dramatic sea ice retreat[END_REF][START_REF] Vihma | Advances in understanding and parameterization of small-scale physical processes in the marine Arctic climate system: a review[END_REF] , indicating ice thickness values that are consistent with in situ measurements and field observations (section 3.2).

Material and Methods

In a layer of ice floating on water, the elastic wavefield contains at least four fundamental guided modes: the quasi-Scholte (QS ), quasi-S 0 (QS 0 ), quasi-A 0 (QA 0 ) and S H0. When the product of the frequency by the thickness of the ice remains under 50 Hz•m, the QA 0 mode is not propagative and the wavefield can be approximated with the combination of the flexural wave, the axial wave, and the shear-horizontal wave [START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF].These are essentially the asymptotic behavior of the guided modes, when the displacement field across the ice thickness is considered linear for the QS mode and constant for the QS 0 mode [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF].

For a homogeneous waveguide, modal dispersion depends only on the product between the frequency and the thickness of the waveguide. In particular, for sea ice, the only dispersive mode under 50 Hz•m is the QS mode. The following thickness estimation procedure relies on a time-frequency analysis of this dispersion. More specifically, we minimize the misfit between the dispersion of this mode in recorded and synthetic signals. This inversion procedure is described thereafter.

Forward Problem: Phase-Shift-Based Propagation Model

Given a Young's modulus, E, and Poisson's ratio, ν, the dispersion of the QS mode is characterized by a time-frequency spectrum which shape corresponds theoretically to a unique combination of ice thickness and propagation distance. A larger propagation distance results in more distorted signals with a dispersion that is specific to the ice thickness. This is illustrated in figure 1 that shows the synthetic waveforms of the QS at a distance of 80 m (figure 1a) and 250 m (figure 1c) from the source, together with their associate short-time Fourier transform (STFT), shown in figures 1b and 1d. Here the source is assumed to produce a 1.5-cycle toneburst with a gaussian shape, and with a centre frequency of 10 Hz, which is representa-tive of the dominant frequency observed in most of the icequakes in our data. The ice layer has a thickness h = 0.65 m, Young's modulus E = 4 GPa and Poisson's ratio ν = 0.33, which is representative of the ice near the seismic array in the Van Mijen fjord [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF].

To simulate the propagation of the QS mode, we use a very basic and efficient model which applies, in the Fourier domain, the appropriate phase-shift to each frequency component of the spectrum. This phase-shift depends on the frequency and is equal to k( f ) × L, where k( f ) is the wavenumber of the QS mode at frequency f and L is the propagation distance. To compute the wavenumber of the QS mode, we use the asymptotic low-frequency model introduced by [START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF], that gives the phase velocity of a time-harmonic flexure wave in an ice layer floating on an infinite water column.

It is noteworthy that sources in the field have complex and varied mechanisms, which have an influence on wave polarization, azimuthal directivity and signals amplitude. This results in signals that are generally quite different from a toneburst. The above-mentioned model by [START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF] is only for calculating phase velocities, which are not dependent on source mechanisms. Hence, source mechanisms are not accounted for in this model. Moreover, the amplitude information can also be modified along the propagation, for instance if the wave encounters a lead or a ridge. Such field-dependent features cannot be accounted for either, in our basic model. Therefore, in order to mitigate their effect through the inversion procedure, the signal in the model is modified as follows:

1. generate a 1.5-cycle toneburst with a gaussian shape and a center frequency of 10 Hz; 2. Fourier transform this signal; 3. replace the amplitude of the spectrum with that of the signal recorded at the geophone; 4. inverse Fourier-transform the modified spectrum.

This results in a pseudo-impulsive, broadband signal with a spectrum which amplitude is more representative of that from the recorded signal. This signal is then propagated between the source position and the receiver via the above-mentioned phase-shifting operation. In comparison with a non-modified toneburst, the modified propagated waveform will have a reduced misfit with the STFT of the recorded signal.

Parameterization of the Problem and Cost Function

In the following, we assume that E and ν can be estimated a priori, and we define a cost function between the data, d, and the output of the model, m(X), based on the time-frequency analysis of the flexural wave, such that

f (d, X) = 1 - 1 N N n=1 corr (STFT{d n }, STFT{m n (X)}) (1)
where corr denotes the correlation coefficient between matrices STFT{d n } and STFT{m n (X)}.

d n is the temporal signal recorded on the vertical channel at station n (n = 1, 2, ..., N) when an impulsive seismic source is activated, for example an icequake. m n (X) is the temporal signal simulated at station n, based on the input parameters, X. The parameters are representative of the wave propagation problem:

• source latitude

• source longitude

• ice-thickness, h

• a time-shifting parameter, ∆t, that accounts for the unknown activation time of the source during the recordings. The aim of this parameter is to "align" the simulated and recorded waveforms. It is the same for all stations. This problem is very well-constrained because there is a one-to-one relationship between the model parameters and the global minimum of the cost function. Hence in theory it should be possible to estimate sea ice thickness based on a single signal recorded by a unique seismic station. In practice, however, this can hardly be achieved because the waveforms are corrupted by ice heterogeneity, thickness variations, the presence of a snow layer, anisotropy of elastic properties etc. Nonetheless, the impact of these uncertainties can be mitigated by including data from a few more seismic stations. This significantly improves source localization, which simultaneously also improves thickness estimation. The resulting inferred ice thickness corresponds to an average over all the direct paths between the source and the stations.

Inverse Problem

To solve the inverse problem, we proceed with Bayesian inference, which provides an ensemble of solutions that fit the data with an acceptable level of likelihood, given the data uncertainty. This ensemble of solutions is represented by the posterior distribution of the model parameters, such that P (X|d) = P (d|X) P(X) P(d) .

(2) P (X|d) is the likelihood function, P(X) is the prior distribution and P(d) is the marginal likelihood function, or the model evidence. The posterior distribution expresses the conditional probability of the parameter values based on evidence from measurements, expressed by the likelihood function, and from prior assumptions, expressed by the prior distribution.

The Markov Chain Monte Carlo (MCMC) algorithm is one of several methods that can be used to evaluate the Bayesian posterior distribution by realizing a Markov chain in the parameter search space. This random walk satisfies the ergodic theorem, which allows the algorithm to converge toward a stationary state that approximates the probability density function (PDF) of the parameters [START_REF] Andrieu | On the ergodicity properties of some adaptative MCMC algorithms[END_REF]. A recurrent problem in Bayesian inference is the difficulty to calculate the marginal likelihood, which is essentially a normalization factor. However, because it is the same for all probabilities, its determination is not necessary in practice, since it can be cancelled by comparing ratios of probabilities rather than absolute probabilities. This is one of the motivations for using methods based on stochastic sampling such as MCMC.

In the present problem, it is assumed that measurement errors are uncorrelated and random, and thus that they can be modeled by a normal distribution. A zero-mean Gaussian likelihood function with variance σ 2 is therefore prescribed:

P(d|X) = exp - ( f (d, X)) 2 2σ 2 , ( 3
)
where σ 2 is the variance associated to the measurement errors. This is a typical likelihood function used in many data fitting problems [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF]. Moreover, it is also assumed for the prior distribution that the model parameters have equal probability over a finite range of values:

• the position of sources, which mainly originate from leads, ridges, or the shore line, is within a distance of 2 km around the centre of the arrays.

• ice-thickness is comprised between 0.1 m and 5 m.

• the phase-shift adjustment parameter is between -t max and +t max , where t max is the size of the time window where the signals to invert are recorded.

MCMC methods generally require a burn-in phase before reaching the posterior distribution. For improved convergence, we precede the MCMC algorithm by a Simulated Annealing (SA) global optimization. The number of iterations is set to 10000 in the SA and to 100000 in the MCMC algorithm. To approximate the PDF of the parameters, we sample 1000 candidates from the posterior distribution of each parameter. Thanks to the ergodic theorem, the mean of the PDF returns the expectation of the parameters, which is also considered here to be the estimated value, because the PDF should be centered around the zone with the highest probability in the parameter space. Another interest of the MCMC algorithm is to provide a confidence in the solution, based on the shape of the PDF.

3 Results

Landfast Ice in the Van Mijen Fjord in Svalbard

The data processed in this section were recorded in the Van Mijen fjord near Sveagruva, in Svalbard [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF]. Figure 2a shows the location of the deployment in the fjord.

The array is located near Sveagruva in a part of the fjord that is surrounded by a moraine and connected to the fjord by a channel (figure 2b). The array contains a total of 247 FairFieldNodal Zland geophones (flat frequency response down to the cutoff frequency of 10 Hz for 1C instruments and 5 Hz for 3C instruments, and an attenuation of -20 dB per decade below), as

shown in red and blue in figure 2c.

Elastic Properties of the Ice

Prior to applying the inversion procedure to the data, the Young's modulus and Poisson's ratio of the ice must be determined. To this end, we evaluate the velocity of the guided QS 0 and S H 0 modes from the noise correlation function (NCF). In passive seismology, the NCF is calculated by correlating the ambient seismic noise (or ambient seismic field) recorded between pairs of stations. It can be shown that it converges toward the impulse response, or Green's function, of the medium [START_REF] Sabra | Extracting timedomain Green's function estimates from ambient seismic noise[END_REF][START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF]. As shown in [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF], this approach can be applied to seismic noise recorded on sea ice to recover the Green's function of the wavefield propagating between the stations of a seismic array. These NCF are shown in figure 3 for the 10 combinations of station pairs as a function of the distance between the stations. They indicate an average velocity of 2170 m/s for the QS 0 mode and 1235 m/s for the S H 0 mode. These can be related to Young's modulus and Poisson's ratio such that [START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF]:

ν = 1 -2 c QS H 0 c QS 0 2
, and

E = ρ × c QS 0 2 1 -ν 2 . ( 4
)
This gives E = 3.8 GPa and ν = 0.35, which is consistent with the values reported in Moreau et al ( 2020), E = 3.9±0.2 GPa and ν = 0.34±0.02, which were estimated from a frequencywavenumber analysis based on the recordings of 52 geophones. The density is set to 900 kg/m 3 , a common value for sea ice. From a rheological point of view, it would be of interest to monitor the evolution of these mechanical properties between 1 and 26 March. Preliminary results indicate that Poisson's ratio keeps a conservative value between 0.32 and 0.33, and Young's modulus remains stable between 3.8 and 4.5 GPa. However this is out of the scope of this paper and is left for separate study.

Next, we apply the inversion method to both active sources (jumps from 1-m height onto the ice), and passive sources (icequakes). The active sources are ideal for a preliminary investigation of the accuracy of the inverted parameters, because they produce signals with a better signal to noise ratio, and also because their exact location is known. We present inversions for these two types of sources, based on signals recorded with 3 and 5 stations.

Inversions of Active Sources For Ice Thickness

On 1 and 26 March 2019, a series of active acquisitions were performed by jumping directly onto the ice from a one-meter height, near the linear arrays to the north, east, south and west of the main array (figure 2c). In total, 16 jumps were performed on both dates, each at a distance of about 0.5 m from the geophones in the linear arrays. Signals were extracted from the continuous recordings in a 5 s time window following the jumps. We refer the reader interested in the exact jumps time to [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF], where they are all listed. The inversion relies on a forward model that is valid when the product of the frequency by the thickness of the ice remains under 50 Hz•m. Since the ice thickness in the fjord was less than one meter, the extracted signals were band-pass filtered in the [1-50] Hz frequency band.

Although the landfast ice in the fjord is more regular than drifting ice in the open sea, it exhibits significant heterogeneities of mechanical properties through the thickness, as well as local thickness variations by up to 20 cm [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF]. These heterogeneities and thickness variations induce an apparent anisotropy in the propagation of the guided waves. For example, figure 4 shows the waveforms recorded at the main array in the two seconds following a jump at the north linear array. The waveforms are sorted with respect to the increasing distance between the stations and the jump location. The 231 stations of the main array are distributed on a square grid. Hence this sorting implies that from one waveform to the next, the propagation direction can be very different (±15 • ) while the propagation distance is almost the same. Although the overall moveout exhibits an increasing arrival time, there are also patterns where arrival times are not increasing monotonically, which does not happen when choosing stations that are all in a same line. This is typical of an apparent anisotropy.

•

Robustness of the Inversion

In order to check the robustness of our inversion procedure with regard to these uncertainties, we performed 30 successive inversions of the same source (jump near the northernmost geophone) at two different dates on 1 and 26 March 2019. Each inversion was performed with a different set of geophones that were randomly selected amongst the 247 available. The manuscript submitted to Journal of Geophysical Research: Oceans results of the inversions are given in table 1 when using 3 and 5 stations for the inversions.

These results show that the method is very robust to the apparent anisotropy, with accurate location of the sources, and thickness estimates with a standard deviation of less than 6.5 cm.

It is also noteworthy that when using 3 stations, we identified two cases where the position of the source was less constrained. This occurred in cases where the three stations were close one another (less than 15 m away), thus creating artificial far field conditions which prevent the azimuthal direction to be constrained. Such configurations were not encountered with 5 stations. Next, we also check the stability of the thickness estimations when applied to sources located nearby each other. To this end, once again the four sources at each linear array are used on both dates. Inversions are performed with 3 or 5 stations (figure 5). When using 5 stations, sources to the north, east, south and west are located very accurately, as shown in figure 5a.

At both dates, the ice thickness estimates are all within a 4 cm range for each of the four source areas (figure 5b). When using 3 stations, results have comparable source location accuracies.

In terms of ice thickness estimation, the values remain within a range of 6 cm, which is slightly larger than with 5 stations, but still very accurate, especially when considering that part of the variability is most likely due to the direction of propagation, and thus to actual spatial variations of thickness. Overall, the inversion can therefore be considered very stable since estimations from sources located in the same area give very similar values. Note the increase of thickness between March 1 and March 26. This was also observed in the field from ice drillings and ground penetrating radar acquisitions, which indicated an increase of 10-15 cm. Moreover, the thickness estimates are consistent with the values found when applying the method introduced in Moreau et al ( 2020). This method is based on a simultaneous inversion of the frequency-wavenumber spectrum of the QS , QS 0 and S H 0 modes.

This gave h = 54 cm with a standard deviation of 3 cm on 1 March, and h = 74 cm with a standard deviation of 5 cm on 26 March. The slight differences between the two methods are due to the fact that, on the one hand, the inversion in [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF] averages thickness variations along a short line of 50 stations (∼ 50 m) via a spatial Fourier transform. On the other hand, the present method averages wave propagation on larger propagation paths (∼ 250 m), each with very different directions between the source and the stations.

Inversions of Icequakes For Ice Thickness

The main objective is to achieve accurate thickness estimations from passive data only, so that long-term monitoring is made possible without the need of human intervention in the field, other than the deployment of the geophones. Therefore, the method is now applied to icequakes recorded on 5 and 24 March 2019. We attribute these icequakes to the presence of cracks that were observed in the field all along the shore line of the fjord. They are most likely produced by the mechanical stress induced by tidal forcing, and they were either closed or partially open, depending on the tidal phase. The waveforms of these icequakes are shown in fig-

ure 6 for three stations, S 1 , S 2 and S 3 (see figure 7), together with the corresponding synthetic waveforms resulting from the inversion with these stations. The average misfit after inversion is of the order of 1% for both icequakes, indicating that the data can be explained remarkably well with the simple forward model based on 4 parameters only.

The results of the inversions are shown in figure 7, when using 5 stations (figure 7a) or 3 stations (figure 7b), both indicating identical locations of the icequakes. The PDF is slightly more dispersed when only 3 stations are used. This was expected, since the more stations for the inversion, the better constrained the location parameters. The icequake on March 5 is located along the shore to the east of the array, and that on March 24 is located along the shore to the west of the array. This is consistent with the directions found from linear beamforming of the data using all stations of the main seismic array. This beamforming is shown in fig-

ures 7c and 7d as slowness of the waveforms versus azimuthal direction. Regarding the thickness, on March 5 the estimated value is 55 cm with a standard deviation of 2.5 cm in the PDF when using 5 stations, and 57 cm with a standard deviation of 3 cm in the PDF when using 3 stations. On March 24, the estimated thickness is 63 cm with a standard deviation of 2.7 cm in the PDF when using 5 stations, and 65 cm with a standard deviation of 3 cm in the PDF when using 3 stations. These values are consistent with those found when using active sources, as well as with estimations from the same icequakes when applying the method introduced in Moreau et al ( 2020), based on an inversion of the wavenumbers of the modes.

Pack ice within the Arctic ocean

The seismic data analyzed in this section were already described in [START_REF] Marsan | Low frequency bursts of horizontally-polarized waves in the Arctic sea-ice cover[END_REF]:

a seismic network was deployed in April 2007 as part of the measurement campaign at the Tara drifting station operating in the framework of the DAMOCLES (Developing Arctic Modeling and Observing Capabilities for Long-term Environmental Studies) project [START_REF] Gascard | Exploring Arctic transpolar drift during dramatic sea ice retreat[END_REF]. Figure 8a shows a part of this network made of four seismic antennae labelled Tromso, Paris, Tartu and Helsinki. Each antenna contains four short-period (1 Hz) vertical seismometers and one broadband Güralp CMG-3ESPC seismometers, installed in a diamond-shaped geometry with the triaxial seismometer in the centre. Figure 8b shows an approximate map of the situation in the field, with many leads around the four antennae, as well as the ice drilling and electromagnetic sounding positions relative to the antennae. The corresponding thickness profile between antennae Tromso and Tartu are shown in figure 8c, together with an histogram of thickness values (figure 8d), which indicates that the average thickness of the level ice is about 2 m.

For our analyses, we use one hour of continuous recordings on 1 May 2007 between 0:00 and 1:00, at antenna Helsinki (figures 8a and 8b). One of these recordings is shown in figure 9a, with a zoom between 0:30:30 and 0:36:00 (figure 9b) and the corresponding spectrogram (figure 9c). Hundreds of icequakes with energy up to 25 Hz are visible. In order to satisfy the validity criterion of the forward model, signals are low-pass filtered under 10 Hz. This allows thicknesses up to ∼ 5 m to be tested while keeping frequency-thickness values under the upper limit of 50 Hz•m.

Elastic Properties of the Ice

Because only one station of the Helsinki antenna is a 3C instrument, the method used in section 3.1 to evaluate Young's modulus and Poisson's ratio from noise interferometry is not possible here. We tried to exploit the horizontal channels of this station by including, in the cost function, the arrival times of the QS 0 and S H 0 modes. This would allow the velocity of the modes to be evaluated, based on the propagation time and the distance from the source, which could be used in turn for a joint inversion of ice thickness and elastic properties. However, this approach was not conclusive due to poor SNR in the data, which prevented the arrival times of these modes to be picked accurately. A potential solution consists in looking for icequakes originating from identical sources (also known as repeaters), via template matching methods. Stacking those waveforms could help tackle the SNR issue. However, the goal here is to demonstrate the potential of the inversion method when applied to a dataset acquired with the appropriate instruments, and there is no doubt that the elastic properties of the ice could be determined from the NCF, for example if the same geophones as those used on landfast ice had been used too. Hence, such ad hoc solutions are not investigated further here, since they are out of the scope of this paper.

For thick pack ice, reports in the literature give Young's moduli that vary between 6.1

GPa and 8.6 GPa, while values for Poisson's ratio are more conservative around ν = 0.33 [START_REF] Hunkins | Seismic studies of sea ice[END_REF][START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF]. In absence of in situ evaluations of the elastic properties,

we use E = 7.2 GPa and ν = 0.33, following the recommendation in [START_REF] Marsan | Sea-ice thickness measurement based on the dispersion of ice swell[END_REF], and the density is set to 900 kg/m 3 . Of course, choosing a different value of Young's modulus would change the value of the inferred ice thickness. However, this change would not be significant because the ice thickness parameter is much more constrained by the curvature of the STFT than it is by Young's modulus. Nonetheless, further investigations are required to quantitatively evaluate this sensitivity. This is left for a separate study, as this represents a considerable amount of inversions and analyses that are out of the scope of this paper.

Inversion of Icequakes For Ice Thickness

The vast majority of the recorded icequakes are not resolved in time, as shown in fig-

ure 9b, where one can see that the waveforms are partly superposed. This prevents a clear dispersion to be extracted from the time-frequency spectrum (figure 9c). The waveforms are also much noisier than those recorded in the Van Mijen fjord. We suggest that this is due to i) instrumentation issues, because SNR varies significantly between the stations, and to ii) a noisier environment, as shown in figure 9c where a low-frequency seismic noise with energy between 0.5 and 3 Hz is always present. Unfortunately, this is a frequency band where part of the useful dispersion information is present as well. This reduces the amount of icequakes suitable for an inversion to less than 1%. Yet, this should not be considered a fundamental problem, given the thousands of icequakes that trigger every day. Moreover, the thickness profile shown in figure 8c indicates significant local variations between 1.6 m and 7 m, due to deformed ice and pressure ridges (Haas et al, 2011). Such large variations likely modify the seismic wave propagation and cause scattering. With all the above-mentioned difficulties, this dataset represents the main challenge for our proof of concept.

To select icequakes suitable for an inversion, the following heuristic criteria were used:

1. the waveforms are clearly resolved in time and well-separated from those of other icequakes. 2. the SNR in the time-frequency spectrum allows the dispersion to be quantified between 0.5 and 20 Hz.

Based on these criteria, three icequakes were selected. The waveforms of two of these icequakes are shown in figure 10, together with the corresponding synthetic waveforms resulting from the inversion with 5 stations. The average misfit after inversion is of the order of 5% to 10%. This indicates that the data can be explained reasonably well with the forward model.

Figure 11 shows the corresponding inversions, performed with 5 stations (figure 11a) and less (figure 11b). The posterior distributions are shown in terms of the icequakes position relative to the array. They all indicate that the icequakes originate from the main lead located to the west of the array. With 5 stations, the estimated thickness was:

• 2.02 m for the icequake recorded at 0:02:06 with a standard deviation of 33 cm in the PDF;

• 2.06 m for the icequake recorded at 0:02:19 with a standard deviation of 32 cm in the PDF;

• 2.04 m for the icequake recorded at 0:28:47 with a standard deviation of 35 cm in the PDF. Note the non-gaussian shape of the PDF. This suggests a local minimum, in the cost function, that is associated to a thickness with relatively high likelihood around 2.5 m.

And with less than 5 stations:

• the inversion for the icequake recorded at 0:02:06 was not conclusive;

• 2.08 m for the icequake recorded at 0:02:19 with a standard deviation of 37 cm in the PDF, with 4 stations;

• 2.03 m for the icequake recorded at 0:28:47 with a standard deviation of 39 cm in the PDF, with 3 stations, although the icequake was not located exactly at the same position as when using 5 stations. Moreover, the PDF of the position parameters exhibits a wider spread. This is most likely a consequence of poor SNR combined with fewer stations used for the inversion.

In absence of thickness measurements around the array, we cannot give a definitive answer regarding the accuracy of these estimates. However, they are consistent with the measurements shown in figures 9b 9c which indicate an average thickness of ∼ 2m for the level ice, west of the main lead. These results are very encouraging, since we demonstrate that, despite poor SNR and significant local thickness variations that scatter wave propagation, icequakes can be located and ice thickness inferred with very good confidence. The standard deviation in the PDF are of the order of 30 cm, which indicates that the confidence in the solution is not as good as for the estimations from the previous dataset. This is a consequence of the larger misfit value (∼ 5-10% instead of ∼ 1%). We are confident that this misfit could be reduced by i) using instruments of the latest generation (for example those used on landfast ice), and ii) using a forward model able to account for the local variations in ice properties, such as a finite element model, but at the cost of significantly increased computational resources. However, we suggest that the actual standard deviation of the ice thickness estimation is much less than that in the PDF, since all inferred values remain within a 6 cm margin.

Concluding Remarks

A method was introduced to estimate the elastic properties of sea ice from seismic noise interferometry, as well as the ice thickness from the time-frequency dispersion of the flexural wave. This approach only requires passive recordings of seismic noise on 3 to 5 triaxial seismic stations (depending on SNR), to provide accurate and robust estimations of sea ice properties, while simultaneously relocating the icequakes. In thin landfast ice, we found an ice thickness between ∼ 55 cm at the beginning of March and ∼ 65 cm at the end of March, with a standard deviation less than 5 cm in the posterior distribution of the thickness parameter. In thick pack ice, we found an ice thickness of ∼ 2.05 m and a standard deviation of ∼ 30 cm in the PDF, but the values found from all inversions remain within a 6 cm interval, which suggests that the actual standard deviation is significantly less. All thickness estimations are consistent with measures made directly in the field (ice drillings, electromagnetic induction, ground penetrating radar).

Our investigations indicate that the proposed monitoring procedure exhibits similar performances when using 3 or 5 geophones, as long as the SNR allows a clear dispersion to be extracted from the waveforms of the icequakes. For example, thickness estimations remain within a 2 cm range in both cases. Of course, in configurations where 5 geophones are used instead of 3, the icequake location is better constrained, and consequently so is the ice thickness. However, the accuracy of the inferred parameters is comparable. Moreover, from the 60 array configurations tested in section 3.1.2, it appears that different array geometries lead to similar results, so long as the stations in the array are not all near one another. In data where SNR is poorer and the extraction of the dispersion is more difficult, using 3 geophones may not always be sufficient for a successful inversion. Hence for optimal monitoring results, we prescribe the deployment of arrays of 4 to 5 geophones per 1 km 2 areas, where the stations configuration is such that it maximizes spatial coverage.

Despite the demonstrated accuracy and robustness of our approach, reducing the misfit between synthetic and field data is essential to reduce the uncertainty in the solution. For this reason, future works should focus on two main directions. The first one is the improvement of SNR, for example via the automatic selection and denoising of icequake waveforms for optimal results. To this end, machine learning-based approaches such as clustering or template matching, have proved very promising. The second one is the use of efficient forward models able to account for local variations of ice properties while keeping the computational cost to an acceptable level.

Our estimations of the ice thickness correspond to an average over the propagation paths between the icequake origin and the seismic stations. It is therefore expected to be applicable at all scales of icequake propagation, that is, any distance up to a few kilometers in the frequency range of the present study. However, an interesting perspective of this work is the study of remote rupture sources, potentially hundreds of kilometers away, by investigating wave propagation at much lower frequency with stations spaced by ∼ 100 km. Most of what is known about sea ice deformation and fracturing comes from large-scale (> 10 km) satellite observations, such as the RGPS (RADARSAT Geophysical Processor System) dataset [START_REF] Kwok | The RADARSAT Geophysical Processor System[END_REF] or from the analysis of Lagrangian trajectories (buoys) [START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analyses[END_REF]. Satellite imagery allows to highlight that large scale deformation is mostly accommodated at the scale of a few days by linear structures of large dimensions (up to the scale of the basin). These structures are likely to be brittle [START_REF] Marsan | Scale Dependence and Localization of the Deformation of Arctic Sea Ice[END_REF]. The temporal resolution of how the rupture occurs along such long structures is still lacking, but could be analyzed by exploiting similar methods as exposed here.

The use of icequakes presents the advantage of being a completely passive approach, without the need of active sources. Given the very large number recorded every hour (>500 both in pack ice and landfast ice), one can fairly expect them to originate from a large range of directions and distances around the stations. This would open the way towards exhaustive, timedependent thickness tomographies of the ice, similar to the velocity tomographies of the crust.

Moreover, the exploitation of the horizontal channels may pave the way towards simultaneous inversions of the ice thickness and elastic properties, potentially with only two stations.

This can be achieved by including, in the cost function, the waveforms of the other two fundamental modes (S H 0 and QS 0 ).

A complementary perspective of this work is to apply the present method to the noise correlation function instead of icequakes. As shown in [START_REF] Moreau | Sea ice thickness and elastic properties from the analysis of multimodal guided wave propagation measured with a passive seismic array[END_REF], the correlation of seismic noise recorded at two stations for a few hours quickly converges towards the impulse response of sea ice between these stations. Since the relative position of the stations is known from GPS coordinates, the source-receiver station distance would be perfectly constrained, allowing for an even more accurate inversion of the ice thickness. The noise correlation function also integrates the contribution of all seismic sources, impulsive or more stationary. It naturally results in more energetic waveforms than those from icequakes, which has potential to extend the scales at which ice thickness can be inferred.

Figure 1 .

 1 Figure 1. (a) Synthetic waveform of the QS mode received by a geophone located 80 m away from an impulsive source in a floating layer of ice with thickness 0.65 m. (b) Short-time Fourier transform of the waveform. (c) and (d): same as (a) and (b) for a geophone located 250 m away from the source.

Figure 2 .Figure 3 .

 23 Figure 2. (a) Location of the seismic array in the Van Mijen fjord near Sveagruva (Svalbard), with (b) a zoom around the array area near Vallunden Lake, a part of the fjord that is surrounded by a moraine and connected to the fjord by a channel. The greyscale shows land which altitude is less than 25 m. All land above 25 m is shown in white to emphasize the shore line. (c) The 247 stations of the array, including the main central array, and the four linear arrays to the north, east, south and west. Red circles are for 1C stations and blue circles for 3C stations. The 5 blue squares indicate the stations used in section 3.1 to calculate the NCF. The large arrowheads indicate the positions of ice thickness measurements.

Figure 4 .

 4 Figure 4. Normalized signals received at all stations of the main array, versus the distance from an impulsive source located in the north linear array. Depending on propagation direction, the moveout of the waveforms is not increasing monotonically, revealing an apparent anisotropy of the wave propagation.

Figure 5 .

 5 Figure 5. (a) Inverted position of all active sources on 1 March (+) and 26 March ( ) 2019. Green circles represent the sources positions and black squares the stations used for the inversions. The four boxes show a zoom of the localized sources near each linear array. (b) Ice thickness estimates from these inversions. (c) and (d) Same as (a) and (b) with 3 stations only

Figure 6 .Figure 7 .

 67 Figure 6. Top: waveforms of the icequakes recorded in the Van Mijen fjord (black) and the corresponding synthetic waveforms (blue), resulting from the inversion with 3 stations (S 1 , S 2 and S 3 : see figure 7). Bottom: PDF of the ice thickness, with the mean value shown as dashed line.

Figure 8 .

 8 Figure 8. (a) Sketch of four the seismic antennae Tromso, Paris, Tartu and Helsinki, installed at the Tara ice station. The inset shows the geometry of one array, with typical inter-sensor distance. (b) Approximate map showing the situation in the field, with the position of the antennae relative to the leads, as reported by the Tara's onboard crew. (c) Ice thickness profile from drilled holes and electromagnetic induction between stations Tromso and Tartu. (d) Histogram of the thickness values.

Figure 9 .

 9 Figure 9. (a) One hour of seismic recording at one station of the HELSINKI antenna. (b) Zoom between 0:30:30 and 0:36:00 and (c) the corresponding short-time Fourier transform. The black arrow indicates a systematic seismic noise with energy between 0.5 and 3 Hz.

Figure 10 .

 10 Figure 10. Top: Waveforms of the icequakes recorded on pack ice (black) and the corresponding synthetic waveforms (red), resulting from the inversion with 5 stations. Bottom: PDF of the ice thickness, with the mean value shown as dashed line.

Figure 11 .

 11 Figure 11. Map showing the result of inversions for icequakes recorded on drifting sea ice at the Helsinki antenna. The diameter and color of the circles are proportional to the likelihood in the posterior distribution of the parameters. The squares indicate the stations used for the inversions. Solid lines show the approximative position of the leads relative to the stations. (a) Inversion of three icequakes with 5 stations. (b) Inversion of the icequake at i) 0:02:19 with 3 stations (black squares), and ii) 0:28:47 with 4 stations (black and red squares). The black dashed curves delimit the PDF of the icequakes position. The dashed curves separate the PDF of icequakes position.

Table 1 .

 1 Results of the 30 inversions of the active source near the northmost geophone (125 m away from

	center of main array)
	date	1 March 2019 1 March 2019	26 March 2019 26 March 2019
	number of stations for inversion	5 stations	3 stations	5 stations	3 stations
	average error on source position	4.6 m	5.2 m	6.5 m	6.3 m
	ice thickness estimation	50 ± 3 cm	52 ± 4.5 cm	67 ± 6.5 cm	68 ± 6.6 cm

• Stability of the Inversion

Acknowledgments

The metadata of the dataset recorded in Svalbard are available at the following doi: http://

dx.doi.org/10.15778/RESIF.XG2019. The full dataset (total size = 3.2 Tb) is available upon request to the corresponding author. A subset can be accessed instantly via a direct download link (25 Gb in size for each day). To request this link, please send an e-mail to resif -dc@ujf-grenoble.fr with your affiliation and IP address. This subset contains the continuous recordings of all stations on 1 and 5 March 2019, from which part of the results presented in this paper were obtained (with 3C stations data downsampled at 500 Hz). The dataset recorded on pack ice in the Arctic Ocean (a few Mb in size) is available upon request to the corresponding author.

ISTerre is part of Labex OSUG@2020. This research was funded by the Agence Nationale de la Recherche (ANR, France), and by the Institut Polaire Franc ¸ais Paul-Emile Victor (IPEV).

The authors would like to thank the associate editor and reviewers for the useful and constructive comments.


Stéphane D ' Alu
email: stephane.dalu@insa-lyon.fr

Oana Iova
email: oana.iova@insa-lyon.fr

Olivier Simonin
email: olivier.simonin@insa-lyon.fr

Hervé Rivano
email: herve.rivano@insa-lyon.fr

Demo: In-flight Localisation of Micro-UAVs using Ultra-Wide Band

In this demo we propose to use the UWB technology on micro UAVs as a way to estimate inter-drone distances. This information is the fundamental building block for a selfmaintaining formation flight of a UAV swarm. We present preliminary results for the computation of distance between UAVs and we show an UAV flying and computing its position in function of three fixed anchors.

Introduction

UAV swarms are usually controlled using high precision motion capture when indoor [START_REF] Preiss | Crazyswarm: A Large Nano-Quadcopter Swarm[END_REF], or Global Navigation Satellite System (GNSS) when outdoor [START_REF] Andrianarison | New strategy of collaborative acquisition for connected gnss receivers in deep urban environments[END_REF]. However, the hardware is expensive, energy-hungry, and for indoor, extra hardware needs to be installed in the infrastructure (so limiting its use in this area). We propose here the use of Ultra-Wide Band (UWB) technology for computing distance measurement between UAVs, which will then be used for estimating their relative position, and maintaining a swarm structure when flying. Combined with the IEEE 802.15.4-2011 PHY layer and its definition of ranging marker, UWB can perform distance measurement with a precision down to 10 cm within a 250 m range. It combines good obstacle penetration, and resilience to both multi-path effects and interference from other wireless technologies [START_REF] Ye | Experimental impulse radio IEEE 802.15. 4a UWB based wireless sensor localization technology: Characterization, reliability and ranging[END_REF].

The main challenges in using UWB for computing the relative position of UAVs in a swarm are:

• Measure the distance between UAVs in flight with good precision. At least three packets need to be exchanged between two UAVs in order to compute the distance between them, while they are both in movement.

• Compute UAV position fast and efficient. Since distance measurements have to be performed with at least three other UAVs in order to compute a position (using multilateration) several packet exchanges have to be made between the UAVs that form the swarm, which can lead to collisions, or network saturation.

• Eliminate the use of any external infrastructure. Computing the precise location of an object usually needs the knowledge of at least three fixed points, known as anchors. The UAVs inside a swarm need to be able to compute their position relative to a leader, without any external help. We present here some preliminary results towards reaching our goal, which cover the first two challenges mentioned above: distance and position computation.

Experimental Setup

In the experiments that we present here we used the Crazyflie micro UAV (weight: 30 g, width: 10 cm)1 that integrates a standard Inertial Measurement Unit (accelerometer, gyroscope, magnetometer, high precision pressure), to which we added the laser ranger for the Z axis, and the Decawave DW1000 UWB radio module2 (see Fig. 1). We differ from the Bitcraze Loco setup, as our goal is to remove fixed anchors and compute distances inside the swarm.

Distance Measurement Between UAVs

In theory, distance between two UAVs could be measured by exchanging one packet and computing its propagation time (difference of time between its transmission and its reception at destination), if they both have the same clock (pico-second precision is required for measuring centimeters). Unfortunately, in our case, the clocks of the UAV present some manufacturing defects and frequency drifts, consequence of the cheap hardware, which makes this approach not possible.

To mitigate this problem, we use a Symmetrical Double-Sided Two-Way Ranging (SDS-TWR) protocol [START_REF] Pelka | Evaluation of time-based ranging methods: Does the choice matter?[END_REF] that performs a round-trip on each side of the measured end points and computes the time difference between the two packets, which increases complexity, and delay between consecutive measurements. Influence of Multi-path. We made a series of experiments in order to confirm the resilience of UWB technology to multi-path interference in our scenario. The only problems that we encountered were an increase in the measured distance values when the UAVs were near the ground level. If we consider the two-ray ground reflection model, the breakpoint where distance measurements start to be altered can be calculated in function of the distance from the UAV to the ground (h t and h r ), and the radio frequency ( f ):

d = 4π h t • h r λ with λ = c f
As we can see in Fig. 2, for a height of 12 cm and a UWB channel 2 ( f 4 GHz), the break point is at ≈ 2.41 m. As the UAVs do no usually need to fly at such low distances, this should not be a problem. Influence of UAV Orientation. As the DWM1000 antenna has its specific radiation pattern 3 and the Crazyflie hardware on which it is attached could affect the radio signal, we wanted to check the impact of the hardware orientation on the distance measurement. We compared the distance distribution, according to different orientations (0°, 90°, 180°, 270°) for UAVs situated 1 meter apart. The results show similar values, with a normal distribution centered near the expected distance (Fig 3). Consequently, the orientation effect is negligible for a micro-UAV of this size.

UAV Localization using UWB

To compute the localization of an UAV in the horizontal plane we need to know the distance to at least three different UAVs (or anchors) whose exact positions are known. Considering the presence of multiple UAVs in a swarm who need to compute their localization, and that using SDS-TWR does not make it possible to perform multiple measurements 3 https://www.decawave.com/sites/default/files/resources/DWM1000-Datasheet-V1.6.pdf 0.8 0.85 0.9 0.95 1 1.05 1. at the same time (compared to a simple packet broadcast), there will be a high number of packets exchanged between UAVs, so a mechanisms has to be put in place for avoiding packet collisions. We propose to use a token-based algorithm that schedules the order in which UAVs perform their set of distance measurements using SDS-TWR (Fig. 4). At the end of the SDS-TWR packet exchange, the UAV piggybacks in the last packet (which contains the calculated distance) its own position and the standard deviation (which will be zero for the anchors). Finally, we use a multilateration based location estimation in which this received information (i.e., the distance, position, and standard deviation) is fed to the Crazyflie's existing Kalman filter, which combined with the IMU information estimates its current position.

Demo Presentation

Our demo consists of at least four UAVs (depending on the in-situ conditions), out of which three will be fixed on the ground and used as anchors forming an orthogonal axis in a 4mx4m area. One or more UAVs will fly and compute its relative position with regard to the anchors (Fig 5). To avoid multi-path influence we fix the anchors at 12 cm above ground level and the altitude for the flying UAVs at 1 m. The video 4 (Fig. 6) illustrates the demo.

Figure 1 :

 1 Figure 1: Crazyflie + UWB

Figure 2 :

 2 Figure 2: Distance accuracy near ground level

Figure 3 :

 3 Figure 3: Distance distribution according to UAV orientation

Figure 4 :

 4 Figure 4: UAV state machine

Figure 5 :

 5 Figure 5: Demo setup with 2 UAVs 3 anchors

Figure 6 :

 6 Figure 6: Demo with 1 UAV and 3 anchors

https://www.bitcraze.io/crazyflie-2/

https://www.decawave.com/product/dwm1000-module/

Video: https://www.youtube.com/watch?v=wJsGsIPM14U

Acknowledgments

This work was supported by the European Commission under the H2020 project BugWright2 (871260): Autonomous Robotic Inspection and Maintenance on Ship Hulls and Storage Tanks.


Matthieu Darnet

Emmanuel Godoy

Daniel Sadarnac

Stephane Gautrais

Current Loop Stability Analysis of a VIENNA-type Three-Phase Rectifier for an Imaging System Power Supply Application

Keywords: VIENNA-type rectifier, stability margins, inner current loop control, input filter, resonance, system modelling

published or not. The documents may come

INTRODUCTION

High-power imaging systems (MRI, X-ray scanner,etc...) are powered from the Hospital power grid through a multi-converters which has the following functions :

• An AC/DC conversion to provide a DC voltage for main load (X-Ray Generators) ; • An isolation from grid; • A DC/AC conversion to provide an AC threephase voltage for auxiliary loads; The AC/DC Converter has to handle large constraints:

• Step power from a few kilowatts to a hundreds kilowatts from main load; • Transparent mode operation despite a wide range of nominal input voltage and input grid impedance.

In order to improve image quality, these systems need accurate dc voltage regulation to increase the power and speed for capturing images.

There are also additional industrial constraints on cost and volume which must be taken into account when designing a new power supply.

Until now, the AC/DC converter is a passive threephase rectifier. And the low-frequency three-phase transformer which provides the isolation has a big volume. Consequently the output DC voltage is unregulated, the power supply has a low efficiency and a large volume, although its cost is low. Thus it limits the available peak power and the slope of pulsed load.

Active three-phase unidirectional rectifier addresses output voltage regulation, high efficiency, low volume. It enables also the use of a smaller high frequency transformer. Among the active rectifier, the Vienna topology has a regulated output, good reliability, good power density, a low total harmonic distortion and is well documented [START_REF] Leibl | Three-Phase PFC Rectifier and High-Voltage Generator for X-Ray Systems[END_REF] [START_REF] Kolar | The essence of threephase PFC rectifier systems[END_REF]. A prototype of an modified Vienna (cf. Fig. 3) has been successfully built by GE Healthcare and validated for several operating points.

In order to validate this topology and its control strategy, this paper investigates the stability margins on a wider range of operating points. The control model of Vienna (cf. Fig. 5) is made of three control loops on input currents, total output voltage and output midpoint voltage. The main stability issues appears in the input currents control loop because it is the fastest loop. Consequently only the current loop will be studied. A single variable approach has been chosen to clearly highlight the main issues regarding resonances and stability margins.

The paper is organized as followed. Section 2 describes a model of the modified Vienna Rectifier and its control strategy. Section 3 draws the stability margins of the current loop over one phase from a linearised averaged model. A sweep of output power, input nominal voltage and grid inductance shows the weakest operating points. Section 4 shows a validation of the model by comparison with measurements of the total output voltage and of one input current under a power load step. Section 5 concludes the main finding of this paper and discusses future work.

RECTIFIER MODELLING AND CONTROL

Fig. 1 shows how the power is supplied to the imaging system. A transformer makes a connection between power grid and hospital's grid. It limits the available power to 150kVA.

Rectifier Topology

Rectifier topology (cf. Fig. 3) is composed by two interleaved Vienna [START_REF] Kolar | The essence of threephase PFC rectifier systems[END_REF] and coupled boost inductances for each phase. This is why it has been named Coupled-Inductances Interleaved Vienna Rectifier (C2IV Rectifier).

A two parallel converters topology has been chosen to divide current constraints on semiconductors by two.

The input phase voltages V a ,V b and V c are ideal three-phase voltage source. The grid line between phase voltages and C2IV is modelled by an inductance L g in series with a resistance R g .

The input filter is composed by three capacitors C x in star connection with their parasitic resistors R x in series. The artificial neutral point of this filter is connected to the ouptut midpoint voltage M for common mode filtering.

Boost Inductances L of phase 1 and 2 are coupled. [START_REF] Costan | Convertisseurs Parallles Entrelacs : Etude des Pertes Fer dans les Transformateurs Intercellules[END_REF] shows the advantages of inductances coupling. Combined with an interleaved command of the two phase of the converters, current ripples in grid and in inductances see their amplitude divided by two and their frequency multiplied by two. Consequently, inductance size can be reduced for the same constraints in current ripple.

The chosen bidirectional switch is a T-type bridgeleg structure (cf. Fig. 4) with two IGBT and two antiparallel diodes. Its main advantage is low cost because of its usage for three-phase inverters.

Table 1 shows the value of hardware parameters of the C2IV prototype.

Rectifier Modelling

Assumptions made to obtain state space equations are identical to what has already been done for Vienna rectifier in [START_REF] Lai | Average Modeling and Control Design for VIENNA-Type Rectifiers Considering the DC-Link Voltage Balance[END_REF] and will not be detailed.

Every identical variables (currents, voltages, switch commands) are put together in vectors to simplify the writing of differential equations.

For each inductance and each capacitor corresponds a state variable as seen in Fig. 3. There are 14 state variables.

Input currents in boost inductances L are shown in eq.1 :

I abc12 = (I a1 I b1 I c1 I a2 I b2 I c2 ) (1)
Output voltages across the output capacitances C are shown in eq.2 :

V c12 = (V c1 V c2 ) (2)
Phase currents in grid inductances L g are shown in eq.3 :

I abc = (I a I b I c ) (3)
Input filter voltages across capacitances C x are shown in eq.4:

V xabc = (V xa V xb V xc ) (4)
There are 6 control variables shown in eq.5, because each bidirectional switch is controlled by only one signal :

D abc12 = (D a1 D b1 D c1 D a2 D b2 D c2 ) (5)
In real-time equations, D abc12 stands for the switching states which take the value 0 or 1. In average equations,D abc12 stands for the duty-cycles .

Perturbation variables are the three phase voltages (cf. eq.6) and the output current (cf. eq.7) :

V abc = (V a V b V c ) (6)
I out = (I o I o ) (7)
Matrices of coupling inductances (cf. eq.8), of duty-cycles (cf. eq.9) and of input currents sign (cf. eq.10) are used to further simplify writings of differential equations.

[L M] =     L 0 0 M 0 0 0 L 0 0 M 0 0 0 L 0 0 M M 0 0 L 0 0 0 M 0 0 L 0 0 0 M 0 0 L     (8) (1 -D abc12 )[i, i] = 1 -Dk i (1 -D abc12 )[i, j] = 0 if i = j with k i ∈ {a1, b1, c1, a2, b2, c2} (9) sign(I abc12 )[i, 1] = 1 i f Ik i > 0 else 0 sign(I abc12 )[i, 2] = -1 i f Ik i > 0 else 0 (10)
The equations of C2IV are finally obtained in eq.11

dI abc12 dt = [L M] -1 • -R L • I abc12 + I 3 I 3 •V xabc -(1 -D abc12 ) • sign(I abc12 ) •V c12 dV c12 dt = 1 C • sign(I abc12 ) T • (1 -D abc12 ) • I abc12 -I out dI abc dt = 1 Lg • V abc -Rg • I abc -(V xabc - 1 3 • ∑ V xabc ) dV xabc dt = 1 Cx • (I abc -[I 3 I 3 ] • I abc12 ) (11)

Rectifier Control

C2IV Rectifier control has a classic cascaded three loops topology, as shown in Fig. 5:

1. An inner loop on input currents I abc12 is used to control the Power Factor (PF) and Total Harmonic Distortion of input currents (THD); 2. An outer loop on total output voltage V dc is used to regulate V dc for a wide range of output power; 3. A third loop on the output midpoint voltage V mid is used to balance the two output voltages.

As one can often find in the literature, the outer loop controller is a PI with the peak input current as output. The third loop controller is a Proportional controller which adds an offset on duty-cycles. These loop are relatively slow.

The inner loop can be made in three different coordinate systems :

1. In dq-frame (I d , I q ) as in [START_REF] Liu | A control scheme of the VIENNA rectifier with unbalanced grid voltage[END_REF] [START_REF] Tang | An improved burst-mode control for VIENNA rectifiers to mitigate DC voltage ripples at light load[END_REF], [START_REF] Lai | Average Modeling and Control Design for VIENNA-Type Rectifiers Considering the DC-Link Voltage Balance[END_REF], [START_REF] Ji | Dual-loop Control for Threephase Vienna Rectifier with Duty-ratio Feedforward[END_REF] and (Liu et al., 2017a) to take advantage of a constant reference; 2. In αβ-frame (I α , I β ) as in [START_REF] Liu | Proportional-Resonant Current Control for VIENNA Rectifier in Stationary αβ Frame[END_REF] to reduce control variables without using a complex PLL like in dq-frame; 3. In abc-frame (I a , I b , I c ) as in [START_REF] Kolar | The essence of threephase PFC rectifier systems[END_REF] and [START_REF] Leibl | Three-Phase PFC Rectifier and High-Voltage Generator for X-Ray Systems[END_REF] for simplicity.

Usually, a PI controller is used in dq-frame. A PI controller in abc-frame has been chosen for this application to simplify controller tunning. Moreover, abcframe could be more robust to a phase loss because each phase is controlled independently. Double carrier-based PWM is usually chosen [START_REF] Liu | A control scheme of the VIENNA rectifier with unbalanced grid voltage[END_REF], [START_REF] Tang | An improved burst-mode control for VIENNA rectifiers to mitigate DC voltage ripples at light load[END_REF], [START_REF] Ji | Dual-loop Control for Threephase Vienna Rectifier with Duty-ratio Feedforward[END_REF] and (Liu et al., 2017a). However a single carrier-based PWM has been chosen for simplicity. It compares duty-cycles from the output of inner loop with a single carrier. This choice is made possible by the use of the absolute value of currents as shown in Fig. 5.

FPGA's capability allows a PWM sampling at 50 MHz. It suppresses the usual PWM sampling delay. Nevertheless the command, once changed by PWM algorithm, must be blocked until the end of half of the switching period to avoid command flickering.

The absolute currents reference in the inner loop is given by eq.12 :

I abc re f = I peak • |I sinus | (12)
The absolute sinusoidal shape is given by eq.13 :

|I sinus | = |V xabc | V xpeak (13)
As the measurements of voltage V xabc can be noisy, an observer has been preferred. The eq.14 shows the relation between input and output voltages and dutycycle for Boost-converter over one phase.

|V xabc | = (1 -D abc12 )[k, k] •V c12 (14)
The output voltages V c12 has been simplified by half of the reference value of the total output voltage

V dc re f 2 .
The peak input voltage V xpeak has been simplified by considering the maximum of the input voltages filtered measurements, which is about 95% of peak input voltage in three-phase. Thus the sinusoidal shape is built like shown in eq.15.

|I sinus | = V dc re f 2 • 1 max(V xabc ) • (1 -D abc12 ) (

CURRENT LOOP STABILITY ANALYSIS

Linearising model of eq.11 around an operating point, the linearised state space of C2IV Rectifier is obtained in eq.16 :

Ẋ = AX + BU Y = CX + DU (16)
with the definitions of vectors in eq.17 :

State Vector X = (

I abc12 V c12 I abc V xabc ) Command Vector U = (D abc12 D mid ) Output Vector Y = (|I abc12 | V mid ) (17)
Once this linearized state space obtained, one can draw the Bode Diagram between duty-cycles and input currents, study resonances and stability margins.

Input Filter Resonance Issue

Fig. 6 shows that the current loop bandwidth does not depend on phase voltage. Indeed, the three curves corresponding to input voltages (311V, -68V, -243V ) are superimposed for frequency above 1kHz.

There are 3 resonances related to hardware input parameters L g , C x , and leakage inductance L f (cf. eq.18) of coupled boost inductances. One zero z 1 (cf. eq.19) and two poles w 1 (cf. eq.20) and w 2 (cf. eq.21).

L f = L -M 2 (18) z 1 1 2π • C x • (L g + L f ) (19)
w 1 1 2π • C x • L f (20) w 2 1 2π • C x • L g •L f L g +L (21)
Resonances at w 1 (20) and w 2 (21) are poorly damped. They are located above 10kHz, near the switching frequency. It implies that the control model cannot efficiently damp these resonances.

Furthermore, resonances move the location of stability margins to higher frequency, nearer to instability.

A passive damping must be added to avoid potential oscillations and increase stability margins.

Switching frequency (cf. f sw in Table 1) limits already the current loop bandwidth at 11kHz. The presence of z 1 limits further the bandwidth at 2kHz for high grid inductance.

Control parameters of Table 2 have been chosen for this application.

The bode diagram of open current loop has been drawn in Fig. 7 combining transfer function of Fig. 6 and transfer function of PI controller (cf. eq.22). A gain margin of 14dB and a phase margin of 36deg are acceptable margins. Fig. 7 also shows that phase margin is calculated after high frequency resonances instead of around 2kHz where gain crosses first time at 0dB. The resonances are clearly not enough damped. As w 1 and w 2 depend mainly on capacitors C x and leakage inductance L f , a parallel damping circuit can be added to one of these components.

H PI (s) = K p 1 + 1 T i s (22)
A RC circuit in parallel with capacitors could be the best good option. Indeed, the resistance R damps the resonance. The capacitor C helps to decrease the damping circuit losses by filtering low frequency harmonics coming from the grid.

Fig. 8 shows the effect of a resistor R x of 250mΩ in series with each capacitor C x (cf. Fig. 3). Resonances are clearly damped. Phase margin is calculated here around 2kHz.

Nevertheless, damping three-phase input filter requires a study in itself and there is an extensive literature on this subject.

For the rest, only a damping resistor of 3mΩ has been added to take into account the parasitic resistors of the real capacitors implemented on the prototype.

Stability Analysis

The variable parameters of Table 1 modify stability margins. To understand where are the less stable operating points, these margins have to be calculated for all the domain.

The process of Fig. 9 have been used to get the stability margins for all the operating points. A simulation model has been used to find a stable operating point for each choice of L g , P out and V in . Then, the model has been linearised around these operating points. Transfer function between duty-cycle and in- The minimum values of gain and phase margins are acceptable.

However delay margins are particularly low. Indeed, a synchronized PWM would have a zero order hold (zoh) at twice the switching frequency (T zoh = 22µs) for interleaved rectifier. This would add a delay at twice the period of the zoh (delay zoh = 11µs), which is greater than the minimum delay on the domain.

Fig. 10 shows that stability margins hardly depend on grid inductance L g . Indeed, leakage inductance L f is very small compared to L g almost everywhere in the domain. In this case, resonances w 1 (cf. eq.20) and w 2 (cf. eq.21) are approximately equal and independent from L g . Fig. 11 shows that margins doesn't depend much on V in , particularly for a low inductance L g . For example, the gain margin range for a sweep of V in varies from less than 1dB to 5dB.

Fig. 12shows that margins depends more on P out . For example, the gain margin range for a sweep of P out varies from 7dB to 10dB.

The curves of phase and delay margins of Fig. 12 show a discontinuity between 40kW and 50kW for (L g ,V in ) = (10µH, 340V). It is explained by the fact that margins are calculated for a frequency above w 2 (cf. Fig. 7) at low power. However depending on operating points, the minimum margin can be calculated for a frequency below w 2 .

Hence C2IV control should stay stable at every operating point. However a better delay margin would be preferable. This can be achieved by the damping of input filter resonances.

Model Validation

A comparison between measurements over a prototype and a Simulink-based Simulation has been made to validate model of eq.11.

The parameters are taken from Table 1. Values of the variable parameters for the comparison case are shown in Fig. 13. A resistive load has been connected at the output of C2IV to make a power step from 1kW to 60kW. It allows a comparison of the system dynamic of the control of output voltage V dc in Fig. 14 and input currents I a1 in Fig. 15.

The absolute error between measurements and simulation data have been drawn to have a closer look at the differences between them. To remove harmonics due to the switching frequency, a filtering have been performed on the absolute error according to eq.23. A moving average over a window of 0.2ms has been used : 23) with N chosen as in eq.24 :

ε avg [ j] = 1 2N + 1 +N ∑ k=-N |X meas [ j + k] -X sim [ j + k]| (
N • T sampling = 1 × 10 -4 s (24) Fig. 14 shows that prototype measurements is very noisy. The error is up to 4.4V corresponding to the minimum of Vdc. This difference corresponds to some neglected losses which are not significant. Fig. 14 highlights that output voltage is not regulated well at no load. However voltage stayed within the +/ -50V safety boundaries. This is acceptable for no load condition. Fig. 15 shows the effects of coupling inductances for both simulation and measures. Current ripple amplitude is varying from nearly 0A to 30A depending on the duty-cycle. As a consequence, currents in boost inductances will enter Discontinuous Conduction Mode (DCM) at low power. In fact, stability margins have been calculated with the hypothesis of Continuous Conduction Mode (CCM) for all power. Thus a stability analysis at DCM operation must be performed for low power.

The mean error is oscillating between 2A and 10A. This error can be explained by the poor synchronization of the high-frequency current oscillations between prototype measurement and simulation data. That does not matter for this study. damped in reality. Nevertheless, these resonances have great chance to be excited. A known damping has to be added to prevent this case to happen. The differences are sufficiently low, which allows to validate the model of C2IV.

CONCLUSIONS

Current loop stability study has highlighted that input filter parameters of C2IV introduce potential resonances in the system and modify stability margins. Although gain and phase margins are acceptable for all operating points, delay margin is very low at low power and at high input voltage. Adding a passive damping is the preferred option to prevent oscillations and move the margins to a more stable location.

Thus stability of the current loop of Vienna-type rectifier has been proved on a wide range of operating points. Weaknesses related to input filter have also been identified.

Future work can improve stability study by using a multi-variable approach, like H∞ methods with small gain theorem and µ-analysis.

Stability of Discontinuous Conduction Mode is another field of study. Coupled Inductances makes this study a challenge for modelling current loop adequately.

C2IV Rectifier show promising robustness performance on all operating points. It can meet the increase demand of pulsed power for Imaging Systems.

Figure 1 :Figure 2 :

 12 Figure 1: power transmission chain between power grid and medical imaging systems.Load profile is shown in Fig.2. Power supply must face with low average power and high peak power.

Figure 3 :

 3 Figure 3: Schematic of the Coupled-Inductances Interleaved-Vienna Rectifier (C2IV) in a pseudo-3D view.

Figure 5 :

 5 Figure 5: Control block diagram of C2IV Rectifier.

Figure 6 :

 6 Figure 6: Bode diagram between the three first duty-cycles D a1 , D b1 , D c1 and their corresponding input currents I a1 , I b1 ,I c1 .

Figure 7 :

 7 Figure 7: Bode diagram of open loop current control.

Figure 8 :

 8 Figure 8: Bode diagram of open loop current control with input filter damping.

  put current have been obtained. Finally the margins of open loop transfer function are calculated and compiled.

Figure 9 :

 9 Figure 9: Process to obtain Stability Margins from the model and control of C2IV.

Fig. 10 ,

 10 Fig.10, 11 and 12 show that the worst operating point considering stability margin is with low grid inductance, high nominal input voltage and low output power. The margins for the point (L g ,V in , P out ) = (10µH, 520Vrms, 10kW) are :(Gain, Phase, Delay) = (8dB, 25deg, 3.5µs)The minimum values of gain and phase margins are acceptable.However delay margins are particularly low. Indeed, a synchronized PWM would have a zero order hold (zoh) at twice the switching frequency (T zoh = 22µs) for interleaved rectifier. This would add a delay at twice the period of the zoh (delay zoh = 11µs), which is greater than the minimum delay on the domain.Fig.10shows that stability margins hardly depend on grid inductance L g . Indeed, leakage inductance L f is very small compared to L g almost everywhere in the domain. In this case, resonances w 1 (cf. eq.20) and w 2 (cf. eq.21) are approximately equal and independent from L g . Fig.11shows that margins doesn't depend much on V in , particularly for a low inductance L g . For example, the gain margin range for a sweep of V in varies from less than 1dB to 5dB.Fig.12showsthat margins depends more on P out . For example, the gain margin range for a sweep of P out varies from 7dB to 10dB.

Figure 10 :

 10 Figure 10: Upper and lower bound of current loop stability margins for a sweep of L g .

Figure 11 :Figure 12 :

 1112 Figure 11: Upper and lower bound of current loop stability margins for a sweep of V in .

Figure 14 :

 14 Figure 13: Load step specification.

Figure 15 :

 15 Figure 15: (Up) Input current I a1 from prototype measurement and simulation data. (Down) Filtered absolute error between prototype and simulation current.

Table 1 :

 1 Simulation Parameters.

	Invariant Parameters	Units	Value
	Output DC voltage V dc	V	800
	Output capacitors C	µF	4240
	Boost inductance L	µH	800
	Boost mutual M	µH	760
	Intern resistance R l	mΩ	100
	Input capacitors C x	µF	10
	Parasitic resistance Rx	mΩ	3
	Grid frequency f grid	Hz	50
	Switching frequency f sw kHz	23
	Variable Parameters	Units Min	Max
	Ouptut power P out	kW	10	140
	Input AC voltage V in	V llrms 342	528
	Grid inductance L g	µH	10	300
	Grid resistance R g	mΩ	10	120

Table 2 :

 2 Control Parameters.

	Control Parameters	Value
	Current Loop PI Gains Proportional Gain K pI Constant Time Integral Gain T iI (s)	5 × 10 -3 4 × 10 -5
	Total Voltage Loop PI Gains
	Proportional Gain K pV Constant Time Integral Gain T iV (s)	2.5 1.3 × 10 -3
	Midpoint Voltage Loop P Gain Proportional Gain K pM	2 × 10 -3