Anders Riber Marklund
email: andersrg@kth.se

Frédéric Michel
email: michel.frederic@cea.fr

Henryk Anglart
email: henryk@kth.se

Anders Riber Marklund

Demonstration of an improved passive acoustic leak detection method on recordings from the Phénix steam generator operating at full power

Keywords: Sodium fast reactors, Acoustic leak detection, Fault detection algorithms, Hidden Markov models

teaching and research institutions in France or abroad, or from public or private research centers.   Demonstration of an improved passive acoustic fault detection method on recordings from the Phénix steam generator operating at full power

Introduction

Sodium fast reactors (SFRs) are fast reactor designs using liquid sodium as coolant. As, almost without exception, all nuclear power plants use steam turbines to drive the generator set, existing SFRs use or have used sodium heated steam generators to transfer energy from the secondary to the tertiary circuit [START_REF]Liquid metal reactors: Experience in design and operation[END_REF]. Sodium heated steam generators however present an inherent risk in the case of a leak between its two circuits. As sodium and water have a quite violent chemical reaction, releasing additional heat and corrosive reaction products, such a leak can damage other steam generator tubes nearby and create a self-accelerating damage process [START_REF]International Atomic Energy Agency, Acoustic signal processing for the detection of sodium boiling or sodium-water reaction in LMFRs[END_REF]. The task of detecting leaks has therefore always been important for SFR steam generator monitoring equipment, generally achieved through the use of chemical hydrogen detection systems with relatively slow response [START_REF] Jeannot | In-sodium hydrogen detection in the Phenix fast reactor steam generator: A comparison between two detection methods[END_REF]. Acoustic leak detection systems have been an interesting alternative for a long time due to their much shorter response times and relative simplicity but have never reached further than the experimental stage. Some overviews of work on the topic can be found in [START_REF]International Atomic Energy Agency, Acoustic signal processing for the detection of sodium boiling or sodium-water reaction in LMFRs[END_REF] and [START_REF] Matta | Supervision and Safety of Complex Systems[END_REF].

Within the ASTRID project [START_REF] Perrin | Cahier des Charges Fonctionnel du Réacteur ASTRID[END_REF], the Commisariat a l' Énergie Atomique et aux Énergies alternatives (CEA) is currently studying the option of using a Brayton cycle energy conversion system for SFRs with nitrogen as tertiary fluid. The steam generators would in this design be replaced by compact sodiumnitrogen heat exchangers. Even if the risk of violent chemical reaction is thereby excluded, hydrodynamic effects affecting both safety and performance of the heat exchanger, as well as the risk of a gas leak over to the primary system still motivate using a leak detection system as part of the continuous monitoring equipment.

The noises of small leaks will be due to aeroacoustic jet noise, fluid-structure interaction, sodium-water reaction (for the steam generator case) and the acoustic effects of bubble populations (absorption, scattering, diffusion and resonance) [START_REF] Yughay | Acoustic effects from water leaking into sodium[END_REF], [START_REF] Hur | Measurement strategy of the water leakage into a low pressure sodium boundary for a liquid metal reactor[END_REF]. The relative importance of these effects will depend on the leak rate, leak geometry, ambient conditions as well as the acoustic wave propagation path from the leak position to the acoustic sensor. Furthermore, neither measurement nor experimental simulation of realistic leak signals in a nuclear power plant are assumed to be possible. These facts taken together motivates aiming for a method that learns to recognize normal operation background noise and detects any deviation from this known noise.

In this work we develop such a method, based on a hidden Markov model concept proposed in [START_REF] Marklund | Passive Acoustic Leak Detection for Sodium Cooled Fast Reactors Using Hidden Markov Models[END_REF]. The performance is assessed by applying it to recordings from argon injection tests performed in the Phenix steam generator at full power. Section 2 contains a description of these tests and the recording data. Section 3 describes the improved hidden Markov model detector method and how detector performance is measured. The results, presented in section 4, starts with a short discussion of what is actually detected in the Phenix recordings. Then, a study of detector parametrizations is made, and the results of this study in detection of synthetic testing signals superposed on the background noise are used to choose a detector parametrization. A final demonstration is then made on the argon injection data as well as on pure background noise. Detection delay, false alarm rate and detectable additive SNR for the proposed method are all estimated.

Discussion of the method and our results as well as conclusions and some ideas for future work are given in section 5 and 6. In 2008 and 2009, a series of tests were carried out in the Phenix SFR plant with the aim of studying feasability of a passive acoustic leak detection system for the steam generators. During these tests, leaks were simulated by performing argon injections in a sodium inlet line to the steam generator. Five acoustic sensors of type ENDEVCO 7704 A were placed on sodium lines connected to the steam generator according to figure 1.

During a first acoustic measurement campaign in late 2008, 600 s of normal full-power background noise was recorded with a sampling rate of 25600 samples/s. Later, in May 2009, the actual recordings of argon injection tests were made, for some reason at a doubled sampling rate. Several injections of various durations and driving pressures were performed. In each test, the injection was started about 35 s after the recording start, injection durations ranged from 5 to 60 s and the recording length was 300 s in most cases. Six of the resulting recordings (including the pure background noise) have been made available by the CEA for use in this article. Characteristics of the analyzed recordings are shown in table 1 and the raw signals are shown in figure 2. Figure 3 shows the corresponding time-frequency evolution up to the lower Nyquist frequency of 12.8 kHz. For the signal on sensor 1, located closest to the injection point, figure 4 shows power spectral density estimates of the signal for background noise as well as for the 5 argon injection tests. The two 8 bar injections can be distinguished towards the end of the sensor 1 signal in figure 2. The other three injections are barely noticeable. Some impulsive behaviour is visible on the signals of sensors 2-4 and at the end of the sensor 5 signal. From figures 3 and 4 it is clear that the injections lead both to acoustic emission and absorption. The injection effects seen in figure 3 also seem to be somewhat delayed for each sensor, indicating that bubbles from the injection are detected as they are transported by the circulating sodium. This effect was also noted in an internal CEA document [START_REF] Descombin | Détection acoustique d'une réaction sodium eau (RSE) dans les générateurs de vapeur (GV) des RNR : Etat de l'art en[END_REF] and will be verified in some more detail in sections 4.1 and 4.4.

Detection method

Improved HMM detector

