Schertzer D Hmco

Pseudo-radar algorithms with two extremely wet months of disdrometer data in the Paris area

Keywords: rainfall, disdrometer, radar algorithm, multifractal

Disdrometer data collected during the two extremely wet months of May and June 2016 at the Ecole des Ponts ParisTech are used to get insights on radar algorithms. The rain rate and pseudo-radar quantities (horizontal and vertical reflectivity, specific differential phase shift) are all estimated over several durations with the help of drop size distributions DSD) collected at 30 s time steps. The pseudo-radar quantities are defined with simplifying hypotheses, in particular on the DSD homogeneity. First it appears that the parameters of the standard radar relations

for these pseudo-radar quantities exhibit strong variability between events and even within an event. Second an innovative methodology that relies on checking the ability of a given algorithm to reproduce the good scale invariant multifractal behaviour (on scales 30 s -few h) observed on rainfall time series is implemented. In this framework, the classical hybrid model ( R Z h -for low rain rates and dp K Rfor great ones) performs best, as well as the local estimates of the radar relations' parameters. However, we emphasise that due to the hypotheses on which they rely these observations cannot be straightforwardly extended to real radar quantities.

1) Introduction

Disdrometers are rainfall point measurement devices that give access not only to rain rates but also to information about the size and velocity of drops falling through the sampling section. They are commonly implemented for research purposes and their operational use is spreading. Weather radars are basically the only devices providing volumetric rainfall data almost instantaneously, but in a less direct measurement. Indeed a radar transmits a wave and actually measures the power of the wave backscattered by the hydrometeors of the atmosphere, and not the rain rate which is the quantity hydrometeorologists are interested in. With the strong simplifying hypothesis that the drop size distribution (DSD) is homogeneous in the radar scanned volume and corresponds to the one measured with the help of disdrometer, it is possible to define pseudo radar quantities based on the pointwise disdrometer measurement. This furthermore assumes incoherent scattering, i.e. there is no clustering of the drops inside radar pulse volumes, whose radial scale is usually in the range of 100 m to 1 km. In other words that drop centres are homogenously distributed. It has been shown that this hypothesis may lead to strong statistical biases [START_REF] Lovejoy | Scalar multifractal radar observers problem[END_REF], Schertzer et al., 2012). These hypotheses enable to partially inverse the retrieval algorithms from radar backscattering to rain rate. This approach was implemented by numerous authors [START_REF] Gires | Multifractal comparison of the outputs of two optical disdrometers[END_REF]; Jaffrain and Berne 2012b; [START_REF] Leinonen | A Climatology of Disdrometer Measurements of Rainfall in Finland over Five Years with Implications for Global Radar Observations[END_REF][START_REF] Ryzhkov | Rainfall Estimation with a Polarimetric Prototype of WSR-88D[END_REF][START_REF] Verrier | Theoretical and empirical scale dependency of Z-R relationships: Evidence, impacts, and correction[END_REF]. Others relied on synthetic DSD functions and ran various simulations using probability distributions for the key parameters characterizing the DSD [START_REF] Anagnostou | Performance evaluation of high-resolution rainfall estimation by X-band dualpolarization radar for flash flood applications in mountainous basins[END_REF][START_REF] Illingworth | The Need to Represent Raindrop Size Spectra as Normalized Gamma Distributions for the Interpretation of Polarization Radar Observations[END_REF][START_REF] Ryzhkov | Rainfall Estimation with a Polarimetric Prototype of WSR-88D[END_REF].

In this paper we use data collected in May and June 2016 by the TARANIS observatory (exTreme and multi-scAle RAiNdrop parIS observatory, [START_REF] Gires | Multifractal comparison of the outputs of two optical disdrometers[END_REF]) of the Fresnel Platform of Ecole des Ponts ParisTech (https://hmco.enpc.fr/Page/Fresnel-Platform/en). This period was extremely rainy over this area, and resulted in severe floods of the Seine River in Paris and its tributaries. The peak flow in Paris was reached on 7 June 2016 with levels not observed since 1982 on the Seine River. Some upstream tributaries of the Paris area, especially the Loing, reached levels greater than the ones of the famous 1910 flood (http://www.developpement-durable.gouv.fr/Point-sur-les-inondations-6-juin.html). This data is used to get insights on various radar algorithms.

The data is presented in section 2, with a highlight on two events exhibiting different features; an almost two day uninterrupted period of light rain (29-31 May 2016), and a short heavy one (17 June 2016). A multifractal analysis of the rainfall data is used to quantify the variability across scales of the measured time series (see Schertzer and Lovejoy 2011 for a recent review). The radar relations are investigated in section 3, after having defined the following pseudo-radar quantities: the horizontal reflectivity Z h , the specific phase shift K dp and the differential reflectivity Z dr . The latter two depend on the oblate shape of drops. The parameters of the power law relations between these quantities and the rain rate are assessed on average over the whole period and also for each event. The variability during an event is studied as well. Finally in section 4, the consequences of the previous findings on rainfall retrieval with radars are analysed. More precisely reconstructed radar observations of rainfall time series (using the relations described in section 3) are compared with actual rainfall time series. An innovative comparison based on multifractal properties is carried out as well as the computation of standard scores commonly used by hydro-meteorologists.

2) Rainfall data and preliminary analysis

2.1) Data collection and computation

The data analysed in this paper was collected between 1 May 2016 and 31 June 2016 by the TARANIS observatory's optical disdrometers which are collocated on the roof of the Ecole des Ponts ParisTech building (East of Paris, France). Two OTT Parsivel 2 disdrometers (Loffler-Mang and Joss 2000, for an initial version; [START_REF] Battaglia | PARSIVEL Snow Observations: A Critical Assessment[END_REF][START_REF] Ott | Operating instructions[END_REF] are the first type of disdrometer used and denoted Pars#1 and Pars#2 hereafter. The two Parsivel 2 are installed perpendicularly. Such optical disdrometer transmits a laser sheet and computes the size and velocity of drops falling through the sampling area by assessing the portion of light occluded before it reaches the receiver. The second disdrometer type is a Campbell Scientific PWS100 [START_REF] Ellis | New laser technology to determine present weather parameters[END_REF], Campbell Scientific Ltd 2012). In this case, it is the light refracted by drops that is analysed. It should be mentioned that a correction of the drops' diameters estimation was implemented to better take into account the oblateness of drops [START_REF] Gires | Multifractal comparison of the outputs of two optical disdrometers[END_REF]. Both devices basically provide a matrix with the number of drop(s) recorded according to classes of equivolumic diameter and terminal fall velocity. From this matrix, it is possible to compute the rain rate R in mm.h -1 (by counting the number of drop(s) that fell during a time interval) and a drop size distribution (DSD) denoted N(D) (in m -3 .mm -1 ). With the aforementioned hypotheses (homogeneous DSD and distribution of the drop centres), the horizontal reflectivity (Z h in mm 6 .m -3 ), the vertical reflectivity (Z v in mm 6 .m -3 ) and the specific differential phase shift (K dp in °.km -1 ) can be estimated as moments of the (local) DSD estimated by the disdrometers:

∫ + - = max min ) ( ) ( 1 1 , ; 2 2 5 4 , D D v h B v h dD D D N m m Z σ π λ (1) [ ] ∫ - = - max min ) ( ) ( ) ( Re 1800 10 3 D D vv hh dp dD D N D S D S K π λ (2)
where D min and D max are the minimum and maximum measured diameters D's (in mm),

v h B , ;
σ (in mm 2 ) the backscattering cross section at horizontal / vertical polarization, Re(S hh/vv ) (in mm) the real part of the forward scattering amplitude at horizontal/vertical polarization, λ the radar wavelength (in mm, for the X-band radar wavelength used in this paper we have λ = 33.3 mm), m the complex refractive index of water (a temperature of 20°C is taken for which m=8.633 +1.289i). The scattering coefficients are obtained with the help of the Python PyTMatrix library (Leinonen 2014), which relies on the T-Matrix code [START_REF] Mishchenko | T-matrix computations of light scattering by nonspherical particles: A review[END_REF]. The following set up is implemented for the computations: (i) an oblate spheroids model for drop shape is taken with an axis ratio -equivolumic diameter relation corresponding to the one implemented in the Parsivel 2 rationale [START_REF] Battaglia | PARSIVEL Snow Observations: A Critical Assessment[END_REF]). In this model the ratio between the vertical and horizontal axis is equal to 1 for D < 1 mm, 0. . This approximation remains valid for D < 3 mm (the linear regression yields a slope equal to 5.95), and a somewhat more complex behaviour is found for greater D. This transition occurs for decreasing drop diameter with decreasing wave length, meaning that it is more relevant to use T-Matrix computations instead of pure Rayleigh approximation at the Xband radar wave length than at the the C-band one. For instance instead of the 3 mm limit found here, the transition occurs at roughly 5 mm at a C-band radar wave length of 53.5 mm).

The term used in the computation of K dp (Eq. 2) exhibits a more complex behaviour with regards to drop equivolumic diameter (Fig. 1.c). Given the simplifying assumption made of spherical drops for D < 1 mm, it theoretically yields a zero anisotropic scattering coefficient on this range and in fact a numerical noise (smaller than 10 -16 software precision used), whereas prolate drops are usually expected [START_REF] Beard | A new model for the equilibrium shape of raindrops[END_REF]. This corresponds to a truncation of K dp , with D min =1 that may introduce a bias for low intensity rainfalls to be discussed below. The small bump for drop ranging from 3 to 4 mm is followed by a power law-like behaviour for greater D's. The differential reflectivity Z dr (no unit) is simply defined as the ratio between the horizontal and the vertical reflectivity:

v h dr Z Z Z = (3)
Observations are carried out with time steps of 30 s. There is no missing data for these two months of observations for the Pasivel 2 s and few ones for the PWS. In some sections of this paper, individual events are studied. An event is defined as a rainy period isolated by at least 15 minutes of dry condition before and after and resulting in more than 0.5 mm of total collected water depth.

The DSD is (at least partly) characterized with the help of the classical total drop concentration N t (m -3 ) and the mass-weighted diameter D m (mm): For better readability of this paper the results will be mainly presented for the device Pars#2, and results for the other two will only be briefly mentioned. Figure 2 Figure 3 presents a scatter plot of the rain rates estimated at 30 s time steps of Pars#1 and PWS vs. Pars#2. It should be mentioned that there is some missing data during heavy rainfall for the PWS device resulting in roughly 7% of points missing. There is no missing data for both Pars#1 and Pars#2. Hence the analysis done in this paragraph is based only on the time steps when all devices are available to allow for a reliable comparison. For these time steps, the cumulative depth is of 181 mm for Pars#1, 191 mm for Pars#2 and 230 mm for PWS. It appears that the estimates for the two Parsivels are similar with a 6% split. The discrepancies with the PWS are larger, with 20% differences in terms of cumulative depths. These differences are great (it is more than 50 % without the correction), and greater than the 10% noticed on a previous study [START_REF] Gires | Multifractal comparison of the outputs of two optical disdrometers[END_REF]). This could be related to the fact that both devices do not work on the same principle, i.e. one estimates drops' size and velocity from occlude light whereas this other uses refracted light.

∫ = max min ) ( D D t dD D N N ( 4

2.3) A focus on two different events

Two events were selected to be studied more precisely throughout the paper. during which great rain rates (Fig. 5.a) were recorded (more than 70 mm.h -1 at 30 second time steps and more than 40 mm.h -1 at 5 min time steps). No significant damage was reported on this area. The DSD exhibits a heavier tail (Fig. 5.c) with D m = 2.08 mm. The number of large drops, i.e. with diameters greater than 3-4 mm, might be overestimated with the Parsivels, as this issue was recently reported by [START_REF] Park | Comparative Evaluation of the OTT PARSIVEL2 Using a Collocated Two-Dimensional Video Disdrometer[END_REF] for intense rainfall. The interpretation that it could be due to the orientation of the Parsivel sampling area with regards to the wind direction is discarded in this case because both perpendicular devices exhibit similar DSD. This possible overestimation could be due to nearly coincident drops intersecting the Parsivel laser causing multiple smaller drops to appear as a larger drop.

2.4) Multifractal analysis of the data

The variability across scales of the measured time series is quantified with the help of the theoretical framework of Universal Multifractals (UM) which has been extensively used to analyse and simulate geophysical fields exhibiting great variability over wide ranges of spatio-temporal scales (see Schertzer and Lovejoy 2011 for a recent review). Only the main elements are presented here, for more details readers are invited to refer to the review.

Let us consider a conservative field ε λ at resolution λ, defined as the ratio between the observation scale l and the outer scale L (λ=L/l). In practice the field measured at the maximum resolution is up-scaled (by averaging consecutive time steps) to obtain it at other resolutions. The first step of a multifractal analysis consists in computing the power spectra. If the field exhibits some scaling feature, then the spectra behaves as a power law with regards to the wave number k:

β - ∝ k k E ) ( (6),
where β is the spectral slope.

If ε λ is a multifractal field, then its statistical moment orders scale with resolution as:

) (q K q λ ε λ ≈ (7),
where K(q) is the moment scaling function that fully characterizes the variability across scales of the field ε λ . The quality of the scaling is assessed with the help of the Trace Moment (TM) analysis. It consists of plotting Eq. 7 in log-log; the slope of the obtained straight line being K(q). In the specific framework of the Universal Multifractals (UM) toward which most multiplicative process converge (this a broad generalization of the central limit theorem, Schertzer and Lovejoy 1987, 1997), K(q) is fully characterized with the help of only two scale invariant parameters C 1 and α. C 1 is the mean intermittency co-dimension and measures the clustering of the (average) intensity at smaller and smaller scales (C 1 =0 for a homogeneous field). α is the multifractality index ( 2 0 ≤ ≤ α

) and measures the clustering variability with regards to intensity level. In this paper, the Double Trace Moment (DTM) technique is used to estimate UM parameters [START_REF] Lavallée | Nonlinear variability and landscape topography: analysis and simulation[END_REF]).

If the studied field φ λ appears to be non-conservative (as it is the case here), then a third parameter H, the non-conservation parameter, should be introduced. It is equal to zero for conservative fields. The TM and DTM analysis should be carried out on the conservative part ε λ of the field written as

H - = λ ε φ λ λ (8)
And H can be estimated with the help of [START_REF] Tessier | Universal Multifractals: theory and observations for rain and clouds[END_REF]):

) 2 ( 2 1 c K H - + = β (9)
Where K c is the moment scaling function of the conservative part.

In practice TM and DTM techniques remain reliable as long as H<0.5. In case of greater H, they should be implemented not on φ λ , but on the underlying conservative field ε λ .

Theoretically, a fractional integration of order H (equivalent to a multiplication by k H in the Fourier space) enables to compute ε λ from φ λ . Here we will use a common approximation [START_REF] Lavallée | Nonlinear variability and landscape topography: analysis and simulation[END_REF]) which consists of taking ε Λ at the maximum resolution simply equal to the renormalized absolute value of the fluctuations of the field i.e.:

) ( ) 1 ( ) ( ) 1 ( i i i i Λ Λ Λ Λ Λ - + - + = φ φ φ φ ε (10)
and then upscaling this field at other resolutions λ.

Analyses are carried out on ensemble average over various samples. Each sample is considered as a realisation of the process and is up-scaled independently. The same methodology as in [START_REF] Gires | Multifractal comparison of the outputs of two optical disdrometers[END_REF] is implemented to select the samples from the available data, i.e. "for each event (i) a sample size is chosen (a power of two, if possible); (ii) the maximum number of samples for this event is computed; (iii) the portion of the event of length equal to the sample size multiplied by the number of samples found in (ii) with the greatest cumulative depth is extracted; (iv) the extracted series is cut into various samples." Sample sizes of power of two is simpler to handle, but implies that some available data are actually not used. More precisely with sample sizes of 64, 128, 256, 512, 1024 and 2048 time steps, the percentage of data actually used of 90, 82, 64, 43, 26 and 11%. Here a size of 64 time steps of 30 s, corresponding to 32 min is used to maximise the amount of data used to estimate UM parameters and get more robust ones. Scaling features will also be discussed with other sample sizes.

UM analysis is mainly discussed for data collected with Pars#2. The quality of the scaling on the available data is very good over the whole range of selected scales (30 s -32 min). See Fig. 6.a for the spectral analysis and 6.b for TM analysis where the quality of the linear regression is very good. Given the great value of the spectral slope (and H, see after), the TM and DTM analyses are implemented not on the field directly but on its. It should be mentioned that when tested with greater sample size, the scaling remains very good (R 2 >0.99) up to 512 time steps (~ 4.3 h). For sample size of 1024 (3 events and 5 samples), there seems to be a break also visible when performing the analysis with 2048 time steps. This limit of 4.3 h is a standard duration for meteorological situations in this area. Given the high quality of the scaling, the DTM analysis is reliable. The UM parameters found are: β = 1.98, C 1 = 0.19, α = 1.64 and H = 0.65. Similar results with same high quality of scaling are found for Pars#1 with β = 2.08, C 1 = 0.19, α = 1.56 and H = 0.70. Due to the missing time steps it was not possible to carry out a similar analysis with the PWS data.

3) Radar relations at various temporal scales

3.1) Presentation of the 3 relations used