Hidden Markov models (HMMs) with Gaussian mixture models (GMMs) as emission probabilitites have been widely used in various speech recognition applications [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. The actual data modeled for human speech is often so-called mel frequency cepstrum coefficients, extracted from the measured audio signal in a way that imitates the functioning of the human ear [START_REF] Leijon | Pattern recognition -fundamental theory and exercise problems[END_REF]. For our application, i.e. fault detection in industrial acoustic noise signals, a new feature calculation scheme based on the Welch PSD estimate [START_REF] Welch | The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short Modified Periodograms[END_REF] was presented in [START_REF] Marklund | Passive Acoustic Leak Detection for Sodium Cooled Fast Reactors Using Hidden Markov Models[END_REF] along with some basic HMM theory. A brief review of this scheme explicitly stating the improvements that are new to this work will now be given.

We denote the PSD estimate of the time discrete signal x(t n ) with X(f j ), where t n are the signal sampling times and f j are the frequency bands. The average PSD of a normal operation signal is denoted X 0 . We monitor the residual XR from this model over sliding time windows of constant length tw, enumerated by m and we propose a componentwise normalization of the residual, i.e.

XR(f

j , m) = X(f j , m) -X 0 (f j ) X 0 (f j ) (1)
Contrary to the normalization used in [START_REF] Marklund | Passive Acoustic Leak Detection for Sodium Cooled Fast Reactors Using Hidden Markov Models[END_REF], this residual will not be indifferent to pure multiplicative changes in the signal, but it can be assumed that such changes can be identified by a parallell algorithm if needed.

In each time window, a vector of PSD component indices ordered by current residual size is also created, i.e.

J k (m) = j k |XR(f j k , m) > XR(f j k+1 , m) (2)
where the integer k enumerates its components. Based on XR(f j , m) and J k (m) we now introduce a modified feature calculation scheme compared to the one of [START_REF] Marklund | Passive Acoustic Leak Detection for Sodium Cooled Fast Reactors Using Hidden Markov Models[END_REF], given by

F(m) = F k (m) = J k (m) + β log j |XR(f j , m)| , k = 1, ..., k max (3)
The second term on the right represents the log of the total power deviation over all studied frequencies while the first term, for each k, points out the frequency band at which the PSD estimate is currently deviating the most, secondmost, thirdmost et c from its average value on the training data. The first term takes on discrete values from 1 to j max . The second term, which takes on real number values, is rescaled to cover approximately the same domain by calculating the factor β according to

β = j max 2σ XR (4)
where

σ 2 XR = E (d XR (m) -µ d XR ) 2 (5)
and

d XR (m) = log j |XR(f j , m)| (6)
i.e. σ XR represents the estimated standard deviation of the second term of equation 3 over the training data.

The relative importance factor of the first and second term of the feature, β, was in our previous work set by an ad-hoc trial and error procedure [START_REF] Marklund | Passive Acoustic Leak Detection for Sodium Cooled Fast Reactors Using Hidden Markov Models[END_REF]. Since we are aiming for a method to be licensed for safety related systems in nuclear power plants, basing the method on well-defined and predictable procedures is of utmost importance. The introduction of equation 4 in this respect represents an important improvement of the scheme. Furthermore, the new normalization of equation 1 has been introduced in order to allow considerably smaller amounts of data in the training sequence while maintaining detector performance. In [START_REF] Marklund | Passive Acoustic Leak Detection for Sodium Cooled Fast Reactors Using Hidden Markov Models[END_REF], training periods of about one minute were needed for a signal sampled at 200 kS/s. Using the scheme presented here, training periods of 20 s are sufficient, even on the signals with a sampling frequency of 25600 S/s, i.e. reduction by a factor of 24 of the amount of training data needed.

We now remind the reader that the model parameters of an HMM are set by applying a so-called training algorithm to an observed feature sequence F(m).

The training is done by local likelihood maximization using the Baum-Welch algorithm [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. When presented with a new feature sequence, it is then possible to calculate a conditional probability that the observed sequence is generated by the trained model by a recursive algorithm. In practice, the logarithm of this conditional probability is used in order to avoid computation problems due to small numbers [START_REF] Leijon | Pattern recognition -fundamental theory and exercise problems[END_REF]. I.e. after training a model λ bg on acoustic normal operation background noise, the following expression is used as a discriminant for fault detection

d HM M (m) = log(P (F(m)|λ bg )) (7)
In [START_REF] Marklund | Passive Acoustic Leak Detection for Sodium Cooled Fast Reactors Using Hidden Markov Models[END_REF], the model quality was studied as a function of the HMM model parameters only by reapplying the trained models to testing data on the same signal type. Here we will instead study the actual detection performance as a function of both the free model, and the free feature scheme parameters. The full set of studied parameters then becomes:

• The sliding window length tw

• The size of the feature vector, k max

• The number of Gaussians in the GMM, G

• The number of states in the HMM, S

Detector performance measures

A standard way of assessing the performance of a discriminant for detection is to calculate a receiver operating characteristic, or ROC, curve [START_REF] Fawcett | An introduction to ROC analysis[END_REF]. A common scalar measure of the detection performance is the area under the ROC curve. For detectors that have maximal area under the ROC curve, this measure can not be used to distinguish their performance. We therefore also use a min-max margin measure, similar to the one used in [START_REF] Marklund | Application of a new passive acoustic leak detection approach to recordings from the Dounreay prototype fast reactor[END_REF]. For discriminants d that are designed to decrease at the onset of detection, this measure is defined by

mmd = min d 1 max d 0 (8)
where d 1 and d 0 are the discriminant outputs in a known detection region and a known normal background region, respectively.

Results

For recordings 2-6, the training period used for each sensor channel extended from 5 to 20 seconds into the signal, for reasons that will be clear later. For recording 1, a training period from 120 to 220 seconds was used. The average background PSD estimate X 0 and the feature sequence F(m) used for training of HMMs was thus calculated on these periods.

Physical interpretation

It is instructive to study the nature of the acoustic signals by some basic change measure before going on to testing of the more sensitive HMM method. By listening to the actual recording from sensor 1 it was found that a tone-like sound appears a few seconds before the response in the higher frequencies. The proposed interpretation of this fact is that the injection itself creates a fairly weak low-frequency whistle and then the argon bubbles, a little later, create a high-frequency "trickle" noise as well as acoustic absorption as they reach the acoustic sensor position. As a simple demonstration, the discriminant defined by equation 6 was calculated for low (0 -2 kHz) and high (21 -26.5 kHz) frequencies on recordings 2-6 from sensor 1. The resulting log power deviations are shown in figure 5. We found that the response of the low-frequency discriminant appears 5-7 seconds before the high-frequency discriminant.

Next, the high-frequency deviation discriminants of recording 2 are shown for all sensors in figure 6. Here, it can be seen that the response is somewhat delayed for each sensor, and extended for sensors 2-5, corresponding to a hypothesis of argon bubbles traveling trough the steam generator and spreading out. The delays are larger from sensor 1 to 2 as well as from sensor 4 to 5. This corresponds well to the positioning of the sensors and the piping configuration (c.f. figure 1) where the bubbles pass through the longer S-shaped tubes only between sensors 1 and 2 and sensors 4 and 5. From sensor 2, at the reheater sodium outlet, to sensor 4 at the evaporator inlet, the bubbles travel through shorter piping, bypassing the superheater S-tubes.

The secondary sodium inventory of the Phenix plant was about 140000 kg per loop and the secondary sodium flowrate was about 800 kg/s. From these values, a sodium circulation period time for the secondary circuit of about 175 s is found. It is however reasonable to assume that some of the secondary sodium remained stationary in a storage tank and did not circulate, thus leading to a shorter period time. The fluctuations seen in the response of the first sensor about 150 s after injection start in figure 6 can therefore probably be attributed to argon bubbles returning after completing one lap in the secondary circuit. If the best detector would have to be found by a parametric study on injection noise data similar to recordings 2-6, the method would not be independent of such data. Instead, an approach where synthetic testing signals are superposed onto the normal background noise will be used here, i.e.

x(t 1 , ..., t 2 ) = x bg (t 1 , ..., t 2 ) + Cx test [START_REF] Descombin | Détection acoustique d'une réaction sodium eau (RSE) dans les générateurs de vapeur (GV) des RNR : Etat de l'art en[END_REF] where C is a constant factor and the testing sigal duration t 2 -t 1 was chosen to be 100 seconds. The first testing signal was a white Gaussian noise and the second one was a stationary signal with stochastic phase, created from a PSD of the background noise, using the same method as in [START_REF] Marklund | Application of a new passive acoustic leak detection approach to recordings from the Dounreay prototype fast reactor[END_REF]. The latter case is in some respect limiting, as it corresponds to detection of a stationary noise which has maximal resemblance to the background. For signals superposed in this way, an advantage is that the resulting signal-to-noise ratio can be controlled by the factor C and straightforwardly calculated from the definition SN R = 10 log 10 P test P bg [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF] where P test and P bg are the signal powers of Cx test and x bg respectively. A measure of the signal power is obtained by summing the PSD components up to the Nyquist frequency.

Putting ourselves in the position of having to parametrize the detector without acces to injection test data we now demonstrate how this can be done using only the synthetic testing signals and a simplistic parameter study. During this study, the background-like noise was added at SN R = -2 dB and the white Gaussian noise at SN R = -28 dB. The min-max distance is shown as function of time window length in figure 7. A local maximum is present at 2560 samples. Now restricting ourselves to parametrizations with tw = 2560 and S > 1, the min-max distance is shown as function of k max in figure 8. The performance can be seen to be split into two groups, corresponding to the two testing signals. For k max values larger than 1 there is no clear trend, but the performance seems to deteriorate at k max = 4, so we choose k max = 2 to promote smaller models. Finally, the performance as function of G and S is shown in figure 9. Larger models are seen to give slightly higher min-max distance. In a future application, more testing signals could be studied and a more sophisticated approach based on the response surface could be used to pick the parametrization with optimal performance. As the amount of training data was in some cases slightly too small for the larger values of G and the spread in performance seem to increase at S = 4, we here chose tw = 2560, k max = 2 G = 3 and S = 3 as the final parametrization. Since tw = 2560 yields 10 sliding windows per second at 25600 samples per second, we will denote this model by HMM10233.

Detectable SNR studies

In another study, the detectable SNR for the HMM10233 was studied, defining "detectable" as cases having a min-max distance larger than 1, i.e maximal area under the ROC. When analyzing the synthetic testing signals, the detectable signal power was found by gradually decreasing the superposition factor C of equation 9 and thereby also the SNR as calculated by equation 10. For the 8 bar argon injection signals, a changed SNR was first estimated on the raw signal according to SN R ∆ = 10 log 10 P ∆ P 0

where

P ∆ = |X 1 -X 0 | and P 0 = X 0 (12)
where X 0 and X 1 are the average PSD estimates in the background and the injection region of the signal respectively. By this method it was found that SN R ∆ ≈ -1.1 dB for the 8 bar injections. The detection performance at lower SNRs was then investigated by superposing additional stationary background noise onto the injection signal, as was done also in a study on injection test recordings from the PFR [START_REF] Marklund | Application of a new passive acoustic leak detection approach to recordings from the Dounreay prototype fast reactor[END_REF]. The resulting min-max distances as function of SNR are shown in figure 10. It is evident that the detectability is extremely dependent on the relative PSD shapes of background and the noise to be detected, a fact that was noted and demonstrated also in [START_REF] Marklund | Development and comparison of spectral methods for passive acoustic anomaly detection in nuclear power plants[END_REF]. Also, the background noise adding procedure used above gradually masks non-stationary fluctuations and creates a more stationary background which is deteriorating HMM detector performance for the 8 bar injection signals. We still found that these signals were detected also for a background noise that was increased by an additional 9 dB stationary noise. Taking into account the SN R ∆ this would correspond to an additive SNR of about -10 dB.

Detector demonstration

The HMM10233 detector was first tested on the pure background noise signals from recording 1. The responses (log conditional probabilities) on all sensors are shown in figure 11. To set a detection threshold, the following procedure was applied: As the output represents probabilities, it was fitted to a beta distribution. To be able to do so for these extremely small values, expressed on logarithmic scale, the following scaling relation was used, denoting the HMM log conditional probability output by d HM M (m)

D(m) = K * e d HM M (m) (13)
where K was a large constant. Then, the scaled output D(m) was fitted to a beta distribution and a threshold corresponding to a false alarm rate of 10 -8 was found by inspection of the cumulative distribution function. The obtained threshold was then inverse scaled back to the output region of d HM M (m). The same was done on all five sensor channels. The typical output level of d HM M (m) was around -150 for all channels when applied to a raw signal with Nyquist frequency of 12.8 kHz. The constant K was set to 10 54 and the resulting thresholds ranged from -307 to -316, indicating robustness of the approach. Due to impulsive events in the background noise, in total 5 false alarms on the 5 * 500 = 2500 discriminant samples outside the training region were observed, corresponding to a false alarm rate of about 0.002. As none of these false alarms lasted longer than 1 second; it would probably be possible to decrease the false alarm rate on the Phenix steam generator background noise to a value much closer to the value of 10 -8 obtained from the -316 threshold on the beta distribution fit.

The HMM10233 was then applied to the signal from sensor 1 for frequencies up to 12.8 kHz, using a threshold of -316, as found to be the minimum in the background noise study. The resulting discriminant is shown in figure 12. Noise changes which were judged to be due to argon bubbles circulating in the secondary circuit are detected abundantly after the first injection.

The HMM10233 was finally applied to the full injection test recordings (2-6) for all sensors and using full frequency data up to a Nyquist frequency of 25.6 kHz, resulting in the responses shown in figure 13.

Discussion

It is of importance to note that neither the ratio of driving argon pressure to the ambient pressure in the sodium line, nor the injector size were in these tests representative for a small leak. It has been demonstrated that apart from a low frequency whistle, the proposed acoustic anomaly detector here basically works as a bubble detector. A realistic small leak (around 1 g/s) would occur through a smaller orifice and with a higher pressure ratio since the secondary side of a steam generator is pressurized at 34-165 bar [START_REF]Liquid metal reactors: Experience in design and operation[END_REF]. The aeroacoustic noise from such a leak will therefore be higher, both in terms of frequency and power, than the injections performed here [START_REF] Zaman | Power spectral density of subsonic jet noise[END_REF], while the volume of bubbles produced will be smaller. Since the high detectability in terms of SNR of the gaussian white noise signal studied here is thanks to significant content in the higher frequency region where plant noise is relatively low, one could assume that real leaks would be detected with performance closer to this testing signal.

For the sodium-nitrogen heat exchanger developed by the CEA, the allowed detection delay can be extended to around one minute [START_REF] Perrin | Cahier des Charges Fonctionnel du Réacteur ASTRID[END_REF], thus making it likely that the proposed method will be sufficient for this application. In this case, using argon as a proxy for real nitrogen leaks is probably justifiable as both gases do not react chemically with sodium. For any unbuilt reactor system however, methods for reliable estimation of the acoustic SNR conditions of a small leak are however lacking. Also, if the detection system is required to function during operating point changes and transients, further development is also necessary. As the HMM concept is originally developed for classification of several different noise sequences (such as different words in speech recognition), it should be well suited also for this task but the resulting anomaly detection performance is not evident.

The procedure of parametrizing the detector used here (one parameter at a time) is rudimentary. The aim is to show that it can be done, using only background noise and relatively basic testing signals, and still arrive at a detector that performs well on real injection data. As stated in section 4.2 the procedure could be improved by using more testing signals and a responce surface method to find an optimal model. The parameter study is however time-consuming (on the order of several hours on the laptop computer used in this work) but in a real plant, it is reasonable to assume it needs to be performed only at first start-up. Re-training of the detector might have to be performed more often but should only take one or two minutes based on the computation times observed here.

Conclusions

A hidden Markov model method proposed earlier for passive acoustic leak detection in sodium fast reactor systems has been improved to eliminate ad-hoc setting of parameters and allow smaller amounts of training data. The improved method is verified to work well in a full power noise environment recorded during argon injection tests in the Phenix plant. The method is based on training a model on known background noise, i.e. without any assumptions on possible leak noise, and the free parameters of the model are set by a parametric study of detection performance for synthetic noises superposed on pure background noise.

The false alarm rate is governed not only by the detection method used, but also by the noise of the plant. For the 5x10 minutes of pure background noise studied here, the only sources of false alarms were impulsive events lasting about one second. By assuming an allowed detection delay of e.g. 3 seconds, it was thereby estimated that a false alarm rate close to the theoretical value of 10 -8 , can be achieved on a single sensor channel. Regarding the minimum signal to noise ratio for detection, it has been shown that the 8 bar injection signals can be detected at conditions corresponding to an additive SNR of at least -10 dB. These results seem reasonable compared to the conclusions of [START_REF] Marklund | Application of a new passive acoustic leak detection approach to recordings from the Dounreay prototype fast reactor[END_REF], and [START_REF] Marklund | Passive Acoustic Leak Detection for Sodium Cooled Fast Reactors Using Hidden Markov Models[END_REF] which were not based on realistic background noise and where the former work suffered from a lack of data preventing a credible false alarm rate estimation.

Possible ways of improving passive acoustic detector performance further would be to specifically monitor high-frequency regions where background noise can be expected to be relatively low (if leak noise can be expected in this region). Multi-sensor methods using adaptive filtering and/or spatial focalisation techniques such as beamforming also represent possibilities for further improvement.

In order to increase the maturity of this research field further, it is also of interest to study methods for estimating SNR conditions of leaks in a given system and to develop detection methods that will work also during operating point changes. Finally, even if the method proposed here is independent on knowledge on possible leak noise, such knowledge will remain of large interest to motivate monitoring of particular frequency regions and thereby increase performance further if needed.

Acknowledgments

This study was performed within the framework of an ongoing collaboration project between KTH and CEA on the instrumentation and safety of sodium cooled reactors, supported by the Swedish Research Council, grant number B0774801.

Figure 1 :

 1 Figure 1: Schematic view of the Phenix steam generator indicating approximate positions of the five acoustic sensors and argon injection system used in the 2008-2009 tests.

Figure 2 :

 2 Figure 2: Raw acoustic signals (acceleration) of all studied recordings with recording breaks indicated by vertical lines and responses of sensors 1 to 5 displayed from top to bottom.

Figure 3 :

 3 Figure 3: Time-frequency evolution of all studied recordings on sensors 1 to 5, displayed from top to bottom. (The common sampling rate and thereby the Nyquist frequency of all recordings is limited by recording 1.)

Figure 4 :

 4 Figure 4: Power spectral density estimates of the recorded signal at sensor 1.

Figure 5 :

 5 Figure 5: Upper panel: Low (solid line) and high (dashed line) frequency discriminants applied to sensor 1, recordings 2-6. Lower panel: Zoom on recordings 5 and 6.

Figure 6 :

 6 Figure 6: High frequency energy deviation discriminant applied to sensors 1-5, recording 2.

Figure 7 :

 7 Figure 7: Min-max distance of synthetic testing signals versus tw.

Figure 8 :

 8 Figure 8: Min-max distance of synthetic testing signals versus of kmax for tw = 2560.

Figure 9 :

 9 Figure 9: Min-max distance of synthetic testing signals versus G and S for tw = 2560 and kmax = 2.

Figure 10 :

 10 Figure 10: Min-max distance as function of additive SNR for WGN (triangles), 8 bar argon injection (diamonds) and background-like noise (squares).

Figure 11 :

 11 Figure 11: HMM discriminants on recording 1 (having a Nyquist frequency of 12.8 kHz). The HMM training region is indicated by a rectangle. The threshold was set by fitting a scaled version of d HM M outside the training region to a beta distribution and requiring a theoretical false alarm rate of 10 -8

Figure 12 :

 12 Figure12: HMM discriminant on recordings 2-6 for sensor 1, using a downsampled raw signal with Nyquist frequency 12.8 kHz in order to be comparable to the background noise study.

Figure 13 :

 13 Figure 13: HMM discriminants on recordings 2-6 using the full raw signals with Nyquist frequency of 25.6 kHz.

Table 1 :

 1 Studied recordings

	Recording Recording Injection P Ar [bar] f s [Hz]
	no.	length [s]	length [s]
	1	600.0	0.0	-	25600
	2	300.0	20.0	7	51200
	3	300.0	20.0	6	51200
	4	300.0	20.0	5	51200
	5	116.6	60.0	8	51200
	6	300.0	60.0	8	51200
	2. Experimental data
	2.1. Argon injection tests at the Phenix SFR

Table 2 :

 2 Studied parametersThe parameters studied are listed in table2. The values for the sliding window length are defined by numbers yielding an integer number of windows on one second of data (at the lower sampling rate of 25.6 kHz). The cases having S = 1 represent pure GMMs, i.e. they lack time evolution modeling and are technically not hidden Markov models.

	Parameter Range	Unit
	tw	800, 1024, 1280, 1600, 2560, 3200 [samples]
	S	1,2,3,4	[-]
	k max	1,2,3,4	[-]
	G	1,2,3,4	[-]
	4.2. Parameter study on synthetic test signals


Arnaud Lemiere

Frédéric Desevedavy

Bertrand Kibler
email: bertrand.kibler@u-bourgogne.fr

Pierre Mathey

Grégory Gadret

Jean-Charles Jules

Christophe Aquilina

Pierre Béjot

Franck Billard

Olivier Faucher

Arnaud Lemière

Frédéric Désévédavy

Frédéric Smektala

Mid

Keywords: chalcogenide glasses, optical fibers, supercontinuum generation, mid-infrared optics

come

Introduction

Many chemical species, as organic compounds and biomolecules, present fundamental vibrational resonances between 2 and 20 µm. The mid-infrared (mid-IR) spectral region is thus considered as a fingerprint region of huge interest for sensing applications [1], [2]. Particularly, greenhouse pollutants gases, such as CH 4 , CO 2 , NH 3 , and NO x , present absorption bands in this spectral region [1], [2]. Current issues related with air quality and climatic modifications drive [3] the interest for developing new sensing devices. Coherent tuneable large band sources emitting in this spectral region are promising for such applications.

Over the past decade chalcogenide-glass based optical fibers have revolutionized the field of mid-IR broadband light sources [4] [5] [6]. Chalcogenide glasses exhibit the suitable wide transparency window and high optical nonlinearity [7], making them the best candidates for midinfrared supercontinuum (SC) generation. Numerous achievements have been accomplished through the development of different kinds of fiber structures, such as step-index, microstructured and tapered fibers, in order to control light guidance according to the pumping conditions. In particular, the advanced methods of microstructured or suspended-core fibers and fiber tapering have been widely used to decrease zero dispersion wavelengths (ZDWs) while enabling enhanced effective nonlinearities. As a result, efficient mid-IR SC generation can be investigated by using relatively shorter pump wavelengths with lower pump powers, thus fitting the current needs of developing compact devices based on fiber laser pumps below 5 µm [8]. However, it is worth mentioning that the broadest SC spectra have been obtained in step-index fibers with larger transmission window (i.e., more robust and less sensitive to aging issues [9]) and by using longer pumping wavelengths delivered by high-energy optical parametric amplifiers [20] operating in the femtosecond regime. Table 1 presents a survey of mid-IR SC spectra obtained in chalcogenide fibers made of distinct glass compositions and index profiles since the ultrabroadband mid-IR SC generated in a step-index fiber by Petersen et al. in 2014 [10]. From Table 1, it clearly appears that, to the best of our knowledge, the broadest SC obtained with As-S or As-Se suspended-core fibers ends up to 9.5µm. A way to allow further mid-IR spectral broadening is to modify the glass system and its transmission window. Another solution is the control of the fiber geometry to manage the group-velocity dispersion. By combining these two features, a suspended-core fiber consisting of dual chalcogenide glasses (Ge-As-Se-Te / Ge-Sb-Se) has been developed, leading to the largest SC obtained in a chalcogenide suspended-core fiber, covering the 1.7-11.3µm spectral region with a pumping wavelength around 5 µm [19]. However, one can note again that step-index fibers based on a similar glass system allow to obtain a larger SC span [23].

Most of glass systems listed in Table 1 contain arsenic as well as antimony which are toxic elements. Especially, they lead to carcinogen oxides when they are exposed to atmosphere [25]. The European REACH (Registration, Evaluation, Authorization and Restriction of Chemicals) [26] regulation aims to limit as far as possible the use of highly toxic substances, and this drives the interest to search for As and Sb-free chalcogenide glasses. From that point of view, the GST (Ge-Se-Te) glass system is an alternative for nonlinear applications in the mid-IR region. GST glass system is known from 1966 [27], but it has been less studied for SC applications though its transmission window is similar to the As-Se glass [20] [21], mainly because the glass formation and fibre drawing is more tricky. Maurugeon et al. successfully drew a Te-rich GST fiber and then performed sensing application between 13 and 16 µm where CO 2 absorbs with the Fiber Evanescent Waves Spectroscopy method [28]. The composition used contains only of a few percent of Se to minimize absorption between 16-18 µm and presents a minimum of losses of 10 dB/m at 10.5 µm. The dependence of the refractive index with respect to distinct compositions of GST glasses was investigated in Ref. [29]. The control of Te and Se contents (i.e., from Se-rich to Terich) in the glass system allow to change the resulting refractive index over the range 2.5-3.4 in the mid-IR. In 2015, Jiang et al. proposed a purification process to obtain a low-loss Se-rich Ge-Se-Te chalcogenide fiber (Ge 20 Se 70 Te 10 ) [30]. The optical loss was found to be 1.4 dB/m at 4.2µm and less than 3dB/m between 1.5 and 6.3 µm. Recently, Sun et al. has synthesized Ge-Se-Te glasses and also studied the evolution of the nonlinear index n 2 as a function of the Te content [31]. It appears that n 2 increases with Te content. Based on this potential of Se-rich GST glass system, we have recently demonstrated the fabrication of optical fibers with distinct index profiles, namely, single-material fibers with or without a suspended core and standard step-index fibers with varying core diameter [20]. In addition, we have evidenced their potential for nonlinear mid-IR devices by means of SC generation spanning from 2 to 14 μm in a step-index profile [20]. In the same work, a preliminary study of suspendedcore fibers is also presented with a SC generation experiment spanning from 2 to 11.5 µm. This result then compete with the state-of-the-art of SC in suspended-core fibers shown in Table 1.

In the present work, our purpose is to present in more detail the development of As-Sb-free chalcogenide suspended core fibers based on the ternary glass system Ge-Se-Te. In particular, the glass composition Ge 20 Se 70 Te 10 is here successfully synthesized, partially purified and drawn to a suspended core fiber. The core optical transmission is measured whereas group-velocity dispersion and light confinement properties are computed and discussed. Finally, to the best of our knowledge, we report on the largest mid-IR extension of SC generation in a suspended-core optical fiber. More specifically, we demonstrate a SC spanning from 3.25 to 13 µm into a 12-µm suspended-core fiber pumped by sub-100-fs pulses at 8.6 µm. Experimental SC results are supported by numerical simulations of the nonlinear pulse propagation.

Glass synthesis, fiber fabrication and characterization

A Ge 20 Se 70 Te 10 glass preform is synthesized by the standard melt quenching method [32]. High purity elements of 5N bulky Ge and Te from Umicore and 5N Se pellets from Alfa Aesar are used. Stoichiometric quantities are weighted and then inserted into a silica ampoule. In this case we performed a slight purification: this ampoule is evacuated at 1.10 -5 mbar for several hours. Then we proceed to a single static distillation (see Fig. 1). This allows removing of refractory elements, such as Carbon or SiO 2 responsible for optical scattering [33] [34]. Finally, the ampoule is sealed and inserted in a rocking furnace for glass synthesis (see Fig. 1). The batch is heated up to 850°C for 10 hours to allow reaction between the starting products, then the melt is quenched in water at ambient temperature and finally annealed at T g -5°C for 12 hours.

The GST glass transmission is measured on a 4.7 mmthick bulk sample, previously polished with SiC papers and diamond suspensions. Mid-IR transmission is measured with a Fourier Transform InfraRed (FTIR) Perkin Elmer spectrometer (Spectrum One). Transmission in the near-IR region is measured with a Perkin Elmer Lambda 900 spectrometer. The transmission T of the glass is presented in The bandgap is measured at around 1 µm which corresponds to an absorption coefficient of 10 cm -1 . We can express the absorption coefficient with the following equation:

with ( 1
)
where I (resp. I 0 ) is the output (resp. input) intensity of the light and L the thickness of the bulk. Assuming Se-H extinction coefficient ε SeH = 1 dB/m/ppm(at) [33][36], we estimate the residual Se-H quantity at 10 ppm(at) which is of the same order of magnitude as previous studies [20][33][34] [37]. In the case of OH related absorption band, OH extinction coefficient is measured in silica as ε OH = 10 dB/m/ppm(wt) [38] and by considering such a value, there is 6 ppm(wt) of OH in our fiber. This value, above that of previous works [20][34] [37], is related to the fact that we only performed a static distillation for the synthesis of our glass, and not a whole purification procedure with various dynamic and static distillation steps [19] [32]. These considerations stands also for Ge-O and H 2 O absorptions [20][35] [37].

Next, a suspended-core fiber is fabricated from the mechanical drilling of a 16 mm outer diameter initial rod [39]. Three holes are drilled with 1 mm in diameter and 30 mm in length. During the drawing process, these holes allow the formation of the air-suspended structure with three fine bridges supporting the core. The drilled preform is indeed drawn to a fiber under a nitrogen gas pressure (1-4 mbar) to control the hole's size as described in the work of Savelii et al. [40]. The drawing process is occurring at a temperature, inside the furnace, of around 380 °C. Several meters of fiber with a suspended-core diameter Φ = 12 µm are obtained from a 6 cm-long preform (see Fig. 3). The refractive index of the glass was calculated using Clausius-Mossoti relation and the refractive index measurements of binary systems GeSe 4 and GeTe 4 found in the literature [41] [START_REF] Mittal | Proc. SPIE[END_REF]. The refractive indices based on bulk reflection measurement curves [29] are not accurate enough for numerical simulations. Here we consider our ternary Ge 20 Se 70 Te 10 glass as a glass obtained from the addition of 0.875 GeSe 4 and 0.175 GeTe 4 [20]. The wavelength dependence of refractive index is plotted in Fig. 4. It is well fitted with a two-pole Sellmeier equation (Eq. 2) and the corresponding Sellmeier coefficients are listed in Table 3.

(2) Group-velocity dispersion and effective mode area of the fundamental mode of the suspended-core fiber are then investigated through numerical simulations based on the fiber geometry derived from the scanning electron microscopy image and by means of a commercial software using a vectorial beam propagation method (RSoft product). Figure 4a shows the wavelength dependence of the dispersion for the fundamental guided mode, for distinct core diameters, as well as the dispersion of the bulk glass. The strong light confinement allowed by the suspended core geometry allows to tune the zero dispersion wavelength (ZDW) well below the intrinsic ZDW of the bulk material located at 6.6 µm. In the present study we focus on a 12 µmdiameter suspended-core fiber; this large core fiber exhibits a single ZDW at 4.6 µm, it also avoids significant light coupling issues in the mid-IR region. One can note that in the case of the simple step-index fiber profile studied in Ref. [20], the effect of the core dimension on the dispersion profile is quite different and the ZDW, when it exists, remains located beyond 5.5 µm whatever the core diameter. The calculated evolution of the effective mode area of our suspended-core fiber is shown in Fig. 4b. In particular, we denote a rather constant confinement of the fundamental guided mode over the large mid-IR wavelength region under study. In the following, we consider the nonlinear coefficient n 2 = 8 × 10 -18 m 2 W -1 corresponding to an average of values reported in Refs. [31][43]; thus the effective nonlinear parameter of our suspended core fiber at 8 µm is found to be γ ~ 0.07 W -1 m -1 .

Experimental configuration for supercontinuum generation

The experimental set-up is shown in Fig. 5. Our suspendedcore fiber is pumped with a high power laser source. The tunable central wavelength (from 3.5 to 11 µm) is delivered by a non-collinear optical parametric amplifier (NOPA) followed by a difference frequency generation (DFG) module. The NOPA is a TOPAS-Prime model, pumped by a chirped pulse amplified (CPA) Ti:sapphire system which delivers pulses centred at 790 nm with 3.2 mJ energy (100fs pulse duration) at 1 kHz repetition rate. The DFG module delivers 65 fs pulses with a maximum average power of 3.5 mW at 8.6 µm (54 MW peak power). The infrared beam is focused in our 40 mm-long chalcogenide fiber with a ZnSe objective (f = 6 mm) to excite the fundamental mode. At the fiber output, the SC light is then collimated with a ZnSe objective (f = 12 mm). The collimated beam is then focused into the slit of a monochromator with a gold coated parabolic mirror (f = 88.9 mm). Our Newport Cornerstone TM130 1/8 m motorized monochromator includes a 150 l/mm diffraction grating used for the 3-13 µm range with a peak efficiency of 85% at the 5 µm blaze wavelength. The output light is then focused onto a liquid nitrogen cooled HgCdTe detector (Hamamatsu P3412-02) with a gold coated mirror (f = 50.8mm). We also made use of several long-pass filters (1.5, 2, 3, 5, 7 and 9 µm) inserted before the detector to avoid higher diffraction orders. Finally the full SC is reconstructed by gathering the filtered spectra.

Results and discussion

In a preliminary study [20], we studied SC generation in a similar suspended-core fiber (with a slightly different composition) with pumping wavelengths (3.5 and 5.5 µm) located in the normal dispersion regime (D < 0) or close to the fiber ZDW. Such wavelengths were first chosen to benefit from the high power delivered by the pump laser. However, we were not able to extend the SC bandwidth as far as in the step-index fiber configuration [20]. Here we present new experimental results, namely using pumping wavelengths centred at 7.5 and 8.6 µm. The corresponding SC spectra obtained in a 4-cm-long segment of 12 µmsuspended-core Ge 20 Se 70 Te 10 fiber are shown in Figs. 6a and6c. We confirm that pumping the fiber further in the mid-IR allows now to push the mid-IR SC edge towards ~13 µm. To date, this corresponds to the largest mid-IR extension of SC generation in a suspended-core optical fiber. We still not reach the fondamental long-wavelength transmission edge at 16 µm (see Fig. 2) since we are limited by the strong absorption band centred at 12.7 µm induced by a combination of different bonds involving oxygen. Moreover, our detection setup is not optimized for wavelengths above 12-13 μm (degraded efficiency of the grating and sensitivity of the MCT detector). In our experiments, the available average powers before the coupling objective are 5 and 3.5 mW at 7.5 and 8.6 µm wavelength, respectively. The low average power and relatively degraded laser beam quality at such pumping wavelengths imply a weak coupling efficiency into our fiber, thus preventing the measurement of output average powers.

To further investigate SC generation (peak power involved, SC bandwidth limitations, etc.) we performed numerical simulations of the nonlinear pulse propagation in our fiber segment. Our modelling is based on the generalized nonlinear Schrödinger equation that includes the full dispersion curve (as shown in Fig. 4), both instantaneous Kerr and delayed Raman nonlinear responses, and the dispersion of nonlinearity [44]. The red areas denote the upper detection limit of our experimental setup, while gray areas are related with the location of the main absorption peak measured in Fig. 2.

The nonlinear fiber parameter was calculated above by means of the simulated effective mode area. For the Raman response function, we used an intermediate-broadening model using convolutions of Lorentzians and Gaussians adapted from spontaneous Raman scattering spectra given in Refs. [31][45]. Our simulations also include the losses of our GST glass measured in Fig. 2. We also considered 80-fs pump pulses, for both pumping wavelengths, which were deduced from the initial pump spectra shown in Figs. 6a and6c. Numerical SC spectra obtained are depicted in Figs. 6b and6d for both pumping configurations. We carrefully took into account the same power ratio between the two pumping wavelengths as in experiments, in particular we used 270 kW peak power at 7.5 µm and 190 kW peak power at 8.6 µm. We observe that our simulations reproduce qualitatively the main SC bandwidth in both cases. We also confirm the twooctave spanning mid-IR light source (from 3.25 to 13 µm) with ~15 nJ femtosecond pulses at 8.6 µm. This main part of the SC is only related to the anomalous dispersion of our fiber, it then results from self-phase modulation and soliton dynamics. The main limitation of SC extension on the long wavelength edge remains the strong absorption band measured at 12.7 µm, however, we reveal that spectral broadening can go through this absorption and reach the multiphonon absorption in the mid-IR near 16 µm. This extra-broadening was not observed experimentally due to the limited bandwidth of our detection system. In additon, one can note that simulated SC spectra are also broader on the short-wavelength edge (up to 2 µm) in contrast to our measurements. As the SC light collection was optimized to highlight the main SC part and the long-wavelength edge, and as possible higher OH absorption occurs around 2.9 µm, all this could explain a lower efficiency to analyze the dispersive waves generated in the normal dispersion regime (i.e., the spectral band on the short wavelength edge). Such issues were also discussed in Ref. [18]. Finally, both recorded SC spectra strongly drop at ~4.5 µm, thus evidencing the position of the ZDW in agreement with numerical predictions.

Conclusion

In conclusion, we demonstrated the fabrication and the application of a short segment of Ge-Se-Te glass-based suspended core fiber for mid-infrared supercontinuum generation that extends from 3.25 to 13 µm. Again this work clearly confirms the potential of these arsenic-and antimonyfree chalcogenide optical fibers for nonlinear optics in the mid-IR. To date, we have reported the largest mid-IR extension of SC generation in a suspended-core optical fiber. Moreover, our numerical simulations highlight the potential of such fibers to reach the multiphonon absorption edge of Se-rich GST glasses.This still requires optimization of our mid-IR detection system and an improved purification process as the one we used in [19]. In comparison to similar step-index GST fibers [20], the development of fibers with a smaller suspended core (higher nonlinearity and blue-shifted ZDW) could be suitable for compact SC sources when considering new femtosecond fiber lasers operating around 3-3.5 µm [46].

Figure 2 .

 2 Figure 2. Multi-phonon absorption is measured around 17 µm.

Figure 1 .

 1 Figure 1. Illustration of (a) distillation under static vacuum and (b) refining step in the rocking furnace.

Figure 2 .

 2 Figure 2. Absorption coefficient of GST glass. The bottom graph is a zoom on the vertical axis for a better observation of absorption bands. Our slight purification process allows to decrease extrinsic absorption bands but not to entirely remove them: especially OH, H 2 O, Se-H, Ge-O, Te-O and Se-O bonds vibration bands are still visible [30][31][35]. However, absorption coefficient due to impurities between 2 and 10 µm are smaller than 0.3 cm -1 which is low enough for nonlinear experiments on short fiber length, see part 4. The absorption due to Se-O, Ge-O and Te-O bonds related with the band centred at 12.7 µm corresponds to an absorption coefficient of 7.5 cm -1 . This value is lower than recent work considering a non-purified glass [30].Table 2 lists the distinct absorption bands measured in our glass.

Figure 3 .

 3 Figure 3. Pictures of the drilled preform (upper left) and the transverse profile (upper right) of the fiber preform. Lower part,

Figure 4 .

 4 Figure 4. (a) Calculated wavelength-dependent dispersion D of the fundamental guided mode in suspended-core fibers with distinct core diameters ϕ = 4, 8, and 12 μm. The solid black line indicates the bulk dispersion of the glass (Ge 20 Se 70 Te 10 ). (b) Calculated wavelength-dependent effective mode area of the fundamental guided mode in the 12 µm-suspended-core fiber.

Figure 4 .

 4 Figure 4. Experimental setup used for supercontinuum generation

Figure 6 .

 6 Figure 6. Experimental (left column) and numerical (right column) SC generation results (solid lines). Initial pump pulse spectra are indicated with dashed lines. The corresponding pumping wavelength is 7.5 µm in subplots (a) and (b), and 8.6 µm in subplots (c) and (d). The red areas denote the upper detection limit of our experimental setup, while gray areas are related with the location of the main absorption peak measured in Fig. 2.

Table 1 .

 1 Recent survey on mid-IR SC generation in various configurations of chalcogenide fibers

	Fiber Composition (core/clad)	Spectral Coverage (µm)	Pump Wavelength (µm)	Year
	As-S suspended-core or tapered fibers
	As 2 S 3 /As 2 Se 3	1.4-7.2	3.25	2017 [11]
	As 2 S 3	1.0-4.8	3.0	2018 [12]
	As 2 S 3 /Ge 20 As 20 Se 15 Te 45	2.1-7.0	4.5	2018 [13]
	As 2 S 3 /Ge 15 Sb 10 Se 75	1.5-9.5	3.0	2018 [14]
	As 2 S 3 /As 38 S 62	1.5-8.3	3.75	2018 [15]
		As-Se suspended-core fibers
	AsSe 2 /As 2 S 5	1.3-5.4	3.39	2014 [16]
	As 38 Se 62	1.7-7.5	4.4	2015 [17]
	As 2 Se 3	2.8-7.1	ZBLAN Cascade	2016 [18]
	Other glass-system suspended-core fibers
	Ge 20 As 20 Se 15 Te 45 /Ge 20 Sb 5 Se 75	1.7-11.3	5	2018 [19]
	Ge 20 Se 60 Te 20	2.0-11.5	5.5	2019 [20]
		As-Se step-index fibers
	As 40 Se 60 /Ge 10 As 23.4 Se 66.6	1.4-13.3	6.3	2014 [10]
	As 2 Se 3 /AsSe 2	2.0-15.1	9.8	2016 [21]
	As 2 Se 3 /As 2 S 3	2.0-12.0	3	2017 [22]
		Other glass-

system step-index fibers

  Ge 20 As 20 Se 15 Te 45 /Ge 20 As 20 Se 17 Te 43 Ge 20 As 20 Se 15 Te 45 /Ge 19 As 16 Se 25 Te 40 /Ge 20 As

		1.5-14.0	4.5	2016 [23]
	(Ge 10 Te 43 ) 90 -AgI 10 /(Ge 10 Te 40 ) 90 -AgI 10 /Ge 10 Sb 10 Se 80	2.0-16.0	7	2017 [1]

20

Se 17

Te 43

Table 2

 2

lists the distinct absorption bands measured in our glass.

Table 2 .

 2 Intensity and position of absorption bands compared to the literature

	Chemical	 (µm)	 (cm -1 )	 (dB/m)	 (dB/m/ppm)	Content (ppm)	Comparison with
	species						literature
	Se-H	4.25	0.025	10.7	1 [33][36]	10.65 ppm(at)	[20][33][34][37]
	OH	2.9	0.141	60.7	10 [38]	6.11 ppm(wt)	[20][32][34][37]
	Ge-O	7.9	0.232	100.7	-	-	[20][35]
	Te-O	15	0.525	228.0	-	-	[20]
	H 2 O	6.3	0.065	28.2	-	-	[20][37]

Table 3 .

 3 Sellmeier coefficients for GeSe 4 , GeTe 4 and Ge 20 Se 70 Te 10

	Glass	A 1	A 2	λ 1 (µm)	λ 2 (µm)
	composition
	GeSe 4	4.813	0.126	0.278	19.738
	GeTe 4	9.943	0.318	0.522	18.298
	Ge 20 Se 70 Te 10	5.225	0.202	0.298	21.606

Acknowledgements

We acknowledge the financial support of the French Investissements d'Avenir program (ISITE-BFC ANR-15-IDEX-0003, SCUVIRA project), the Council of the French Région Bourgogne Franche-Comté, and the European program FEDER, through the grant "Jeune Chercheur Entrepreneur" scheme.


Evgenia Kartsaki

Bruno Cessac

Gerrit Hilgen

Evelyne Sernagor

Probing retinal function with a multi-layered simulator

☞

  of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. Cell ü A set of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. ü A set of parameters that constrain the cell's evolution : e.g. conductance, reversal potential, membrane capacitance, characteristic time of a channel's activity Cell ü A set of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. ü A set of parameters that constrain the cell's evolution : e.g. conductance, reversal potential, membrane capacitance, characteristic time of a channel's activity ü A function controlling the cell's evolution with a differential equation Cell ü A set of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. ü A set of parameters that constrain the cell's evolution : e.g. conductance, reversal potential, membrane capacitance, characteristic time of a channel's activity ü A function controlling the cell's evolution ü Isyn : A synaptic input corresponding to synaptic connections with other cells Cell ü A set of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. ü A set of parameters that constrain the cell's evolution : e.g. conductance, reversal potential, membrane capacitance, characteristic time of a channel's activity ü A function controlling the cell's evolution ü Isyn : A synaptic input corresponding to synaptic connections with other cells ü Iext : An external input corresponding either to a visual input or the electric current provided by an electrode Visual or Electrical Input Basic building blocks -synapse Synapse ü Chemical or electrical (gap junction) ü A set of variables that evolve in time : e.g. conductance ü A set of parameters : e.g. synaptic weight Design a local circuit with specific connectivity patterns ☞ Deploy it to the whole retina Source: From Gabriel Luna, Steve Fisher and Geoff Lewis, retinal cell biology group at UC Santa Barbara Neuroscience Research Institute Connectivity graph 𝐾 𝑥, 𝑦, 𝑡 = 𝐾 ' 𝑥, 𝑦 𝐾 ( (𝑡) +(𝑆 * K)(x,y,t) 𝑦, 𝑡 = 𝐾 ' 𝑥, 𝑦 𝐾 ( (𝑡) +(𝑆 * K)(x,y,t) 𝐼 345 = -𝑔 345 𝑉 9:3; -𝐸 345 𝐼 =>? = -𝑔 =>? 𝑉 9:3; -𝑉 9@A 𝑦, 𝑡 = 𝐾 ' 𝑥, 𝑦 𝐾 ( (𝑡) +(𝑆 * K)(x,y,t) 𝐼 345 = -𝑔 345 𝑉 9:3; -𝐸 345 𝐼 =>? = -𝑔 =>? 𝑉 9:3; -𝑉 9@A D 𝑉 -𝑉 D + 𝐼𝑠𝑦𝑛 + 𝐼𝑒𝑥𝑡 + 𝐼 IJK 𝐾 𝑥, 𝑦, 𝑡 = 𝐾 ' 𝑥, 𝑦 𝐾 ( (𝑡) +(𝑆 * K)(x,y,t) 𝐼 345 = -𝑔 345 𝑉 9:3; -𝐸 345 𝐼 =>? = -𝑔 =>? 𝑉 9:3; -𝑉 9D 𝑉 -𝑉 D + 𝐼𝑠𝑦𝑛 + 𝐼𝑒𝑥𝑡 + 𝐼 IJK 24 25

Acknowledgements

Prof. Evelyne Sernagor Dr. Gerrit Hilgen Dr. Bruno Cessac AMDT team

General structure

Cell type 1

Cell type 2

Visual or Electrical Input

Inter-type connectivity


Evgenia Kartsaki

Bruno Cessac

Gerrit Hilgen

Evelyne Sernagor

Probing retinal function with a multi-layered simulator

Retinal neurons inventory

Source: Masland et al., 2001 Leverhulme funded project: "A novel approach to functional classification of retinal ganglion cells" ☞ Understand their role in population encoding of complex visual scenes ☞ Achieve full characterisation of retinal ganglion cells sub-classes based on:

• pharmacogenetics • large-scale retinal electrophysiology • immunochemistry Coombs et al., 2006Baden, Berens, Franke et al., 2016Martersteck, Hirokawa, Evarts et al., 2017 o

oo☞

  Novel and interdisciplinary o Combines shared gene expression with physiology and parts of anatomy o Identify RGCs sharing a common gene expression by combining pharmacogenetics, large-scale high-density electrophysiology and immunohistochemistry Aims o Identify RGCs sharing a common gene expression by combining pharmacogenetics, large-scale high-density electrophysiology and immunohistochemistry o Isolate subgroups from that RGC pool based on responses to different stimuli Aims o Identify RGCs sharing a common gene expression by combining pharmacogenetics, large-scale high-density electrophysiology and immunohistochemistry o Isolate subgroups from that RGC pool based on responses to different stimuli Unravel how RGCs in these subgroups respond and interact to basic and complex visual scenes Aims o Identify RGCs sharing a common gene expression by combining pharmacogenetics, large-scale high-density electrophysiology and immunohistochemistry o Isolate subgroups from that RGC pool based on responses to different stimuli Unravel how RGCs in these subgroups respond and interact to basic and complex visual scenes of body, head, or eye of the observer motion of objects within the scene Global motion Object region: Excitatory input to the ganglion cell (G) from many fast bipolar cells (B). Background region: Inhibitory input to the ganglion cell (G) from amacrine cells (A). Retinal circuit related to motion Inspired from Baccus et al., 2008 DREADDs on Ganglion cells (RGCs) -If excitatory, the neuronal activity will increase -If inhibitory, the neuronal activity will decrease CNO effect on the circuit Inspired from Baccus et al., 2008 DREADDs on Amacrine cells (ACs): -If excitatory, the neuronal activity of G will decrease -If inhibitory, the neuronal activity of G will increase CNO effect on the circuit Inspired from Baccus et al.makes the interpretation of functional results from DREADD-expressing RGCs more challenging. of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. Cell ü A set of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. ü A set of parameters that constrain the cell's evolution : e.g. conductance, reversal potential, membrane capacitance, characteristic time of a channel's activity Cell ü A set of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. ü A set of parameters that constrain the cell's evolution : e.g. conductance, reversal potential, membrane capacitance, characteristic time of a channel's activity ü A function controlling the cell's evolution with a differential equation Cell ü A set of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. ü A set of parameters that constrain the cell's evolution : e.g. conductance, reversal potential, membrane capacitance, characteristic time of a channel's activity ü A function controlling the cell's evolution ü Isyn : A synaptic input corresponding to synaptic connections with other cells Cell ü A set of variables evolving in time and characterizing the cell's evolution : e.g. membrane potential, probability that a ionic channel of a given type is open etc. ü A set of parameters that constrain the cell's evolution : e.g. conductance, reversal potential, membrane capacitance, characteristic time of a channel's activity ü A function controlling the cell's evolution ü Isyn : A synaptic input corresponding to synaptic connections with other cells ü Iext : An external input corresponding either to a visual input or the electric current provided by an electrode Visual or Electrical Input Basic building blocks -synapse Synapse ü Chemical or electrical (gap junction) ü A set of variables that evolve in time : e.g. conductance ü A set of parameters : e.g. synaptic weight Design a local circuit with specific connectivity patterns ☞ Deploy it to the whole retina Source: From Gabriel Luna, Steve Fisher and Geoff Lewis, retinal cell biology group at UC Santa Barbara Neuroscience Research Institute

Source: ETH Zürich & RavaAzeredo da Silveira et al., 2011 Macular


Estimation of reference curves for fetal weight

SANDIE FERRIGNO a , MYRIAM MAUMY-BERTRAND b a Inria team BIGS Nancy -Grand-Est, Lab IECL, Université de Lorraine, France b Lab IRMA, Université de Strasbourg, France

OBJECTIVES

Reference or standard curves are required in many biomedical problems. Values which lie outside the limits of these reference curves may indicated the presence of disorder.

Data are from the French EDEN mother-child cohort [START_REF] Heude | Cohort profile: The eden motherchild cohort on the prenatal and early postnatal determinants of child health and development[END_REF]. We are studying fetal weight (in grammes) that depends on the gestional age (GA) in the second and the third trimesters of mother's pregnancy.

Some classical parametric and semi-parametric methods as polynomial [START_REF] Royston | How to construct 'normal ranges' for fetal variables[END_REF][3] and LMS method [4] [START_REF] Cole | Smoothing reference centile curves : the lms method and penalized likelhood[END_REF] or nonparametric methods as kernel estimation [START_REF] Gannoun | References curves based on non-parametric quantile regression[END_REF] are used to construct these curves. However, some of them requires strong assumptions.

EDEN COHORT [1]

EDEN is a mother-child cohort study investigating the prenatal and early postnatal determinants of child health and development.

2002 pregnant women were recruited before 24 weeks of amenorrhoea in two maternity clinics from middle-sized French cities (Nancy and Poitiers). From May 2003 to September 2006, 1899 newborns were then included.

Detailed information on phenotypes and exposures were collected through questionnaires from preg-nancy until 8 years and through clinical examinations of the mother and the child.

The main outcomes of interest are fetal (via ultrasound) and postnatal growth, adiposity development, respiratory health, atopy, behaviour and bone, cognitive and motor development.

The main early-life determinants studied are maternal nutrition, diet, eating behaviour, environmental pollutants and socioeconomic and psychoemotional factors.

We acknowledge the commitment of the EDEN mother-child cohort study group and more particularly Barbara Heude which allowed us to use the data from the cohort. An α% age-related reference interval is the range bounded by two centil curves which encompasses α% of the data at each age.

DATA AND REFERENCE INTERVAL

The estimate of the reference curves will be carried out from 30 GA because before we do not have enough data.

POLYNOMIAL METHOD [2][3]

Polynomial regression is one of the most common parametric approach for modelling growth data espacially during the prenatal period. It is based on the asumption that at each t (GA) the measurement of weight has a normal distribution and that the mean and standard deviation vary smoothly with GA. The C α (.) centile curves is estimated using

C α (t) = m(t) + Z α σ(t), (1)
where Z α is the normal equivalent deviate of size α corresponding to a particular centile and m(.) and σ(.) are the mean and standard deviation, respectively, at the required GA. In this study, these two functions are respectively estimated by polynomials of order three and one.

LMS METHOD [4][5]

The LMS methods assumes that at each age the weight follows a normal distribution after a Box-Cox transformation. Therefore, data can be summarised by three GA-dependent functions that will be estimated and smoothed to obtain the location parameter curve M (t), the coefficient of variation curve S(t) and the value of the power needed to normalise the data at each age, L(t). These three curves are fitted using cubic splines and by maximising the penalised likelihood. They together allow any centiles to be calculated

C α (t) = M (t) (1 + L(t) * S(t) * Z α ) 1/L(t) , L(t) = 0 C α (t) = M (t) exp(S(t) * Z α )), L(t) = 0

NONPARAMETRIC METHOD [6]

Nonparametric method is used in equation ( 1) to estimate the GA-specific mean function m(.) and the standard variance function σ 2 (.). These estimates are obtained using Nadaraya-Watson kernel estimation. mn (.) and σ2 n (.) are calculated from n realizations {(X i , Y i ), i = 1, . . . , n} from the random variables X (GA) and Y (weight) by

mn (t) = n i=1 K X i -t h Y i n i=1 K X i -t h , σ2 n (t) = n i=1 K X i -t h Y 2 i n i=1 K X i -t h -m2 n (t),
where h and K(.) are respectively a bandwidth parameter and a kernel function.

Figure 1 :

 1 Figure 1: Scatter plot of fetal weight vs. GA.

Figure 2 :

 2 Figure 2: 3rd, 10th, 50th, 90th, 97th centile curves.

Figure 3 :

 3 Figure 3: 3rd, 10th, 50th, 90th, 97th centile curves.

Figure 4 :

 4 Figure 4: 3rd, 10th, 50th, 90th, 97th centile curves.


Mingxiao Ma

Abdelkader Lahmadi

Isabelle Chrisment

Demonstration of Synchronization Attacks on Distributed and Cooperative Control in Microgrids

Demonstration of Synchronization Attacks on Distributed and Cooperative Control in Microgrids

Mingxiao Ma, Abdelkader Lahmadi, Isabelle Chrisment Université de Lorraine, CNRS, Inria, Loria, F-54000 Nancy, France Email: {firstname.lastname}@loria.fr we have only studied its impact through theoretical analysis and Matlab based simulations. In this paper, we demonstrate the practical feasibility and the impact of the MaR attack using a hardware platform modeled after a simplified microgrid.

II. ATTACK MODEL AND EXPERIMENTAL SETUP In our previous work [START_REF] Ma | On the Impact of Synchronization Attacks on Distributed and Cooperative Control in Microgrid Systems[END_REF], we proposed a distributed and cooperative control structure for a microgrid system. Specifically, we use quadratic droop control as the primary control layer and cooperative control as the secondary control layer. Simulation results have proved that this design is more realistic and reliable than traditional control schemes as proposed in [START_REF] Teixeira | Voltage control for interconnected microgrids under adversarial actions[END_REF] and [START_REF] Ma | Voltage control in distributed generation under measurement falsification attacks[END_REF]. In this work, and for experimental reasons, we simplify the microgrid model presented in [START_REF] Ma | On the Impact of Synchronization Attacks on Distributed and Cooperative Control in Microgrid Systems[END_REF] and build a hardware platform to validate attacks on communication network used by the control layers. As shown in Figure 1, the simplified microgrid consists of N DG units and a sparse communication network. In this microgrid platform, we characterize the network topology by a directed graph (digraph) G = (V, E, A G ) with a nonempty finite set of N nodes V = {v 1 , v 2 , . . . , v N }, a set of edges E ⊆ V × V, and the associated adjacency matrix A G . Specifically, DG units and their corresponding communication links are considered as the nodes and edges of the communication graph. As shown in Figure 1, only one leader node (DG Unit 1 ) gets access to the the reference V * and all other nodes will synchronize their corresponding reference values with the leader node while communicating with their neighbouring nodes.

For each DG Unit i , the required reference value V * i (t) and voltage measurement V s i (t) are transmitted through the Abstract-Microgrid platforms are experiencing an increasing integration of distributed and cooperative control coupled with communication links to guarantee safe and reliable operations of their distributed power generators (DG). However, they are subject to several threats and attacks because of the widely use of communication networks in their control layers. In this demo paper, we demonstrate the practical feasibility and the impact of a novel attack, named measurement as reference attack (MaR), on the distributed control mechanism of a microgrid, where an attacker replaces the reference values with measurements during its synchronization operations. We build in this work a hardware platform modelled after a simplified m icrogrid w ith D G units based on Raspberry Pi and Arduino boards, DC motors for power generation and light bulbs as electrical loads. We study two variations of the MaR attack using man-in-the-middle (MITM) and malware techniques to demonstrate and validate its impact on the microgrid synchronization and its voltage stability which affect end user's electrical devices.

I. INTRODUCTION

Modern power distribution systems relying on microgrids will have an increasingly high integration of distributed power generators (DG). A microgrid is a cluster of DG units that coordinates the inter-connected DGs by interfacing them through power electronic devices such as the voltage sourced inverter (VSI) [START_REF] Katiraei | Micro-grid autonomous operation during and subsequent to islanding process[END_REF]. To ensure safe and reliable operations of microgrids, a hierarchical control structure using supervisory control and data acquisition (SCADA) systems is applied [START_REF] Khodayar | Integration of high reliability distribution system in microgrid operation[END_REF] to collect data from remote facilities and send back control instructions to those facilities. In these systems, communication networks play an important role to transmit data between their distributed controllers and generators [START_REF] Isozaki | On detection of cyber attacks against voltage control in distribution power grids[END_REF].

However, using communication networks at the different control layers of microgrids introduces challenges regarding their security where they become targets of cyber attacks. The well-known power blackout caused by the malware BlackEnergy in Ukraine during 2015 has proved that cyber attacks could cause a major blackout. Thus, it is of great importance to study potential vulnerabilities and attacks of these systems and design mitigation or detection schemes against the highrisk threats [START_REF] Teixeira | Voltage control for interconnected microgrids under adversarial actions[END_REF]. In a previous work [START_REF] Ma | On the Impact of Synchronization Attacks on Distributed and Cooperative Control in Microgrid Systems[END_REF], we have designed a novel attack, named Measurement as Reference attack (MaR), where an attacker is able to affect the voltage stability of the microgrid by replacing the reference voltage values with the measurement value while they are exchanged between distributed and cooperative controllers. In the previous work, communications links. As shown in Figure 2, each DG is composed of of Raspberry Pi and Arduino Mega 2560 board acting as the primary controller, two motors A and B acting as a voltage generator, and one light bud acting as an electrical load. The Raspberry Pi is the main control center of the DG unit which sends commands to the Arduino board through USB port to control the voltage generator by using a PI algorithm. The generated voltage will power the light bulb, and its fluctuation affects the level of brightness. For the system described in Figure 1 and the experimental setup depicted in Figure 2, we implement our MaR (Measurement as Reference) attack that affects the voltage stability and reference signal synchronization as we studied it theoretically in [START_REF] Ma | On the Impact of Synchronization Attacks on Distributed and Cooperative Control in Microgrid Systems[END_REF]. We assume that the attacker has knowledge about the distributed and cooperative control structure where it targets the communication links between DG units or the controller station to maliciously replaces the reference signal with the measurement of the previous node. This attack is "naturally stealthy" because usually the reference value and the voltage measurement are close enough and it would be difficult for traditional threshold based detectors to identify malicious data replacement without causing false alarms.

Without loss of generality, the MaR attack targets the communication link between DG units i and i + 1, where the attacker accesses and manipulates the exchanged data by replacing the reference signal V * i+1 (t) for DG unit i + 1 with the voltage measurement V s i (t) of DG unit i , witten as:

V * i+1 (t) = V s i (t).
III. MAR ATTACK DEMONSTRATION In this section, we demonstrate two types of attack techniques, man-in-the-middle attack (MITM) and malware, to implement the MaR attack as described in Section II.

A. MaR with a man-in-the-middle attack

In this scenario, we demonstrate the MaR attack when the attacker uses a MITM technique to introduce itself between two communicating DGs neighbours. We assume, in this scenario, that the communications between DGs are not encrypted and using the TCP protocol. As shown in Figure 3, only neighbouring nodes exchange data through the communication link between them. The attacker can introduce itself into a conversation between two DG units and acts as relay/proxy, impersonates both nodes and gets access to the reference and measurements values that the two nodes are sending to each other.

B. MaR with a malware based attack

In the second scenario, we demonstrate the MaR attack by using a malware technique to infect a DG, as shown in Figure 3. In this demonstration, we craft a malware that installs itself on the network hosting the DGs nodes. It scans the network looking for open SSH or Telnet ports on the DGs devices. When it founds an open port on one of them, it uses a dictionary of credentials to brute-force the access to the device. After success, it installs itself and replaces the control program with a new version that replaces the reference value with the measurement. In this scenario, we mainly show the feasibility of the MaR attack by infecting a DG when the communications between DGs are encrypted.

  hal.archives-ouvertes.fr/hal-02389307 Submitted on 2 Dec 2019

Fig. 1 .

 1 Fig. 1. A microgrid control platform consisting of primary and secondary controllers. Neighbouring DG units can exchange data through sparse communication networks, e.g., transmitting reference signals and voltage measurements (denoted by the superscript s).

Fig. 2 .

 2 Fig. 2. Each DG unit is composed of one Raspberry Pi, one Arduino board, two motors and one light bulb to demonstrate the effect of voltage instability on the level of brightness.

Fig. 3 .

 3 Fig. 3. MaR attack implementations using either MITM or malware techniques.

ACKNOWLEDGMENT

This work has been funded by the French Government under grant FUI 23 PACLIDO (Protocoles et Algorithmes Cryptographiques Légers pour l'Internet Des Objets).