In this paper we study and discuss three radar relations (see references in the following section where parameters' values are discussed):

(i) The R Z h -relation: b h aR Z = (11)
(ii) The dp K Rrelation:

d dp cK R = (12)
(iii) The

dr h Z Z R - - relation: g dr f h Z eZ R = (13)
Here it should simply be mentioned that some authors express dr Z in dB (ex. Illingworth and Blackman 2002) while others keep it in linear scale (Leinonen 2012, Fuigueras i Ventura 2012) as we did here. We focus on these relations because they are the most commonly used, especially the first two. Other ones can be found in the literature such as the

dr dp Z K R - - (Ryzhkov

3.2) "Climatic" analysis

First we consider all the available data to compute average values of the radar relations' parameters. In this study they will be denoted "climatic" values keeping in mind that they were obtained only with two months of data which is limited for such a statement. The plots yielding to the parameters' estimates of the R Z h -and dp K Rrelations for the Pars#2 data are shown in Fig. 7. The quality of the regression is good with coefficient of determination greater than 0.8. Estimates are reported in Table 1 for all radar relations and disdrometers. They are consistent with those commonly reported in the literature.

With regards to the

3.3) Event based analysis

In this section the radar relations are estimated independently for each event with Pars#2 data. Similar results are found for the other devices. They are computed only for events with more than 10 points available to carry out the linear regression and whose cumulative depth is greater than 3 mm. These conditions result in 17 events studied for the R Z h -and

dr h Z Z R - -
relations and 16 for the dp K Rrelation. Figure 8.a displays b vs a for the R Z hrelation. No tendency was found between the event general features (duration, total rainfall amount…) and the parameters of the radar relations. It appears that the ranges of values is rather great and greater than the ones reported by Jaffrain and Berne (2012) in their analysis of 36 events over the Lausanne area (Switzerland). It is of the same order of magnitude as those found by [START_REF] Matrosov | X-Band Polarimetric Radar Measurements of Rainfall[END_REF] who found a ~ [94-624] and b ~ [1.34-2.36] when analysing 15 rainfall events observed with a transportable X-band radar. No trend is visible with this data. The variability is smaller for the exponent b than for the pre-factor a, with coefficients of variation (CV) respectively equal to 14% and 33%. Here there is no obvious trend between c and d, maybe a slight positive correlation (as in Jaffrain and Berne 2012 reported) but nothing of a statistical significance. Similarly to the R Z h -relation, CV is smaller for the exponent d than the pre-factor c (28% vs. 41%). The variability between events seems slightly greater for the parameters of the dp K Rrelation than for the ones of the R Z h -relation. It is likely to be due to a greater sensitivity of the dp K Rrelation than the R Z h -one to the disparities in terms of DSD. For these selected events we find CV=35% for D m and 34% for N t , which is only slightly smaller than the values found for all the events. Figure 8.c, 8.d and 8.e respectively display the e (CV = 40%) , f (CV = 9.3%) and g (CV = 53%) parameters of the

dr h Z Z R - -
relation vs. one another. There simply seems to be a negative correlation trend between e and f. The CV are greater than for the other relations except for f.

3.4) Intra-event analysis

In this section, we go a step further and investigate the variability of the R Z h -and dp K Rradar relations at the highest available temporal resolution with the help of data collected with Pars#2 device. In order to achieve this, the parameters a, b, c and d are computed for each time step considering a moving window of 40 time steps (i.e. 10 minutes before and 10 minutes after). Again they are only computed if 10 points are available to perform the linear regression. This is done for the two selected events of 29-31 May 2016 and 17 June 2016.

With regards to the 29-31 May event, Fig. 9.a displays the computation of a and b ( R Z hrelation) for the whole event. Figure 9.c and d shows the temporal evolution of the "local" (i.e. with the moving window) of respectively a and b; while a scatter plot of b vs. a is in Fig. 9.b. Note that the colour of the points in Fig. 9.b corresponds to the time within the event at which a and b are estimated. It means that two points with similar colours correspond to the same moment of the storm. Figure 10 displays the same information but for the dp K Rradar relation. It should be noted that the number of time steps for which it is possible to compute c and d (only 6% of time steps for the May event) is more reduced than for a and b. This is due to the condition 2 . 0 ≥ dp K °.km -1 which requires significant rain rates. This explains the lack of information found on Fig. 10.c and 10.d. Figures 11 and12 are the equivalent of respectively 8 and 9 but for the 17 June event. For this event the dp K Rcan be fitted on most of the event.

It can be seen that at the event scale, the linear regression illustrating the radar relations are good meaning the underlying radar relations are reliable. Second it appears that the variability observed between events (see section 4.3) is also present within an event. Indeed the temporal evolutions of radar parameters exhibit some quick variability, even with the use of the moving window technique. Figure 13 displays the temporal evolution of a and b for the 29-31 May event for the three co-located disdrometers. It highlights that parameters a and b exhibit a very local variability, which if taken into account could improve rainfall retrieval with radars. Furthermore although with some differences (local variability, sampling variability, device functioning), it appears that similar patterns are found. It suggests that the temporal variability of a and b is associated with a variability of the underlying physical processes and not instrumental or sampling uncertainties. For the 17 June event a negative correlation between the parameters of a given radar relation ) is visible for the first 50 minutes. Such obvious correlations are not so visible for the May event which is longer. Let us note that for the R Z h -relation one can distinguish some correlation on small portions of the event (some "lines" on Fig. 9.b with similar colours); suggesting that such correlation would only be valid for short periods of time, typically few tens of minutes.

An attempt was made to establish a correlation between the mass weighted diameter D m and the radar parameters. It was computed in a similar way with a moving window for each time step. Results are presented in Table 2. They are contrasted and do not enable to obtain strong conclusions. Limited correlation is found for the R Z h -relation. For the dp K Ran anticorrelation is noted with c, while contrasted results according to the event are found for d, suggesting a lower dependence on D m . Similar results were found for Pars#1.This issue should be further investigated to yield more robust conclusions. 2: Correlation between the radar parameters and the mass weighted diameter (D m ) for the two selected events with Pars#2 4) Consequences on rainfall retrievals with radars.

4.1) Methodology

In this section the portion of radar algorithms converting filtered radar quantities to rainfall rates are mimicked with the help of disdrometer measurements. More precisely: (i) Pseudo-radar quantities (Z, Z dr , K dp ) are estimated with the help of the instantaneous DSD and Eqs. 1 and 2. The rain rate is computed from the raindrops collected. (ii) The parameters of the radar relations are estimated at various temporal scales ("climatic", event, and local) with the help of log-regressions, as described in section 3. (iii) Rain rates that would have been obtained with radars are simulated by implementing the radar relations. Four of them are tested; the three studied ones in previous sections and an additional hybrid one:

- b h rad a Z R / 1       = (14)
-

d dp rad cK R = (15)
-

g dr f h rad Z eZ R = (16)
-A rather standard hybrid model which consists in using a R Z hrelation for low rain rates When "local" estimates of the radar relations' parameters are not available the "event" ones are used. When "event" estimates are not available, "climatic" ones are used. All computations are carried out at 30 s time steps. The issue of rain profiling algorithm using all the values along a full radar beam are not explored here [START_REF] Anagnostou | High-Resolution Rainfall Estimation from X-Band Polarimetric Radar Measurements[END_REF][START_REF] Park | Correction of Radar Reflectivity and Differential Reflectivity for Rain Attenuation at X Band. Part I: Theoretical and Empirical Basis[END_REF][START_REF] Testud | The rain profiling algorithm applied to polarimetric weather radar[END_REF]) and only local algorithms are mimicked. These radar relations, which are used in all radar algorithms based on pre-defined scattering relations, are sensitive to the choice the drop shape model and the discrepancies with the unknown true shapes yield some significant uncertainty on the rain rate retrieval. It should be highlighted that the developed approach enables to artificially avoid this very complex problem by simply using the same to compute the reconstructed pseudo-radar quantities and the radar relations. Obviously this cannot be the case with true radar data since the drop shape is then unknown.

With these algorithms, we obtain 12 (4 types with either "climatic", "event" or "local" parameters for the radar relations) pseudo-radar rainfall time series (R rad ) for each event and they are compared with the rainfall time series obtained directly from drop measurements (R). The first step of the comparison consists in performing the same multifractal analysis that was done in section 3 and checking whether similar results are obtained. In a more standard way, we also computed scores commonly used for such tasks (Diss et -The Normalized Bias (NB) whose optimal value is 0:

1 - = R R NB rad -
The correlation coefficient (corr) which varies between -1 and 1 and whose optimal value is 1:

( )( ) ( ) ( ) ∑ ∑ ∑ ∀ ∀ ∀ - - - - = i rad i rad i i i rad i rad i R R R R R R R R corr 2 , 2 ,
-The Nash-Sutcliffe model efficiency coefficient (Nash), which varies between ∞ and 1 and whose optimal value is 1:

( ) ( ) 2 2 , 1 ∑ ∑
∀ ∀ - - - = i i i i i rad R R R R Nash
-The Root mean square error (RMSE), which varies between 0 and ∞ + and whose optimal value is 0:

( ) N R R RMSE i i i rad 2 , ∑ ∀ - =
-The Slope and Offset of the orthogonal linear regression. The optimal values are respectively 1 and 0.

Where R rad and R correspond respectively to the reconstructed pseudo-radar data and rain drop derived rainfall. <> denotes the average. Time steps (index i in the previous formulas) of all the events are used in the sum for each indicator.

4.2) Overall results

Results are presented for the data collected with Pars#2 device. Similar ones are found for the other disdrometers. Figure 14 displays the scatter plot of R rad vs. R for the 12 algorithms tested. Table 3 shows the output of the multifractal analysis performed on the various fields. For instance if we consider the R Z halgorithm with "climatic" parameters, we find α=1.84 and C 1 =0.25 while for the actual series we have α=1.64 and C 1 =0. [START_REF] Krajewski | Radar -Rainfall uncertainties : Where are We after Thirty Years of Effort[END_REF]. It means that the series obtained with the pseudo radar algorithm exhibits too strong multifractality and intermittency with regards to the actual one, resulting in unrealistic stronger extremes. Standard comparison scores are presented in Table 4. These results enable to sort the various algorithms in terms of performance.

Considering a fixed temporal scale for the computation of the radar relations' parameters (i.e. "climatic", "event" or "local") it appears that the hybrid algorithm is the one performing best. Indeed the UM parameters estimates are closer to actual values for the time series computed with this algorithm and the standard scores are better. This confirms the results found by Figueras i Ventura et al. ( 2012) who analysed at the hourly time scale radar data from France comparing them with rain gauge data. Then follow respectively the

dr h Z Z R - - and R Z h - algorithm. The behaviour of the dp K R -
relation is trickier to analyse. Indeed in terms of scores, it seems to perform similarly or best than the other algorithms except for NB and total cumulative depth. This is related to the behaviour for "small" intensity for which there is a strong negative bias. For instance with the climatic average, we have NB equal to almost -0.3 for time steps when the rain rate is smaller than 10 mm.h -1 whereas it is lower than 0.1 with rain rates greater than 10 mm.h -1 (which is visible on Fig. 14). Fitting it using the whole range of available rain rates yields very low regression quality. This is consistent with the more complex behaviour of the scattering coefficient for computing K dp vs. D (Fig. 1.c). Indeed the relation becomes clearer only for rather large drops which are visible during heavy rainfall periods, whereas the aforementioned truncation of drops with D <1 mm cannot be forgotten for low intensity rainfalls For each algorithm it appears that the performance improves when going from using "climatic" to "event" to "local" estimates of the radar relations' parameters. This effect was expected and is now quantified. The R Z h -algorithm is especially sensitive to this effect and its performance strongly improves when a and b parameters are tuned more locally. With "local" estimates, its performances become comparable with the

dr h Z Z R - -
ones which is not the case with "climatic" estimates. [START_REF] Matrosov | X-Band Polarimetric Radar Measurements of Rainfall[END_REF] noticed this effect for the R Z halgorithm at the event scale. This improvement could somehow be expected from Fig. 1.a-b which displays the backscattering coefficient as a function of the equivolumic drop diameter. Indeed, the relation being straightforward, once parameters of the R Z h -are fitted with a local DSD which limits its influence on the variability (i.e. the main source of uncertainty), it is not surprising that the R Z hperforms well. This improvement is less pronounced for the hybrid algorithm mainly because its performances are already good with "climatic" estimates of radar parameters, and the hybrid algorithm somehow already tunes the parameters according to rain rates by using two different relations. Although it is practically complicated these results make a case for the use (when possible) of tuned radar relations' parameters, if not in real time at least during re-analysis. The rather good performances of the simple R Z halgorithm show at first glance that when no dual-polarization radar can be installed, the combination of a disdrometer used to estimate in real time the parameters of the R Z hrelation and a single polarization radar can yield good rainfall estimates. However, one should keep in mind that tuned parameters rely on DSD measurement at a single point not taking into account its variability in space.

In this paper, multifractal techniques are used to assess the quality of pseudo-radar algorithms The closeness of UM parameters assessed on radar time series to the ones obtained on the actual rainfall time series may help to assess the performance of the algorithms. The use of such techniques is very interesting because it intrinsically validates the outcome across all studied scales and is not limited to a single resolution as it is the case for the standard scores also computed in this paper. Such practical applications of scaling techniques were already developed to assess the quality of high frequency long rainfall time series [START_REF] Hoang | Assessing the high frequency quality of long rainfall series[END_REF]). The analyses carried out here are another example. 4: Standard scores of the rain rates computed via reconstructed pseudo-radar data and algorithms (R rad ) vs. the one obtained from direct drop measurements (R) (total cumulative depth = 204 mm) for the data collected with the help of Pars#2 device.

Radar algorithm Parameter estimation

Radar conversion

β R 2 α C 1 H b h rad a Z R / 1       = 1.

4.3) Focus on the two selected events

In this section, we illustrate more precisely the general results discussed quantitatively in the previous section on the performance of the various radar algorithms. This is done by looking into more details at the time series for the two selected events and Pars#2 data. Figure 15.a displays the temporal evolution of the rain rates obtained via reconstructed pseudo-radar data and the three R Z halgorithms ("climatic", "event", and "local") along with the direct estimation of R. The cumulative depth is shown in Fig. 15.b, where it is visible that the "climatic" estimates yield the worst cumulative depth (an underestimation in this case). A closer look at the temporal evolution of rain rates confirms the improvement brought by the "local" (in time) estimates of the radar relations' parameters. For example during the pseudo-peak slightly before 4.7 h (Fig. 15.c) the rain rate is below 5 mm.h -1 and the use of either the "climatic" or "event" relations yields strong overestimation, while the local parameters enable to better represent actual rain rates. Similar conclusion is found around the peak at 11.9 h, with an underestimation by the "climatic" and "event" based algorithms. The hybrid models yield similar results, but results are not shown here because of the actually limited number of time steps for which a local dp K Rcan be computed (6%). Figure 16 (left column) displays the temporal evolution of the rain rates computed via reconstructed pseudo-radar data and the 4 algorithms with "climatic", "event", and "local" parameters along with the direct estimation of R. Figure 16 (right column) exhibits the corresponding cumulative depth. Again the "local" estimates of radar relations' parameter provide the best estimates. In this case the "climatic" parameters yield an overestimation expect for

dr h Z Z R - - algorithm.
It is especially visible for the R Z halgorithm. The improvement brought by the "event" based and "local" estimates are visible on the successive peaks (Fig. 16.a).

5) Conclusion

In this paper, the rainfall data measured by disdrometers of Ecole des Ponts ParisTech during the two extremely wet months of May and June 2016 are analysed. At first, this rainfall data exhibit very good quality of scaling behaviour over scales ranging from 30 s to few hours.

Classical relations between rain rate and the horizontal reflectivity, the differential phase shift and the differential reflectivity (

R Z h -, dp K R - , dr h Z Z R - -
) are investigated with the help of pseudo-radar quantities reconstructed from the disdrometer data. This reconstruction requires that both the DSD and the drop centre distribution are homogeneous. In this simplistic framework, the characteristic parameters of these power-law relations are computed at various scales: "climatic", "event" and "local" with moving time window. It appears that these parameters exhibit a strong variability from one event to the other and even within an event.

Finally these radar relations are used to reconstruct rain rate time series that would be observed with the help of weather radars. In addition to the previously mentioned relations, the standard hybrid model relying on a R Z h -relation for low intensities and a dp K Rrelation for large intensities is tested. An innovative methodology based on assessing the quality of an algorithm on its ability to reproduce the scaling behaviour of actual time series was developed. This ensures that results are valid across scales. The main conclusion is that the hybrid algorithm is the one performing best for the pseudo-radar quantities, whatever the temporal scale of the computation of the radar relations' parameters. However, the disdrometer algorithm estimating the drop shape unfortunately truncates the computation of the pseudo K dp at D min =1 mm and therefore the contribution of the prolate drops. It also appears that the performances of all algorithms significantly improve by using parameters computed over shorter period of time, but this unfortunately corresponds to less predictive value. Not too surprisingly, the famous R Z hrelation tuned with "local" parameter estimates exhibits very good performance, but this seems to be a somewhat artificial result from the power-law-like behaviour of the radar cross-section σ(D). Overall, the analyses presented in this paper call for extended studies that would take into account (including in real time) the spatial variability of the DSD and the possible clustering of drops, as well as the contribution of the prolate drops, that could not be addressed in this paper.

  7 for D > 5 mm, and linearly related to D with a slope equal to 0.075 in for intermediate D. (ii) Drop orientation is the same as in Leinonen et al. (2012); i.e. drops are partially aligned and a normal distribution (mean and standard deviation respectively equal to 0° and 7°, in agreement with the findings of Bringi et al. 2008) characterizes the angle of the symmetry axis. Figure 1 displays typical curves of the scattering coefficients with respect to the particular diameter. As expected B σ is smaller for vertical polarization than horizontal ones and the differences get greater with increasing D as drops get more oblate (Fig. 1.a). B σ vs. D is plotted in log10-log10 plot in Fig. 1.b. It enables to point out a slight interest of using T-Matrix computations code rather than rely on a pure Rayleigh approximation in which 6

2 )

 2 Data for May and June 2016

Finally

  R Z h -relationship we find a = 195 and b = 1.78 for Pars#2. The value of a is significantly greater for Pars#1 with similar value of b. For PWS, values found for a and b are respectively greater and smaller than with Pars#2 data. There is a wide variety of values reported in the literature for these parameters. Typical ones are a=200 and b=1.6 (the widely used parameter set since Marshal and Palmer, 1948) or a=300 and b =1.4 (which are the Nexrad standard values in USA, Fulton et al. 1998). Matrosov et al. (2002) found values of a = 250 and b = 1.68 in the X-band using experimental raindrop size distributions. Jaffrain and Berne (2012b) reported values ranging from 220 to 320 for a and from 1.4 to 1.6 for b after analysing 36 rainfall events with disdrometer data. Hence it appears that the values found here for a and b are located respectively in the lower range and upper range of reported values in the literature. This apparent slight discrepancy is actually not surprising. Indeed Verrier et al.(2013) showed that a and b are dependent on the scale at which they are computed. Given that they are computed here at 30s which is much smaller than the minute or few minutes commonly used, the differences are expected. Values are in agreement with those found by Verrier et al. at this temporal scale. It should also be reminded that one should be careful when comparing values because they are also dependent on the regression method (mainly the choice to express quantity in log or not, and the regression type) used to extract them(Campos and Zawadski, 2000; Jaffrain and Berne 2012; Verrier et al. 2013). For the dp K Rrelation, we find c = 12.9 and d = 0.75 for Pars#2. Very close estimates are found with Pars#1 data and close ones with PWS data. Typical values for the X-band are c = 19.9 and d = 0.85 (Beard and Chuang, 1987) or c = 23.2 and d = 0.79 (Brandes at al. 2002). Jaffrain and Berne (2012b) found also in the X-band c ~ [11-14] and d ~ [0.69-0.75] according to the studied event. Leinonen (2012) found in the C-band c=21.0 d =0.72, meaning that the parameter c would be similar to the one found here, had the study be done in X-band. Anagnostou et al. (2010) reported values of c = 19.26 and d = 0.85 in the X-band. Matrosov et al. (2002) found values of c = 12.3 and d = 0.81 in the X-band. Hence the values obtained in this analysis are in agreement with the ones found in the literature. 0.023, f = 0.77 and g = -4.05 for Pars#2. Comparable values are retrieved with the other devices. This relation is less used by the community than the other two. Leinonen et al. (2012) reported values of e = 0.012, f = 0.822 and g = -2.28. Ryzhkov et al. (2005) found e ~ [0.0067-0.0159], f ~ [0.73-1] and g ~ [-4.73, -1.03] according to the drop shape model, and the type of data used (experimental or modelled DSD).

Figure 8 .

 8 b displays d vs. c for the dp K Rrelation. Again the range of obtained values is wider than in Jaffrain and Berne (2012b).

  (a and b Fig 11.b or c and d Fig 12.b

Figures:Figure 1 :Figure 2 :

 12 Figures:

Figure 3 : 2 Figure 4 :

 324 Figure 3: (a) Scatter plot of the 30 s time step rain rate for Pars#1 vs. Pars#2. (b) Same as (a) but for PWS vs. Pars#2

Figure 5 :

 5 Figure 5: For the 17 June 2016 event with device Pars#2. (a) Temporal evolution of the rain rate in mm.h -1 with 30 s time steps (solid black, the observation time step) and 5 min time steps (dashed red). (b) Cumulative rainfall depth in mm. (c) Temporal evolution of the DSD N(D).

Figure 6 :

 6 Figure 6: For the data collected with Pars#2 (a) Spectral analysis, i.e. Eq. 6 in log-log (b) TM analysis, i.e. Eq. 7 in log-log. Samples of size 64, corresponding to 32 min, are used.

Figure 7 :

 7 Figure 7: For the Pars#2 device, estimation of the exponents of the power-law relations R Z h -(a) and dp K R -(b) by regression in a log-log plot, considering all the available events.

Figure 8 :

 8 Figure 8: For the Pars#2 device, plot of the radar relations' exponents vs. one another for the R Z h -(a), the dp K R -(b) and the

Figure 9 :

 9 Figure 9: R Z hrelation for the 29-31 May 2016 event with Pars#2. (a) Computation of the parameters a and b at the event scale (eq. 10 in a log-log plot). (b) Scatter plot of b vs. a estimated for each time step (20 minutes moving window). The colour of the points corresponds to the time within the event at which the parameters are estimated (c) Temporal evolution of a. (d) Temporal evolution of b.

Figure 10

 10 Figure 10:

Figure 11 :

 11 Figure 11: R Z hrelation for the 17 June 2016 event with Pars#2. (a) Computation of the parameters a and b at the event scale (eq. 10 in a log-log plot). (b) Scatter plot of b vs. a estimated for each time step (20 minutes moving window). The colour of the points corresponds to the time within the event at which the parameters are estimated (c) Temporal evolution of a. (d) Temporal evolution of b.

Figure 12

 12 Figure 12:

Figure 13 :

 13 Figure 13: R Z hrelation for the 29-31 May 2016 event with the three disdrometers. (a) Temporal evolution of a during the full event. (b) Temporal evolution of b during the full event. (c) and (d) Zoom of (a) and (b) during a shorter period.

Figure 14 :

 14 Figure 14: Scatter plot of the rain rates computed via reconstructed pseudo-radar data and algorithms (R rad ) vs. the one obtained from direct drop measurements (R) for the data collected with the help of Pars#2 device.

Figure 15 :

 15 Figure 15: For the 29-31 May 2016 event and Pars#2 data. (a) Temporal evolution of the rain rates obtained via reconstructed pseudo-radar data and the three R Z halgorithms ("climatic", "event", and "local") along with the direct estimation of R. (b) Corresponding cumulative depth. (c) Zoom of (a) for a shorter period of time. (d) Zoom of (a) for another shorter period of time.

Figure 16 :

 16 Figure 16: For the 17 June 2016 event and Pars#2 data. (a) Temporal evolution of the rain rates obtained via reconstructed pseudo-radar data and the three hybrid algorithms ("climatic", "event", and "local") along with the direct estimation of R. (b) Corresponding cumulative depth.

  et al. 1995), but they are not explored here.

	The parameters a, b, c, d, e, f, and g are estimated using the measured DSD at the observation
	time step of 30 s. For the	R Z h -relation, an orthogonal regression is carried out for time
	steps with	R	≥	0	.	2		mm.h -1 . An orthogonal regression minimizes the orthogonal distance from
	the data points to the fitted line, contrary to the ordinary linear regression which minimizes
	the vertical distance. For the	R -	K	dp	relation, an orthogonal regression is carried out for time
	steps with			≥	0	.	2				R	-	Z	h	-	Z	dr	relation, a linear regression minimizing the
	quadratic error considering only the time steps with	R	≥	0	.	2	mm.h

dp K

°.km -1 . For the -1 is used. Computations are carried out only when more than 10 points (i.e. time steps) are available.

Table 1 :

 1 "Climatic" parameters of the radar relations obtained with the disdrometer data used in this study.

	Radar relation	Parameter	Pars#1	Pars#2	PWS
	Z = h	aR	b			a b	243 1.76	195 1.78	240 1.52
	R =	cK	dp	d			c d	12.6 0.80	12.9 0.75	14.7 0.82
	R =	eZ	f h Z	dr	g	e f	0.018 0.80	0.023 0.77	0.012 0.87
										g	-3.95	-4.05	-4.52

  al., 2009; Emmanuel et al., 2012; Figueras I Ventura et al., 2012; Krajewski et al., 2010; Moreau et al. 2009):

Table 3 :

 3 Results of the multifractal analysis performed on the rain rates computed via reconstructed pseudo-radar data and algorithms (R rad ) vs. the one obtained from direct drop measurements (R) for the data collected with the help of Pars#2 device.

	Climatic												69	0.98	1.84	0.25	0.59
		rad R =	cK	dp	d			1.93	0.99	1.44	0.20	0.63
		rad R =	eZ	h	f	Z	dr	g	1.41	0.99	2.05	0.23	0.44
		Hybrid					1.81	0.99	1.76	0.21	0.60
	Event	rad R	=	  	h a Z	1   	/	b	1.58	0.99	1.79	0.24	0.51
		rad R =	cK	dp	d			1.97	0.99	1.42	0.19	0.63
		rad R =	eZ	h	f	Z	dr	g	1.55	0.99	1.91	0.22	0.49
		Hybrid					1.97	0.99	1.60	0.19	0.65
	Local	rad R	=	  	h a Z	1   	/	b	1.74	0.99	1.69	0.22	0.56
		rad R =	cK	dp	d			1.94	0.99	1.47	0.19	0.63
		rad R =	eZ	h	f	Z	dr	g	1.82	0.99	1.72	0.20	0.59
		Hybrid					1.96	0.99	1.61	0.19	0.65
	For the actual time series			1.98	0.99	1.64	0.19	0.65

Acknowledgements

Authors greatly acknowledge partial financial support form the Chair "Hydrology for Resilient Cities" (endowed by Veolia) of Ecole des Ponts ParisTech,


C Lebas

S Sahuguede

A Julien-Vergonjanne

P Combeau

L Aveneau

Infrared and Visible links for medical Body Sensor Networks

Keywords: LiFi for medical applications, infrared transmissions, channel modelling

Our previous studies focused on channel simulation and performance evaluation of optical wireless links for medical body sensor networks. This allowed us to increase our expertise in this field and to propose here a full optical wireless bidirectional system named as LiFi communication system for medical monitoring applications. The full duplex bidirectional communication is based on an infrared uplink and visible downlink. The studied scenario considers a patient wearing both an infrared transmitter connected to a medical sensor and a visible light receiver. The infrared receivers and visible light emitters are included in light fixtures at the ceiling. The study relies on channel modelling considering patient body features, locations and movements. Using ray-tracing technique associated to Monte-Carlo method, we determine the overall performance in terms of outage probability. We apply our methodology considering various medical sensors and show that it is possible to satisfy the overall quality of service in terms of bit error rate and data rate.

I. INTRODUCTION

Optical Wireless Communication (OWC) in the nanometerrange from ultra-violet to InfraRed (IR) has the potential to become a complement or alternative to RF technology [START_REF] Uysal | Optical Wireless Communications -An Emerging Technology[END_REF]. Recently, OWC has gained popularity in the visible range domain named Visible Light Communication (VLC), due to the dual use of Light Emitting Diode (LED) devices for both illumination and wireless communications. This technology, named as LiFi communication offers many benefits such as unlicensed spectrum, high bandwidth, simple and cheap frontends and no electromagnetic interference. This latter benefit is important when considering portable medical monitoring devices incorporating wireless transmission technologies. Indeed, the use of RF technologies around the human body can be problematic in terms of impact on health, security and confidentiality. Therefore, OWC technology can be a solution to provide low cost, efficient and safe communications for medical devices, between medical equipment or between patients and monitoring stations [START_REF] Murai | A novel visible light communication system for enhanced control of autonomous delivery robots in a hospital[END_REF][START_REF] Cheong | Hazardless Biomedical Sensing Data Transmission Using VLC[END_REF][START_REF] Song | Indoor hospital communication systems: An integrated solution based on power line and visible light communication[END_REF][START_REF] Yy | Mobile health-monitoring system through visible light communication[END_REF][START_REF] Chevalier | Optical Wireless Links as an Alternative to Radio-Frequency for Medical Body Area Networks[END_REF][START_REF] Cahyadi | Patient monitoring using Visible Light uplink data transmission[END_REF][START_REF] Al-Qahtani | A non-invasive remote health monitoring system using visible light communication[END_REF][START_REF] Dhatchayeny | EEG biomedical signal transmission using visible light communication[END_REF][START_REF] Adiono | Visible light communication system for wearable patient monitoring device[END_REF].

Involved in OWC field for medium and low-rate medical applications, we have developed skills and dedicated tools to address the modelling of indoor optical wireless channels applied to Body Sensor Networks (BSN) [START_REF] Chevalier | Optical Wireless Links as an Alternative to Radio-Frequency for Medical Body Area Networks[END_REF]. By using robust mathematical models (Monte-Carlo, Metropolis, etc.), our team has developed RaPSor (Ray Propagation Simulator) software

which is an open and extensible tool, based on the Netbeans platform for the modelling of Infrared and Visible links [START_REF] Behlouli | MCMC Methods for Realistic Indoor Wireless Optical Channels Simulation[END_REF]. We have also designed portable OWC devices to carry out experimental test-beds for low data rates infrared uplinks exploiting indoor diffuse reflections [START_REF] Bas | Theoretical and Experimental Approach for the Design of an Optical Wireless Physical Activity Monitoring System[END_REF].

In this study, we propose a full optical BSN based on bidirectional LiFi communication system using IR technology and VLC for indoor health monitoring. A diffuse IR transmitter ensures the uplink between a wearable medical sensor and receivers included in VLC-ready light fixtures. Visible downlinks from the luminaires allow sending medical data to a VLC receiver also carried by the monitored patient. In the study, we focus on parameters affecting IR and VLC channel behavior related to the presence of patient body wearing the IR/VLC transceivers and to the locations and movements inducing orientation variations of the IR emitter and VLC receiver. System consumption and thus emitted power is also of main concern in order to meet power standards and extend usage time.

We investigate these issues by simulation using RaPSor and discuss the overall IR / VLC performance considering various medical monitoring constraints and lighting conditions. In particular, we determine whether it is possible to establish VLC transmissions with very low power levels so that the light is quite not perceptible by the human eye.

The paper is organized as follows. The system description we investigate is presented in section 2. Section 3 focuses on the theoretical analysis of IR uplink and VLC downlink channels considering patient body presence and movements. Section 4 reports the study of overall performance before conclusion in section 5.

II. SYSTEM DESCRIPTION

To evaluate the overall performance of IR/VLC for health monitoring, we consider the scenario illustrated in Fig. 1.

The studied area has dimensions of (5×5×3) m 3 (length, width, height) which is a typical one. Four LED based luminaires are located at the ceiling (coordinates are indicated in TABLE I) oriented towards the floor. The LED radiation patterns are modeled by a Lambertian one with an order m=1 that corresponds to a half-power angle of 60°. In the studied scenario, the visible downlinks are established between LEDs and a VLC receiver integrated into a tablet for example worn by a patient moving in the environment. To represent this scenario we will model the patient body considering random positions in the environment. In addition, the receiver is located at a distance of 20cm from the body and 1.2m from the floor. As the patient is moving, it is possible that the orientation of the receiver varies. Receiver orientation can be defined through polar and azimuthal angles. Consequently, we will consider random values of these two angles taken in the ranges defined in TABLE I, which also lists the values taken for responsivity 𝑅, Field-Of-View (FOV) and physical area of the photodetector which are typical values for VLC receivers. For simplicity and without loss of generality we will not model furniture.

TABLE I. VLC system parameters

Concerning IR uplink transmission, we suppose an emitter integrated into a medical device worn by patient at the shoulder (1.5m from the floor), and four IR receivers located at the same coordinates as the luminaires (see TABLE II which also reports the typical parameters taken for IR receivers such as responsivity 𝑅, FOV, physical area). In addition, we use one of the simplest combining techniques which is Switching Combining (SwC) selecting only once the information if one or several receivers correctly receive it.

Besides, because patient movements, the wearable emitter orientation can vary, which is taken into account through polar and azimuthal angles (see TABLE II). Note that as the sensor is worn at the shoulder, the variation range of azimuthal angle is restricted compare to VLC receiver case. Furthermore, it is supposed that the IR emitter pattern is Lambertian with an optimal m order equal to 2 that is an half-power angle of 45° [START_REF] Behlouli | Efficient simulation of optical wireless channel, Application to WBANs with MISO link[END_REF].

Among the specifications of medical remote monitoring applications, it can be noted that the bit rates are generally less than 10 Mbps but with a rather high quality of service [START_REF] Chevalier | Optical Wireless Links as an Alternative to Radio-Frequency for Medical Body Area Networks[END_REF]. For this low data rates and considering typical impulse responses length and mean delay spreading, we can neglect inter-symbol interference in both IR and VLC cases. The channels are characterized by the optical DC gain, defined as :

𝐻(0) = ∫ ℎ(𝑡)𝑑𝑡 = 𝐻 0 +∞ -∞ (1)
where ℎ(𝑡) represents the channel impulse response.

Patient mobility induces different link configurations related to patient positions and orientations of VLC receiver and IR emitter. Consequently, the optical DC gain is a random variable. In the following, we study the channel statistical behavior taking also into account the human body model. scenario by using stochastic method of Monte Carlo, associated to ray tracing algorithm [START_REF] Behlouli | MCMC Methods for Realistic Indoor Wireless Optical Channels Simulation[END_REF][START_REF] Behlouli | Efficient simulation of optical wireless channel, Application to WBANs with MISO link[END_REF]. For all the simulations, in order to manage tradeoffs between calculation times and accuracy, we consider three reflections per optical beam, which is a classical approach for non-directed Line-Of-Sight (LOS), or diffuse transmissions. Furthermore, the wall material is plaster surface with reflectivity fixed to 0.73 and the floor material is plastic one with a reflectivity set to 0.18 [START_REF] Lee | Indoor Channel Characteristics for Visible Light Communications[END_REF].

TABLE II. IR system parameters

Due to patient mobility inside the room, we will consider a set of impulse responses h(t) of IR and VLC links obtained with ray-tracing simulations. A number of 1000 random positions has been taken since it has been verified that it permits to have convergent results. The angle variation step for the orientation changes of the IR transmitter and the VLC receiver is set to /6. . The channel behavior will be then characterized by means of H0 statistical distribution in terms of probability density function (PDF) and cumulative density function (CDF).

Our objective is to analyze the link reliability from CDF. It will be higher the more the 𝐶𝐷𝐹(𝐻 0 ) will be small for a given 𝐻 0 since the goal is to reach high DC gain 𝐻 0 values with a high probability.

B. IR channel

To study the impact of the human body on the IR channel, we have to model its geometry and its reflectance properties. It is important to take into account the width of the body because of possible blocking effects since receiver is supposed to be at the shoulder; however, level of details is not important for evaluating channel gain when considering low data rates [START_REF] Bas | Human body impact on mobile visible light communication link[END_REF][START_REF] Behlouli | Impact of Physical and Geometrical Parameters on Visible Light Communication Links[END_REF]. We thus consider the body shape as illustrated in Fig. 2 but with a given width of 20 cm. For the reflectivity, since there is a large range of possible values depending on the skin and the clothes, we have taken two extreme values i.e. ρb1=0.1 and ρb2 =0.9.

These two cases are compared in the following to the reference case where neither the body nor the variations of emitter orientation are taken into account, as summarized in Fig. 2. To determine the DC gain we define 𝐻 𝑖 as being the DC gain between given position and orientation of the IR emitter and the receiver 𝑅 𝑖 with i varying from 1 to 4. In order to take into account the SwC method, we define 𝐻 0-𝐼𝑅 which corresponds to the best DC gain considering the four active receivers:

𝐻 0-𝐼𝑅 = Max[𝐻 1 , 𝐻 2 , 𝐻 3 , 𝐻 4 ] (2)
As the emitter position is random, 𝐻 0-𝐼𝑅 follows a statistical distribution 𝑝(𝐻 0-𝐼𝑅 ). The CDF of the gain 𝐻 0-𝐼𝑅 can be thus expressed as:

𝐶𝐷𝐹(𝐻 0-𝐼𝑅 ) = ∫ 𝑝(𝐻 0-𝐼𝑅 )𝑑𝐻 0-𝐼𝑅 𝐻 0 -∞ (3)
We have reported in Fig. 3 the evolution of 𝑃𝐷𝐹(𝐻 0-𝐼𝑅 ) and 𝐶𝐷𝐹(𝐻 0-𝐼𝑅 ) for all the studied cases defined in Fig. 2. It can be seen that the presence of the body and variations in the orientation of the transmitter affect the statistical distribution and the values of DC gain. This impact is all the more negative as the reflectivity of the body is weak.

C. VLC channel

For the VLC channel, we consider the same body shape as presented for IR but we do not model the width of the body. Actually, the VLC receiver is here supposed to be in front of the body and not over the body as previously for IR emitter.

We compare the same cases as for the IR channel described in Fig. 2 but here orientation variations are related to VLC receiver instead of IR emitter. Cases description:

We can see in Fig. 4 compared to Fig. 3 that the impact of the body and the orientation changes is much more significant for the VLC downlink than for the IR uplink. Another difference is that the effect is important even if the reflectivity of the body is high.

To compare VLC and IR link behaviors, we have plotted in Fig. 5 the two PDF corresponding to case (b)with orientation changes and body reflectivity ρb1=0.1 -for both channels gains 𝐻 0-𝑉𝐿𝐶 and 𝐻 0-𝐼𝑅 . We observe that VLC downlink globally overcomes IR uplink since channel gain values are higher most of the time. However, performance also depends on the emitted power of both links and the receiver characteristics. In the next part, we will investigate the global performance of the dual link IR/VLC to discuss the possibility of establishing downlinks with a VLC emission power that is sufficiently low to consider that the lighting is not very noticeable. This is important because the application must be able to operate whatever the lighting level in the room. For the IR uplink system, the challenge is also to reduce as much as possible the emitted power to decrease the consumption of the device carried by the patient.

IV. SYSTEM PERFORMANCE ANALYSIS

A. Outage probability definitions

As DC gain randomly varies because of mobility, the Signal to Noise Ratio (SNR) also follows a statistical distribution. Thus, we focus on outage probability which is a common metric to evaluate the link performance. It is defined by [START_REF] Goldsmith | Wireless Communications[END_REF]:

𝑃 𝑜𝑢𝑡 = p (𝑆𝑁𝑅 < 𝑆𝑁𝑅 0 ) (4)
Where SNR0 represents the threshold SNR needed to support a given bit data rate 𝑅 𝑏 . Thus, Pout can be seen as a function of SNR0 and is impacted by average emitted power 𝑃 𝑡 , data rate and ambient noise characterized by its power spectral density 𝑁 0 .

Considering On-Off Keying (OOK) modulation, SNR can be expressed as:

𝑆𝑁𝑅 = 𝑃 𝑡 2 𝐻 0 2 𝑅 2 4𝑁 0 𝑅 𝑏 (5)
The photodiode at the reception converts the optical signal to an electrical signal adding thermal and shot noise. Generally, shot noise is the most preponderant due to the ambient photocurrent 𝐼 𝐵 whose value depends on wavelength. A conventional value of 𝑁 0 is given by:

𝑁 0 ≈ 2𝑞𝐼 𝐵 (6)
Were q is the electron quantum charge. We will consider in this study typical values of 𝐼 𝐵 around 5100µA for VLC and 200µA for near-infrared wavelengths [START_REF] Moreira | Characterisation and modelling of artificial light interference in optical wireless communication systems[END_REF].

From the distribution of 𝐻 0 for each channel established previously, we will consider the distributions 𝑝(𝑆𝑁𝑅 𝐼𝑅 ) 𝑝(𝑆𝑁𝑅 𝑉𝐿𝐶 ) obtained using ( 5) and so we can define two outage probabilities 𝑃 𝑜𝑢𝑡-𝐼𝑅 and 𝑃 𝑜𝑢𝑡-𝑉𝐿𝐶 corresponding to the link performance of IR and VLC channels:

𝑃 𝑜𝑢𝑡-𝐼𝑅 = p (𝑆𝑁𝑅 𝐼𝑅 < 𝑆𝑁𝑅 0-𝐼𝑅 ) (7)
𝑃 𝑜𝑢𝑡-𝑉𝐿𝐶 = p (𝑆𝑁𝑅 𝑉𝐿𝐶 < 𝑆𝑁𝑅 0-𝑉𝐿𝐶 )

Where 𝑆𝑁𝑅 0-𝐼𝑅 , respectively 𝑆𝑁𝑅 0-𝑉𝐿𝐶 , is the SNR needed on the IR link, respectively VLC link, to support a given bit rate 𝑅 𝑏-𝐼𝑅 in the uplink, respectively 𝑅 𝑏-𝑉𝐿𝐶 in the downlink.

In order to analyze the performance of the whole bidirectional communication system, we define global outage probability traducing that there is a failure either on one link or on the other. The global outage probability Pout is thus the probability that the SNR drops below the threshold value on any of the link. We can express it as:

𝑃 𝑜𝑢𝑡 = 𝑝 ( 𝑆𝑁𝑅 𝐼𝑅 < 𝑆𝑁𝑅 0-𝐼𝑅 ∪ 𝑆𝑁𝑅 𝑉𝐿𝐶 < 𝑆𝑁𝑅 0-𝑉𝐿𝐶 ) = 𝑝 ( 𝑆𝑁𝑅 𝐼𝑅 < 𝑆𝑁𝑅 0-𝐼𝑅 ) + 𝑝( 𝑆𝑁𝑅 𝑉𝐿𝐶 < 𝑆𝑁𝑅 0-𝑉𝐿𝐶 ) -𝑝 ( 𝑆𝑁𝑅 𝐼𝑅 < 𝑆𝑁𝑅 0-𝐼𝑅 ). 𝑝(𝑆𝑁𝑅 𝑉𝐿𝐶 < 𝑆𝑁𝑅 0-𝑉𝐿𝐶 ) = 𝑃 𝑜𝑢𝑡-𝐼𝑅 + 𝑃 𝑜𝑢𝑡-𝑉𝐿𝐶 -𝑃 𝑜𝑢𝑡-𝐼𝑅 . 𝑃 𝑜𝑢𝑡-𝑉𝐿𝐶 (9)
Here, Pout is a function of a couple (SNR0-IR, SNR0-VLC). For each link, they depend on average emitted power (𝑃 𝑡-𝐼𝑅 , and 𝑃 𝑡-𝑉𝐿𝐶 ), data rate (𝑅 𝑏-𝐼𝑅 and 𝑅 𝑏-𝑉𝐿𝐶 ) and ambient noise power.

Here, Pout is a function of a couple (SNR0-IR, SNR0-VLC). For each link, they depend on average emitted power (𝑃 𝑡-𝐼𝑅 , and𝑃 𝑡-𝑉𝐿𝐶 ), data rate (𝑅 𝑏-𝐼𝑅 and 𝑅 𝑏-𝑉𝐿𝐶 ) and ambient noise power.

In the following analysis, we will consider that for the studied medical BSN, we have the same required data rate on both links: 𝑅 𝑏-𝐼𝑅 = 𝑅 𝑏-𝑉𝐿𝐶 = 𝑅 𝑏 .

In addition, we assume that the targeted BER are the same for the two channels. So, considering the same OOK modulation for both links this means that 𝑆𝑁𝑅 0-𝐼𝑅 = 𝑆𝑁𝑅 0-𝑉𝐿𝐶 .

B. Global performance analysis

Our objective is to study the possibility of establishing bidirectional full optical IR/VLC link with VLC downlink emission power sufficiently low to consider that the light is not turned on or is not very noticeable. This is important because the application must be able to operate whatever the lighting level in the room including night conditions. For the IR uplink system, the challenge is also to reduce as much as possible the emitted power to decrease the consumption of the worn device in order to have the greatest possible lifetime.

As an example, we focus on the minimal powers for IR and VLC links (𝑃 𝑡-𝐼𝑅 and 𝑃 𝑡-𝑉𝐿𝐶 ) needed to achieve a given global outage probability value using [START_REF] Dhatchayeny | EEG biomedical signal transmission using visible light communication[END_REF], for 𝑅 𝑏 = 1 Mbps and 𝑆𝑁𝑅 0-𝐼𝑅 = 𝑆𝑁𝑅 0-𝑉𝐿𝐶 = 16.1 dB (𝐵𝐸𝑅 ≈ 10 -10 for OOK). This corresponds to classical BSN requirements for medical applications that should support bit rates in the range of 10 kbps to 1Mbps with very low BER [START_REF] Julien-Vergonjanne | Optical Wireless Body Area Networks for Healthcare Applications[END_REF]. To illustrate more precisely VLC power corresponding to real use cases, instead of using 𝑃 𝑡-𝑉𝐿𝐶 we have illustrated results in terms of corresponding average lighting level Emean in lux, within the room at the height of 1.2m [START_REF] Hussain | Visible Light Communication System Design and Link Budget Analysis[END_REF].

All the couples of the minimal values (Emean ; 𝑃 𝑡-𝐼𝑅 ) are plotted in Fig. 6 for respectively 𝑃 𝑜𝑢𝑡 = 0.1%, 1% and 10% when for both links the body reflectivity is set to ρb2=0.9. Fig. 6. Evolution of the minimal needed Pt-IR as a function of the associated minimal illumination Emoy in order to respect BER=10 -10 -ρb2=0.9

Obviously, it can be seen that to guarantee a certain overall quality of service for the IR / VLC transmission, if the IR power is decreased, the lighting level of the room must be increased.

Concerning the value of 𝑃 𝑡-𝐼𝑅 that can be chosen, , we can remark that it is always below the eye safety constraints of around 300mW for IR source with half-power angle of 45°. Moreover, the minimum possible value of 𝑃 𝑡-𝐼𝑅 is as much lower as the target 𝑃 𝑜𝑢𝑡 is high. We can also remark that for this lowest IR power value, the corresponding VLC power must highly be increased, illustrating the tradeoff between both power constraints. On the other side, the minimal value of Emean as the range of lighting values also depend on the 𝑃 𝑜𝑢𝑡 .

Same results have been obtained for the other value of reflectivity coefficient ρb2=0.1 and the conclusions are quite the same for 𝑃 𝑡-𝐼𝑅 . However, Emean required values are significantly increased to satisfy the targeted 𝑃 𝑜𝑢𝑡 . The lowest limit values of 𝑃 𝑡-𝐼𝑅 , respectively Emean, corresponding to a maximal average lighting, respectively a maximal IR power, are reported in TABLE III for different medical application types and for the two extreme values of ρb.

These results highlight that, for data rates lower than 1Mbps, specifications with low IR emitted power can be satisfied, which is of great importance for consumption. In addition, we show that low VLC powers for downlink are sufficient so that the illumination is almost imperceptible by human eye. For some cases (in bold) it may not be possible to respect IR and/or VLC constraints respectively for eye safety and night watch use [21]. These limitations will mostly depend on body's reflectivity value and the targeted Pout.

V. CONCLUSION

The results presented in this article illustrate the impact of two very important parameters for the performance of a communicating mobile full-optical wireless system with transmitters and receivers worn by a person. It is the presence of the body and its optical reflectivity properties as well as variations of direction of the transceivers induced by the movements. This study was conducted for a bidirectional system involving an IR uplink and a VLC downlink in indoor environment. Furthermore, we consider health monitoring applications.

0.02 0.07 0.12 0.17 0.22 3 3.5 First, we have analyzed the statistical behavior of both IR and VLC channel thanks to ray-tracing simulations, allowing modeling LOS and Non-LOS links. The results showed that it is important for the IR or the visible channels, to model the presence of the body and the changes of orientations related to the movements of the body. On the other hand, the impact of the reflectivity of the body is much more significant on the VLC downlink than on the IR uplink for which the transmitter is on the shoulder.

Then, the channel statistics were used to determine the overall performance in terms of outage probability jointly considering IR uplink and visible downlink transmissions.

In particular, we focused on a scenario where the data rates and the required quality of service for both links were identical. The analysis of this overall performance was done by examining the power-related issues for infrared and visible and seeking to minimize both emitted optical powers.

We have pointed out that for bit rates lower than 1 Mbps, it is possible to follow the medical data of a mobile patient with low power emitted in IR, which is very important for reasons of consumption since the system is worn. In addition, we have shown that low VLC downlink powers are possible, which means that in this case the illumination is very little perceptible to the human eye. Some use cases such as the video application, however, require a high transmitted IR power, which does not satisfy the eye safety, or a too high level of illumination if the targeted performance is high. However, our study concludes that for most medical sensors, the proposed system is able to satisfy the overall quality in terms of bit error rate and data rate, which makes the proposed Li-Fi communication system promising.

Fig. 1 .

 1 Fig. 1. Environment description for the bidirectional full optical link

Half-power angle 𝜑 1

 1 conceptsIn order to study the IR and VLC channels, we used a raytracing based simulator developed in XLIM laboratory. It allows determining the impulse response h(t)

CaseFig. 2 .

 2 Fig. 2. Shape of the body for IR and VLC modeling and cases description

Fig. 3 .

 3 Fig. 3. 𝐻 0-𝐼𝑅 Probability Density Function and Cumulative Density Function

Fig. 4 .

 4 Fig. 4. 𝐻 0-𝑉𝐿𝐶 Probability Density Function and Cumulative Density Function

Fig. 5 .

 5 Fig. 5. 𝐻 0-𝑉𝐿𝐶 and 𝐻 0-𝐼𝑅 Probability Density Function for case (b)

  ) VLC channel case (b) IR channel

Table III .

 III Examples of performance for medical sensor applications

		ρb2=0.1		ρb2=0.9
	𝑃 𝑜𝑢𝑡	Emean (Lx)	𝑃 𝑡-𝐼𝑅 (mW)	Emean (Lx)	𝑃 𝑡-𝐼𝑅 (mW)
		Accelerometer, blood pressure, SPO2, temperature
		Rb <10kbps; BER<10 -10
	10 -1	0.48	4	0.3	2.6
	10 -2	1.4	20	0.45	13.5
	10 -3	2.4	25.3	0.55	18
		Electrocardiogram (ECG); Rb=72 kbps; BER<10 -10
	10 -1	1.3	12.6	0.9	7.6
	10 -2	3.7	56	1.2	38.3
	10 -3	6.5	67.8	1.4	50
		Electromyogram (EMG); Rb=1.5 Mbps; BER<10 -10
	10 -1	5.9	48	4.1	32
	10 -2	16.8	236	5.5	145
	10 -3	29.6	301	6.7	220
		Video; Rb=10 Mbps; BER<10 -3
	10 -1	7.8	55	5.4	38
	10 -2	22.2	331	7.2	223
	10 -3	39.2	430	8.9	293