Florian Monsef

Andrea Cozza

Limitations of the equivalence between spatial and ensemble estimators in the case of a single-tone excitation

Keywords: numbers: 43.55.Cs, 43.55.Br, 43.55.Gx

The ensemble-average value of the mean-square pressure is often assessed by using the spatialaverage technique, underlying an equivalence principle between spatial and ensemble estimators. Using the ideal-diffuse-field model, the accuracy of the spatial-average method has been studied theoretically forty years ago in the case of a single-tone excitation. This study is revisited in the present work on the basis of a more realistic description of the sound field accounting for a finite number of plane waves. The analysis of the spatial-average estimator is based on the study of its convergence rate. Using experimental data from practical examples, it is shown that the classical expression underestimates the estimator uncertainty even for frequencies greater than Schroeder's frequency, and that the number of plane waves may act as lower bound on the spatial-average estimator accuracy. The comparison of the convergence rate with an ensemble-estimator shows that the two statistics cannot be regarded as equivalent in a general case.

I. INTRODUCTION

In room acoustics, the interest of statistical approaches lies in the ability of predicting statistical properties of the sound field independently of the fine details of a specific configuration. Indeed, the configuration of a room and the position of the source impact on the modal topographies or equivalently on the directions and the amplitudes of the plane waves composing the sound field. If one considers an ensemble of rooms with approximately the same volume and the same reverberation time, one should recover the same estimated moments in accordance with the statistical model.

Statistical models have been developed for more than a half century. The most common one, called the idealdiffuse-field model, is based on random wave theory [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] . The model considers, on the one hand, the reverberant part of the sound field, that is the sound field measured at some distance from the source and, on the other hand, that the modal overlap is sufficiently great. It has been shown that a minimum average overlap of three modes is needed [START_REF] Schroeder | On frequency response curves in rooms. comparison of experimental, theoretical, and monte carlo results for the average frequency spacing between maxima[END_REF] , from which ensues Schroeder's frequency. Other studies [START_REF] Lyon | Statistical analysis of power injection and response in structures and rooms[END_REF][START_REF] Davy | The relative variance of the transmission function of a reverberation room[END_REF][START_REF] Weaver | On the ensemble variance of reverberation room transmission functions, the effect of spectral rigidity[END_REF] , based on a modal approach, aimed at deriving models that could be applied to a more general class of systems, i.e., to configurations where a low modal overlap can be assumed.

In practice, the validation of any model depends on estimation techniques based on the use of samples from which estimates of the statistical moments (generally the first two moments) must be derived. The random samples can be obtained mainly by two different methods. The first method consists in sampling the sound presa) florian.monsef@lss.supelec.fr sure in the central volume of a room [START_REF] Davy | The relative variance of the transmission function of a reverberation room[END_REF][START_REF] Waterhouse | Statistical properties of reverberant soundfields[END_REF][START_REF] Jacobsen | The ensemble variance of pure-tone measurements in reverberation rooms[END_REF] , for a given configuration of the sound-field topography. The estimation is based on a single realization of the sound-field topography. The estimation method using spatial samples will hereafter be referred to as the spatial-estimator method. The second method consists in using a single spatial sample for an ensemble of different sound-field topographies [START_REF] Balachandran | Random sound field in reverberation chambers[END_REF][START_REF] Jacobsen | The coherence of reverberant sound fields[END_REF] that can be obtained electronically, by the use of a varying tone signal [START_REF] Kuttruff | Room acoustics[END_REF] , or mechanically, using a large rotating vane [START_REF] Balachandran | Random sound field in reverberation chambers[END_REF][START_REF] Jacobsen | The coherence of reverberant sound fields[END_REF] . This latter is particularly well-suited for single-tone excitation which is the case of interest in the present study. Hereafter, this estimation method will be referred to as an ensemble-estimator method.

Using the squared value of the sound pressure related to each samples set, two estimations of the mean-square pressure can then be computed. Although never demonstrated, to the best of our knowledge, the estimates resulting from the spatial-and ensemble-average estimators are often regarded as equivalent. This equivalence principle implies a similar asymptotic convergence rate of both estimators. The equivalence assumption is directly related to the mixing (and consequently ergodic) property of dynamical systems [START_REF] Joyce | Sabine's reverberation time and ergodic auditoriums[END_REF] . Using a classical wave theory approach, an analogous mixing property is found within the ideal-diffuse-field description [START_REF] Witting | Modelling of diffuse sound field excitations and dynamic response analysis of leightweight structures[END_REF] . This equivalence assumption has been widely used to optimize spatialaveraging techniques, studied in detail by Schroeder [START_REF] Schroeder | Spatial averaging in a diffuse sound field and the equivalent number of independent measurements[END_REF][START_REF] Schroeder | Effect of frequency and space averaging on the transmission responses of multimode media[END_REF] , Lubman [START_REF] Lubman | Spatial averaging in a diffuse sound field[END_REF][START_REF] Lubman | Spatial averaging in sound-power measurements[END_REF][START_REF] Lubman | Effectiveness of continuous spatial averaging in a diffuse sound field[END_REF] and Waterhouse [START_REF] Waterhouse | Discrete versus continuous space averaging in a reverberant sound field[END_REF] forty years ago, to provide the effective number of discrete independent spatial samples for practical experiments. However, in the case of a single-tone excitation, these works are restricted to frequencies above Schroeder's frequency.

The aim of the present paper is to study theoretically, in the case of a harmonic excitation and considering only the reverberant part of the sound pressure, wether the spatial-and ensemble-average estimators can be regarded Spatial and Ensemble estimators 1 as equivalent. To this end, considering the mean-square pressure, expressions of the first two moments of the estimators are derived. The ensemble and spatial estimators are shown to be both unbiased but with different convergence rates. As a result, the confidence intervals related to both estimators cannot be regarded as identical.

To derive these results, a plane-wave model ensuing from a modal expansion is used to express the sound pressure. Recalling on the one hand that a mode can be expanded into a strictly finite number of plane waves [START_REF] Morse | Theoretical Acoustics[END_REF] , and considering on the other hand a finite number of overlapping modes, the model used includes a discrete number of plane waves. Hence, this approach ensures a realistic description of the field and is presented in section II. Using this sound-field model, we establish in section III a general expression of the spatial-average estimator variance. The accuracy of the estimation is assessed on the basis of its convergence rate, and it is shown to present a lower bound directly related to the number of plane waves composing the sound pressure. In section IV we relate the number of plane waves to the room parameters, providing a predictive model applied to practicalroom examples. We show that the lower bound of the convergence rate has a significant impact on the accuracy of the estimated value even for frequencies above Schroeder's frequency.

II. HARMONIC SOUND FIELD MODEL

A. From modal theory to discrete plane-wave expansion

In any enclosure with reflective boundaries (room, theater, reverberation chamber,...), a sound pressure induced by a source signal is made up of standing waves related to the normal modes that form a complete set of orthogonal functions [START_REF] Morse | Theoretical Acoustics[END_REF] . Accordingly , a harmonic sound field can be expressed as a sum of M dominant overlapping modes at a given frequency. Considering a harmonic excitation, the sound pressure will be expressed in the frequency domain. The resulting complex pressure p ω (r) can then be stated as

p ω (r) = M ∑ m=1 α m p m (r)ψ m (f ), (1)
where p m (r), α m and ψ m (f ) correspond to the modal topography, the modal weight and the frequency response of the m th -mode, respectively. As it will be presented thereafter, the spatial-average estimator consists basically in computing a spatial integral in which the modal topographies are needed. The computation is achievable if explicit expressions of the modes are available which is only the case in very specific configurations. This points out a serious limitation of the modal expansion in (1) as a way to study the accuracy of the spatial-average estimator.

Nevertheless, in the case of a standard room, a mode can be expanded into a finite number of plane waves. Consequently, for a finite modal overlap, the sound field can be rigourously expressed by means of a finite number of plane waves. To the best of our knowledge, this is the first analysis in room acoustics introducing such a model, allowing a more realistic description of the sound field. Indeed, it has been known for a long time that the infinite number of plane waves assumed in the ideal diffuse model is an unrealistic concept [START_REF] Lyon | Needed: a new definition of diffusion[END_REF] and is not physically justified. In practical terms, the finite number of plane waves accounts for the finite number of possible paths (as illustrated in Ref.10), or periodic orbits [START_REF] Chinnery | Experimental visualization of acoustic resonances within a stadium-shaped cavity[END_REF] , inducing cavity resonances and stationary waves. The resulting complex pressure p ω (r) expressed in (1) can then be restated as

p ω (r) = N ∑ p=1 γ p e -jk0 kp •r , ( 2
)
where N is the number of plane waves, r is the position at which the pressure is sampled and k 0 = ω/c is the propagation constant for the homogeneous medium filling the cavity, whit ω the angular frequency of the harmonic excitation signal and c the sound speed. The p th -plane wave propagates along the direction kp , with a complex amplitude γ p . We want to stress that the plane-wave model introduced in (2) should be regarded as a modal approach expressed in the Fourier basis. In other words, using (1) or ( 2) is totally equivalent. The main advantage of ( 2) over (1) however, is that the exponential term provides an explicit spatial dependence well suited to the study in section III of the spatial estimator. This will allow thereafter a rigorous derivation of the statistical moments of the spatial estimator.

B. Ensemble statistics

In practice, the plane-wave sets of parameters { k} = { k1 , . . . , kN } and {γ} = {γ 1 , . . . , γ N } cannot be measured, the only observable being the pressure. As commonly done in statistical mechanics and in the modeling of complex systems, the use of (2) allows studying sound field statistics in a simple way, as soon as the two plane-wave sets of parameters { k} and {γ} are treated as random variables.

Hence, to make the derivation of our approach simpler, we adopt the following assumptions on { k} and {γ} 1. the two sets will be assumed to be independent 2. for each set, variables identified by different indexes p in (2) will be regarded as independent and identically distributed (iid).

Although the expression of the sound pressure in (2) holds in a general case, i.e., for the sound field near the source and at some distance from the source, the hypotheses made on the sets { k} and {γ} implies physically a certain degree of diffusion of the sound field. Accordingly, the model that will be established in the present study is meaningful only for the reverberant part of the field with a certain degree of diffusion.

Spatial and Ensemble estimators 2

The mean-square pressure being the quantity of interest, we use (2) to express the square of the pressure as follows

p 2 ω (r) = N ∑ p=1 |γ p | 2 + N ∑ i̸ =j γ i γ * j e -jk0( ki -kj )•r , (3)
where * stands for the complex conjugate. In order to derive the statistical properties of (3), we start by expressing its ensemble-average value µ p 2 ω (r) defined as

µ p 2 ω (r) = E [ p 2 ω (r) ] = N ∑ p=1 E [ |γ p | 2 ] + N ∑ i̸ =j E[γ i γ * j ]E [ e -jk0( ki -kj )•r ] , (4)
where E[•] stands for the ensemble-average operator. It is noteworthy that we have taken advantage of the independence assumption between the sets {γ} and { k} to split the ensemble operator in the last term in (4). One of the properties sought in room acoustics applications, e.g., power measurements, is the statistical uniformity of the sound field. In other words, µ p 2 ω is expected to be independent of the position r. One way to justify this property in some previous work [START_REF] Schroeder | Effect of frequency and space averaging on the transmission responses of multimode media[END_REF] , is to consider {γ} as a set of zero mean complex random variables, such as

E [γ * i γ j ] = E [γ * i ] E [γ j ] = 0 ∀i ̸ = j. (5)
In fact, this assumption is not necessary, as long as the sampling position is in a region where the ensemble of the plane waves cover uniformly angles of arrival over 4π steradian. It is worth noting that this property holds even for a finite number of plane waves. Assuming the { k} to be iid leads to have the ensemble-average value of the complex exponential term in (4) to be equal to zero. Therefore,

µ p 2 ω (r) = N ∑ p=1 E [ |γ p | 2 ] =N E [ |γ p | 2 ] . (6)

III. SPATIAL ESTIMATOR

In order to characterize in detail the properties of this estimation method, we make use, for convenience, of a continuous spatial average which provides the following estimated value μp 2

ω μp 2 ω = 1 M Ω ∫ Ω p 2 ω (r) dr = ⟨p 2 ω (r)⟩ Ω , ( 7
)
where ⟨•⟩ Ω stands for the space averaging operator as taken over the sample region Ω.

In order to derive a more explicit expression of the spatial-average estimation, we substitute (3) into (7).

Accordingly,

μp 2 ω = 1 M Ω N ∑ i=1 N ∑ j=1 γ i γ * j ∫ Ω e -jk0( ki-kj )•r dr = N ∑ i=1 N ∑ j=1 γ i γ * j ⟨e -jk0( ki-kj )•r ⟩ Ω . ( 8
)
Organizing the set of the plane-wave weights {γ} into the γ column vector yields the following quadratic form μp 2 ω = γ H P γ, (9) where H stands for the Hermitian operator and P is a matrix with elements P ij defined as

P ij = ⟨e -jk0( ki-kj )•r ⟩ Ω . ( 10
)
In the case of a 3D cartesian sampling volume, the elements of P can be expressed as,

P ij = 1 M Ω ∫ Ω e -jk0( ki-kj )•(xn1+yn2+zn3) dx dy dz, ( 11
)
where n1 , n2 and n3 are the Cartesian unit vectors related to the x, y and z components, respectively. If we consider a square parallelepiped whose sides are referred to as a 1 , a 2 , a 3 , respectively, and if the origin is taken at the middle of the sampling volume Ω, we can easily restate (11) as

P ij = 3 ∏ v=1 sinc ( π a v λ ( ki -kj ) • nv ) , ( 12
)
where sinc (•) is the sine cardinal function, λ is the wavelength in the filling medium. We note that the elements of matrix P are a function of the random variables kikj .

Beyond the mathematical convenience provided by the concise form of (9), the matrix P is of particular interest in understanding the influence of the non-diagonal terms of P on the estimator accuracy. To illustrate our statement, suppose that P reduces to the identity matrix I. A perfect estimator is obtained since μp 2 ω converges to the ensemble average given in (6). The scenario ensuring this convergence is however unrealistic since corresponding to an infinite sampling volume as illustrated by relation (12). Indeed, if a v tends to infinity the cross-terms P ij tend to zero.

In practice however, the sampling volume being finite, non-diagonal terms of P will appear, i.e., P ̸ = I, inducing inevitably estimation errors. In other words, the analysis of the spatial-estimator accuracy lays on the analysis of the eigenvalues of P .

In order to study the accuracy of the estimate μp 2 ω , we need to compute the moments of the spatial estimator with respect to ensemble statistics.

A. Statistics of the error induced by a spatial estimator

Spatial and Ensemble estimators 3

The accuracy of the spatial estimator lies in the ability of the corresponding estimate to provide the true ensemble-average value. Accordingly, the quantity

ϵ = μp 2 ω -µ p 2 ω µ p 2 ω ( 13
)
provides a measure of the relative error of the spatialaverage value with respect to the ensemble-average value.

In order to assess the accuracy of the spatial estimator the two statistical moments of ( 13) must be studied. The first moment, i.e., the mean value, informs on the asymptotic convergence of the spatial estimator towards the ensemble-mean value, i.e., informs on the biased or unbiased nature of the estimator. The second moment, i.e., the variance, is crucial since it assesses the confidence interval of the estimated value, i.e., indicates the accuracy of the estimate. Before proceeding, we express the relative error ϵ in a different form for a more convenient mathematical derivation. We start by expressing ( 13) in a more explicit form. Using ( 6) and ( 9) yields

ϵ = μp 2 ω µ p 2 ω -1 = γ H P γ N E [|γ p | 2 ] -1. ( 14
)
It is convenient to introduce the random vector

v = γ √ E [|γ p | 2 ] , (15)
which is made up of iid normalized random variables, i.e., E

[ |v i | 2 ] = 1, ∀i ∈ [1, N ]. That allows recasting (14) in a simpler form ϵ = 1 N v H P v -1. ( 16
)
As explained in the previous subsection, a convenient way of analyzing the accuracy of the spatial estimator is to introduce the eigenvalues of P using the following change of basis

P = XΛX -1 , ( 17
)
where Λ is a diagonal matrix composed of the eigenvalues λ i of P . Consequently,

ϵ = 1 N w H Λw -1, (18)
where w = X -1 v is the vector of the plane-wave weights in the new basis. X being a unitary matrix, the probability density function of the random variables v i and w i is the same. Moreover, by computing the covariance matrix of the set {w}, it can easily be shown that the statistical properties assumed for the v i random variables are conserved for the w i .

B. Mean of the error

The statistical description of the relative error deals with a multivariate statistics analysis. To perform such analysis we need to consider a joint probability density function p

( {v}, { k}
) . Recalling the independence assumptions between the two sets of random parameters {v} and { k}, the joint probability function can be factored into the following way p ( {v}, { k}

) = p v ({v}) p k ( { k} ) , ( 19
)
where p v ({v}) and p k ( { k}

) are the marginal probability density function of the random sets {v} and { k}, respectively.

Hereafter, we will adopt the notation E v [ϵ] related to the conditional average of the random function ϵ conditioned by the random set {v}, defined as

E v [ϵ] = ∫ ϵ ( {v}, { k} ) p k ( { k} ) d{ k}. ( 20
)
The ensemble average µ ϵ of the error can then be derived using nested conditional averages applied to the two sets of concern in (18) as follows

µ ϵ = E [ϵ] = E k [E w [ϵ]] . ( 21
)
From expressing the average at the inner level, it follows that

E w [ϵ] = 1 N E w [ w H Λw ] -1 = 1 N N ∑ i=1 λ i E [ |w i | 2 ] -1. (22)
Reminding that X is a unitary matrix yields

E [ |w i | 2 ] = E [ |v i | 2 ] = 1 ∀i ∈ [1, N ] . ( 23
)
Consequently,

E w [ϵ] = 1 N N ∑ i=1 λ i -1. ( 24
)
The sum over the λ i is equal to Tr(P ), the trace of P . Besides, the main diagonal of P being made of ones, it follows that

Tr(P ) = N ∑ i=1 λ i = N, (25)
which implies in a clear way that

E [ϵ] = 0. ( 26
)
We can see that the estimator μp 2 ω is still unbiased for a finite number of plane waves, provided that these latter are sufficiently stirred to justify of a statistical approach that assumes the independence between the sets of random variables. The unbiased nature of the spatial estimator provides a necessary but not sufficient condition to consider spatial-averaging and ensemble-averaging measures as mean-ergodic processes [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF] . Indeed, one needs to study the variance of the error to assert the equivalence properties in a more general frame.

Spatial and Ensemble estimators 4

C. Variance of the error

The mean of the error being null according to (26), the variance σ 2 ϵ of the error reduces to E

[ ϵ 2 ]
. Adopting much the same approach as previously, the variance can be expressed using nested conditional averages as

σ 2 ϵ = E k [ E w [ ϵ 2 ]] . ( 27
)
From ( 18) we express ϵ 2 as follows

ϵ 2 = ( 1 N w H Λw -1 ) 2 = ( 1 N N ∑ i=1 λ i |w i | 2 ) 2 + 1 - 2 N N ∑ i=1 λ i |w i | 2 . ( 28
)
At the inner level of (27) the conditional average on the {w} set can be developed as

E w [ ϵ 2 ] = E w   ( 1 N N ∑ i=1 λ i |w i | 2 ) 2   + 1 - 2 N N ∑ i=1 λ i E w [ |w i | 2 . ] (29)
Recalling that the components w i are normalized random variables, the use of ( 25) yields the simpler form

E w [ ϵ 2 ] = 1 N 2 E w   ( N ∑ i=1 λ i |w i | 2 ) 2   -1. (30)
To obtain a more explicit expression, we expand the term between square brackets as

( N ∑ i=1 λ i |w i | 2 ) 2 = N ∑ i=1 λ 2 i |w i | 4 + ∑ i̸ =j λ i λ j |w i | 2 |w j | 2 . ( 31
)
This allows recasting (30) as

E w [ ϵ 2 ] = 1 N 2   ν 4 N ∑ i=1 λ 2 i + ν 2 2 ∑ i̸ =j λ i λ j   -1, (32)
where the ν n terms are the following conditional moments defined as,

ν n = E w [|w i | n ] = E [|w i | n ] . ( 33
)
In order to make (32) more tractable, we use algebraic properties stating that

N ∑ i=1 λ 2 i = Tr ( P T P ) = ||P || 2 F = N ∑ i=1 N ∑ j=1 P 2 ij , ( 34
)
where ||P || F stands for the Frobenius norm of matrix P Taking the square of (25)

(Tr(P )) 2 = N ∑ i=1 λ 2 i + N ∑ i̸ =j λ i λ j = N 2 , ( 35
)
and using (34), we can show that ∑ i̸ =j

λ i λ j = N 2 -||P || 2 F . (36) Whence, E w [ ϵ 2 ] = ||P || 2 F N 2 ( ν 4 -ν 2 2 ) . ( 37
)
In order to express the variance of the error, we still need to compute the conditional average on the { k} set, leading to

σ 2 ϵ =E k [ ||P || 2 F N 2 ( ν 4 -ν 2 2 ) ] = 1 N 2 ( ν 4 -ν 2 2 ) E k [ ||P || 2 F ] . ( 38
)
We can note that the conditional moment on { k} only applies on the square Frobenius norm since this latter is the only term in (37) depending on the plane-wave directions. Reminding the structure of P , we can restate (34) as

||P || 2 F = ∑ i,j P 2 ij = N + ∑ i̸ =j P 2 ij ( 39
)
The non-diagonal elements P ij , denoted by a function Φ( kikj ), inherit the iid property from the { k} random set. Accordingly, we can expand the conditional average

E k [ ||P || 2 F ] as E k [ ||P || 2 F ] = N + N (N -1)E k [ Φ 2 ( ki -kj ) ] . (40)
Finally, the variance of the error can be expressed as

σ ϵ 2 = [ N -1 N G Ω (a) + 1 N ] ( ν 4 -ν 2 2 ) (41)
with,

G Ω (a) = ∫ ∫ Ω Φ 2 ( ki -kj )p k( ki )p k( kj )d ki d kj , ( 42
)
which is the only term depending on the dimensions of Ω, where a is a vector gathering the sampling-area dimensions a v used in (12). Using a Monte-Carlo method and assuming a uniform angular distribution, the evolution of G Ω is plotted in Fig. 1, where a line of length a (dashed line), a square of side a (dotted line) and a cube of side a (solid line) have been considered, respectively. As shown in Appendix A, the term ν 4 -ν 2 2 in (41) is the relative-ensemble variance σ 2 r,p 2 of the mean-square pressure, defined as,

σ 2 r,p 2 = E [ ( p 2 ω (r) -µ p 2 ω (r) ) 2 µ 2 p 2 ω (r) ] . ( 43
)
If no spatial averaging is performed, then G Ω (0) = 1.

It follows that the variance σ ϵ 2 of the relative error reduces logically to σ 2 r,p 2 , i.e., to the relative variance of Spatial and Ensemble estimators 5

the mean-square pressure in a single point. The relativeensemble variance has been widely studied in the literature using wave theory [START_REF] Jacobsen | The ensemble variance of pure-tone measurements in reverberation rooms[END_REF][START_REF] Jacobsen | Sound power emitted by a Spatial and Ensemble estimators 11 pure-tone source in a reverberation room[END_REF] and modal theory 3-5 models and will not be analyzed in detail herein. All these models have been validated experimentally, and show that the relative-ensemble variance σ 2 r,p 2 is inversely proportional to the number of overlapping modes M s and tends to one for an ideal diffuse field, i.e, for an infinite number of plane waves.

Concerning spatial averaging, the modal models of Davy [START_REF] Davy | The relative variance of the transmission function of a reverberation room[END_REF] and Weaver 5 included the influence of multiple independent receivers, i.e., of independent spatial samples, in the overall relative variance. Accordingly, if we consider a single fixed source and N R independent receivers, the mean estimated value μp 2 ω of p 2 ω (r) can be deduced using an arithmetic mean, i.e., a discrete version of (7), as

μp 2 ω = 1 N R NR ∑ i=1 p 2 ω,i , (44)
which is the traditional mean-estimator. Using modal models [START_REF] Davy | The relative variance of the transmission function of a reverberation room[END_REF][START_REF] Weaver | On the ensemble variance of reverberation room transmission functions, the effect of spectral rigidity[END_REF] , it follows that

σ ϵ 2 = 1 N R σ 2 r,p 2 , ( 45
)
where 1/N R is referred to as the convergence rate. Within the modal approach, it appears that the more (independent) receivers we take the more accurate the estimated mean value will be. This of course unrealistic because a spatial correlation exists; but as explained in section II modal models cannot include easily the effect of spatial field distribution into the moments of spatialaverage estimates.

Using the plane-wave model presented herein, (41) can interestingly be restated as follows

σ ϵ 2 = G(a, N )σ 2 r,p 2 , ( 46
)
where G(a, N ) is the convergence rate of the spatial estimator expressed as

G(a, N ) = N -1 N G Ω (a) + G P W (N ), (47)
with

G P W (N ) = 1/N .
Note that on a single point, i.e., for a = 0 the convergence rate G(0, N ) = 1 as expected, independently of the number N of plane waves.

It turns out that the convergence rate of the spatial estimator obtained in (47) is conditioned by two terms. The first term, as expected, is related to the influence of the sampling term G Ω (a) whose inverse corresponds to an equivalent number of independent samples. The number of plane waves has very little influence on the first term of (47). On the contrary, the second term depends only on the number of plane waves.

The influence of the number of plane waves in the spatial-averaging convergence rate has never been predicted to the best of our knowledge by any study. The study of the accuracy of the spatial-average estimator consists then in analyzing the convergence rate given in (47) which can be simplified as,

G(a, N ) ≃ G Ω (a) + G P W (N ), (48)
since the term (N -1)/N can be reasonably neglected. Accordingly, to study which of the two terms dominates in the overall estimator variance, the knowledge of σ 2 r,p 2 is not needed since this latter is common. Consequently, we shall focus our analysis in the next section, on a comparison of the two terms of (48).

IV. ACCURACY OF THE SPATIAL ESTIMATOR

A. Existence of a lower bound

We are interested in comparing the two terms in (48) in order to see if the contribution of the plane waves may be significant in some cases. This is of relevant importance because the presence of the 1/N term in the variance expression acts as a lower bound in the spatialaverage estimator accuracy. Indeed, the number of plane waves is related to the room configuration, but is independent of the sampling area Ω. Consequently, if the sampling area is larger, the G Ω term will decrease, as seen in Fig. 1, but leaving the 1/N term unchanged. In practical terms, raising the sampling area will not diminish continuously the degree of uncertainty deduced from the spatial-average-estimator variance.

To stress the significative risk of this limitation, we study in the next subsection practical room cases for which an assessment of the two contributing terms, i.e., of G P W and G Ω , is performed.

Spatial and Ensemble estimators 6

B. Influence of room parameters

In order to assess quantitatively the role of G P W , we consider real room configurations. We use the experimental data presented in Ref. 7 for the case of a small room (40 m 3 ) and a large reverberation room (245 m 3 ) with large fixed diffusers. Both rooms are essentially rectangular. For the small room, lightly damped and damped cases are considered. The corresponding Schroeder's frequencies are 500 Hz and 330 Hz, respectively [START_REF] Jacobsen | The ensemble variance of pure-tone measurements in reverberation rooms[END_REF] . For the large room very little damping and large damping cases have been considered. The corresponding Schroeder's frequency are of 310 Hz and 200 Hz, respectively [START_REF] Jacobsen | The ensemble variance of pure-tone measurements in reverberation rooms[END_REF] . The measurements of the reverberation time T 60 were performed on a frequency range [100 Hz : 3.2 kHz] and is shown in Fig. 2.

The configuration of a room has a direct impact on the modal topographies but also on the frequency response and consequently on the number of plane waves composing the sound field. In order to fruitfully use the experimental room parameters we need to establish a way to relate them to the number of plane waves in the following way. At a frequency f , we may consider β plane waves per resonant mode. In order to assess the number M of modes that are found within the average resonant -3dB bandwidth B M , we use the simplest form of Weyl's approximation [START_REF] Weyl | Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Differentialgleichungen (mit einer Anwendung auf die Theorie der Hohlraumstrahlung)(Asymptotic distribution of the partial differential equations eigenvalues (with an application to black body theory))[END_REF] of the modal density. The number of plane waves involved in the sound field can then be approximated as

N = β • M = β n(f )B M ≃ β 4πV f 2 c 3 B M , ( 49
)
where n(f ) is the modal density (in modes per Hertz), V is the volume of the room and c the sound speed.

Recalling that the reverberation time T 60 is related to the average bandwidth of a mode B M , such as B M = 6.91/π • T 60 yields the following expression of the plane-200 300 400 500 600 700 800 900 10 wave term in ( 48)

G P W = c 3 T 60 27.64 β V f 2 , ( 50
)
providing a relation allowing to assess the lower-bound term using the room reverberation time T 60 and the room volume V . For the four room configurations, Fig. 3 shows the evolution of G P W as a function of frequency (lower x -axis), and the evolution of the spatial G Ω term as a function of a/λ in the 2D case for 3 < a/λ < 10 (upper x -axis) where a value of 3.3 meters has been assumed for the parameter a. As already said, the number of plane waves per mode is hard to estimate. The rooms being essentially rectangular we can expect approximately eight plane waves per mode. In order to consider an upper bound we use a value of ten plane waves per mode. We stress that this value being used for all the modes on average, the G P W term may not be regarded as overvalued.

For the case of the small room Fig. 3a shows that the plane-wave term dominates even for frequencies above Schroeder's frequency, whereas for the large room cases described in Fig. 3b, the two terms are of the same order of magnitude. These results can be regarded as a further proof to the well-known fact that the acoustics of small rooms are in a way more complicated than those of large ones. In the present case, it is an illustration that, at a given frequency, the smaller the room the worse the accuracy of the spatial estimate will be; not only because of the relative-ensemble variance as commonly known, but Spatial and Ensemble estimators 7 4. GP W (bold line) and GΩ (solid line with circles) as a function of the acoustical volume V λ of a room. GΩ is related to the 2D-section of the sampling volume for which two cases are considered : the sampling volume fits the entire room volume (α = 1) and only 10% of the room volume is used in the second case (α = 0.1). For GP W , frequencies below and above Schroeder's frequencies have been chosen (100Hz and 1kHz, respectively) and three values of reverberation times given in Fig. 2 were used, namely the minimum value of 1 s, the maximum value of about 6 s and an in-between value of 3 s. also due to the spatial-averaging effect assessed herein. In both cases of Fig. 3, the finite number of plane waves has an influence on the overall convergence rate not only for frequencies near Schroeder's limit, as intuitively expectable, but in fact for a frequency range whose upper limit would be approximately ten times Schroeder's frequency.

At this stage a comparison between G P W and G Ω (a) has not been performed on the basis of a common parameter. To this end, we consider that the 3D-sampling zone, referred to as Ω 3D , is a portion α of the room volume V , such as Ω 3D = αV . In order to include the influence of the frequency (or equivalently the wavelength λ) at which the measurement is performed, we introduce the acoustical sampling volume defined as Ω 3D λ = Ω 3D /λ 3 . Equivalently, we define the acoustical volume of the room defined as V λ = V /λ 3 , yielding

G P W = T 60 f 27.64 β V λ . ( 51
)
The computation of G Ω requires the dimensions of the rooms. These latter being not known, we consider the following geometrical ratios 3 : 4 : 5, related to the width a 1 , the height a 2 and the length a 3 of the sampling volume, respectively. These ratios account for the absence in practice of degenerate eigenfrequencies.

For the sake of plotting clarity, we restrict the analysis to the a 2 a 3 section, referred to as Ω 2D , of the sampling volume Ω 3D . To compute the resulting G Ω term, two cases of Ω 3D are considered : the case for which the sampling volume fits the total room volume, i.e., for α = 1, and the case for which the sampling volume repre-sents 10% of the room volume, i.e., for α = 0.1. Accordingly, Fig. 4 shows the resulting 2D-G Ω term for α = 1 and α = 0.1, and the plane-wave term. For this latter, frequencies below and above Schroeder's frequency were considered, 100 Hz and 1 kHz, respectively, and three values of reverberation times were chosen on the basis of the experimental data provided by Fig. 2. The selected values correspond to the minimum of 1 s, to the maximum of about 6 s and to an in-between value of 3 s well suited to the large damped room and small lightly damped room cases.

We can see in Fig. 4 that at a fixed frequency (1kHz in the present case), the decrease of the reverberation time induces a lower convergence rate as stated by (51). If the frequency f varies, so does the acoustical volume since this latter is a function proportional to the cube of f . If we refer to the large and small rooms cases, taken at 1 kHz, respective acoustical volumes of 6000 and 991 are found. Using Fig. 4 for α = 0.1, we find that G P W = G Ω = 10 -2 for the small room; for the large room, G P W = 2 10 -3 and G Ω = 3.5 10 -3 . This illustrates the influence of the plane-wave term, since it contributes, for the small room case, to 50% of the overall convergence rate. TABLE I. For sampling volumes corresponding to 10% (α = 0.1) and 100% (α = 1) of the small-and large-room volumes, respectively, relative margin of errors mr,ϵ and the corresponding 2D sampling areas Ω2D are reported. The number Neq of equivalent independent spatial samples is also reported for α = 0.1, assumed as a more realistic case.

Configuration

mr,ϵ Ω2D Neq

small room α = 1 20% 31m 2 - α = 0.1 28% 6.8m 2 100 large room α = 1 10% 105m 2 - α = 0.1 15% 22m 2 250
In practice, the spatial-average estimate expressed in (7) can be regarded as a random Gaussian variable whose relative standard deviation σ ϵ conditions the confidence interval, i.e., the accuracy of the estimator. We define the margin of error m ϵ value as the half-width of the confidence interval. If we consider that 95% of the cases fall into this latter, the margin of error is approximately m ϵ = 2σ ϵ . Our concern dealing with the comparison of the two terms of (48), the estimator efficiency can be described by the following relative margin of error m r,ϵ = m ϵ /σ r,p 2 = 2 √

G Ω + G P W , plotted in Fig. 5 as a function of the acoustical volume of the room, for the 1D (triangles), 2D (circles) and 3D (squares) cases. For the plane-wave term a reverberation time of 3 s at 1 kHz has been used. We consider the cases for α = 1 (solid line) and α = 0.1 (dashed line) in order to assess the influence of the sampling area on the uncertainty range related to the spatial estimator efficiency. Spatial and Ensemble estimators 8 If we refer to the practical cases described previously for which G P W = G Ω = 10 -2 for the small room and G P W = 2 10 -3 and G Ω = 3.5 10 -3 for the large room, we find the respective relative margins of error of 20% and 10% for α = 1, and 28% and 15% for α = 0.1.

Regarding the G Ω term as an inverse number N eq of equivalent independent discrete spatial samples, 100 and 250 equivalent samples are found for the small and large room, respectively, for α = 0.1. All these values are summed up in Tab. I where the equivalent sampling 2Dareas Ω 2D have also been reported providing the order of magnitude of the surface related to the relative margins of error previously mentioned.

Interestingly, if the α = 0.1 case of the spatial estimator were replaced by an ensemble estimator (using a rotating vane or moving the source and receiver positions as in Ref.7) based on N s independent samples, we would need, on the basis of the same margins of error, i.e., 28% and 15%, 51 and 177 samples for the small and large room, respectively. For the small room case, to compensate the plane-wave term G P W , a spatial-averaging process needs twice more independent samples than an ensemble estimator. It is worth stressing that in the case of an ensemble estimator, each realization corresponds to the case of G(0, N ) = 1, i.e., the number of plane waves has no influence.

The method presented herein can be then extended in a general way to any configuration. It provides a practical theoretical tool for predictions and estimation procedures applied to the mean-square pressure in room acoustics.

In the general model derived in section III, we did not make any assumption on the degree of diffusion of the field unlike the works dedicated to spatial averaging in the litterature [START_REF] Schroeder | Spatial averaging in a diffuse sound field and the equivalent number of independent measurements[END_REF][START_REF] Lubman | Spatial averaging in a diffuse sound field[END_REF][START_REF] Lubman | Spatial averaging in sound-power measurements[END_REF][START_REF] Lubman | Effectiveness of continuous spatial averaging in a diffuse sound field[END_REF][START_REF] Waterhouse | Discrete versus continuous space averaging in a reverberant sound field[END_REF] . Indeed, these latter have been carried out within the frame of the ideal diffuse field hypoth-esis. To the best of our knowledge, no other studies have been presented in a more general case. Accordingly, it is interesting to analyze how the general expression, given in (46), becomes in the particular case of an ideal diffuse field.

C. The ideal diffuse field case

Relative variance expression

In the case of an ideal diffuse field, i.e., for an infinite number of plane waves, using (47), we can show that the convergence rate G(a, N ) reduces to G Ω (a). The relative-error variance of the spatial estimator given in (46) can then be expressed as lim

N →∞ σ ϵ 2 = G Ω (a) σ 2 r,p 2 . ( 52
)
The relative-ensemble variance of the mean square pressure σ 2 r,p 2 can be analyzed in a straightforward way. Indeed, the complex plane-waves weights {w} in the new basis result from an infinite sum of the {γ} random variables. Accordingly, using the central limit theorem, we can state that the real and imaginary parts of the set {w} follow a normal distribution. Consequently, the elements |w i | 2 follow a chi-square distribution with two degrees of freedom whose moments are such that

E[|w i | 4 ] E[|w i | 2 ] 2 = ν 4 ν 2 2 = 2. ( 53
)
Hence, recalling that E[|w i | 2 ] = 1, the relative variance σ 2 r,p 2 reduces to unity as commonly accepted for the ideal diffuse field case.

We can then express (52) in the following simpler form lim

N →∞ σ ϵ 2 = G Ω (a), (54)
which is the expression given by Schroeder [START_REF] Schroeder | Spatial averaging in a diffuse sound field and the equivalent number of independent measurements[END_REF] . In Schroeder's work [START_REF] Schroeder | Spatial averaging in a diffuse sound field and the equivalent number of independent measurements[END_REF] a Gaussian assumption is made on the {γ} set distribution. Note that, in the model derived in section III such assumption is not necessary to recover the expression of the error variance in the case of an ideal diffuse field.

Estimators equivalence

In order to analyze if a spatial-and ensemble-average estimators are equivalent, one needs to compare their convergence rates.

To this end, we introduce the ensemble-meanestimated value of the mean-square pressure µ p 2 ω defined as,

µ p 2 ω = 1 M r Mr ∑ i=1 p 2 ω,i , (55)
where p 2 ω,i is the mean-square pressure corresponding to the i-th realization of the ensemble statistics and M r is Spatial and Ensemble estimators 9 the number of independent realizations. A realization can either be a position of a rotating vane or a pair of source/receiver positions. If the total number M r of realizations are independent, µ p 2 ω is a random variable whose centered and normalized form, as done in (13), follows a normal distribution of zero mean and of variance σ 2 r,p 2 /M r . The spatial-average (centered and normalized) estimate ϵ follows also a normal distribution with a zero mean. However, the variance and/or convergence rate is not identical in all the spatial-sampling configurations as shown in Fig. 1. Indeed, we observe different slopes preventing a general equivalence principle between spatialand ensemble-average estimates as the confidence interval is not identical. Accordingly, the claim in Ref.17 that the spatial-average estimator provides the "true estimate" is not correct in a strict sense although the 1D/2D/3Dspatial-variances G Ω tend asymptotically to zero.

In order to carry on a more accurate comparison of the convergence rates of μp 2 ω and µ p 2 ω , we focus on the different slopes of G Ω for the 1D, 2D and 3D cases which turn to be 0.5, 0.4 and 0.26 per octave, respectively. For the 1D case, the slope is in accordance with the asymptotic expansion of the expression of Ref. 13. For the 2D-square and 3D-cube cases no analytical expressions being available, we have estimated the following asymptotic expressions

G Ω (a) ≃      0.5/(a/λ) in the 1D case 0.2133/(a/λ) 1.8 in the 2D case 0.0740(a/λ) -2 in the 3D case (56) valid for a/λ > 2.
These expressions, derived for the first time for the 2D and 3D cases to the best of our knowledge, allow to see the impact of the dimensionality of Ω on G Ω slopes.

The difference of slopes, linked to the difference of exponent power value in the asymptotic expressions of G Ω , is induced by the spatial correlation of the field in the room. The correlation effects on the spatial variance have already been studied in the cases of a line, a circle and a disk sampling areas [START_REF] Lubman | Spatial averaging in a diffuse sound field[END_REF][START_REF] Lubman | Spatial averaging in sound-power measurements[END_REF][START_REF] Lubman | Effectiveness of continuous spatial averaging in a diffuse sound field[END_REF] . For the mean-square pressure, the spatial autocorrelation function R(k 0 ∆r) was shown to be [START_REF] Lubman | Effectiveness of continuous spatial averaging in a diffuse sound field[END_REF][START_REF] Chu | Note on the independent sampling of meansquare pressure in reverberant sound fields[END_REF] ,

R(k 0 ∆r) = ( sin(k 0 ∆r) k 0 ∆r ) 2 , ( 57
)
where ∆r is the distance between two points. Using (57), a zero correlation is found for distances of λ/2, or entire multiples. On the basis of this property, the ensuing idea of using a discrete, rather than a continuous, spatial-average estimator has been studied by Waterhouse 18 who was interested in studying the different possible arrangements in 1D, 2D and 3D to have uncorrelated discrete samples. It can be shown that the 1D equivalent discrete averaging estimator, i.e., sampling along a line on positions apart from λ/2 (or multiples), is the only case for which uncorrelated samples can be obtained. Indeed, in the 2D and 3D cases no geometrical configuration can provide a zero correlation between each pair of points. Consequently, in the 1D case, doubling the number of (independent) samples, i.e. doubling a/λ, will reduce the estimator variance σ 2 ϵ of a factor 2, confirming the slope obtained in the 1D case in Fig. 1.

To compare G Ω with the convergence rate 1/M r of the ensemble estimator, we need to use a discrete spatialaverage estimator in order to have the number samples as a common parameter.

It has been shown [START_REF] Waterhouse | Discrete versus continuous space averaging in a reverberant sound field[END_REF] that the relative-error variance of the discrete spatial-average estimator can be fairly approximated by G Ω , as long as the geometric parameter a is greater than 3/2 λ. Accordingly, assuming a sufficient number of samples allows using (56) to describe properly the asymptotic convergence rate of the discrete spatial-average estimator. In order to compare the convergence rates of the two estimators, we need to introduce a number of spatial samples whose arrangement will impact inevitably on the accuracy of the spatial-average estimate.

We choose a regular arrangement of the samples separated by a distance λ/2. This arrangement allows, on the one hand, to have samples easily placed and on the other hand, to have a zero-correlation in many cases. If N s is the number of samples in the 1D case ensuing from the regular arrangement, N 2 s and N 3 s will be the samples obtained in the 2D-square and 3D-cube cases, respectively. Moreover, if a length a is assumed for the 1D case, the integer part of 1 + a/(λ/2) correspond to the number of samples. Equivalently, for a sufficient value of a, we have approximately 2a/λ samples along a line or along the side of a square and a cube of length a. Consequently, using the asymptotic approximations of G Ω we can compute the following asymptotic convergence rates Using the results of G Ω shown in Fig. 1, we plot in Fig. 6 the variance of the relative error of both estimators as a function of the number of samples.

G Ω (N s ) ∝      N -1
The 1D case and the ensemble estimator being based on independent samples, they converge with the same rate. Whereas, the 2D and 3D cases are affected by spatial correlation. Consequently, their convergence rates are lower. Accordingly, for the scenarios where the ideal diffuse field model can be assumed as reliable, the ofteninvoked equivalence assumption between spatial-average (in the 2D and 3D cases) and ensemble-average estimators does not hold.

V. CONCLUSION

The present study showed that the spatial-and ensemble-average estimators cannot be regarded as equivalent, since their convergence rates differ. The results were obtained using a wave-theory approach based on a finite number of plane waves; the plane-wave model is equivalent to a modal expansion, while allowing to explicitly express the spatial dependence of the sound-field distribution, well-suited to the study of the spatial estimator.

For an infinite number of plane waves the spatialaveraging estimates coincide with the results presented forty years ago using an ideal-diffuse field model. In the general case, the present work points out theoretically that the number of plane waves induce a limitation of the spatial-estimator accuracy, emphasizing the uselessness of over-sampling. This effect has never been predicted to the best of our knowledge. On the basis of practical cases, the results show that this effect impacts the convergence rate for frequencies about ten time Schroeder's frequency.

Although the model is physically justified and theoretically reliable, it would be interesting to establish this limitation on the basis of experimental measurements. These studies would provide an experimental evidence of the maximum accuracy that can be expected.

FIG. 1 .

 1 FIG.1. Evolution of the spatial term GΩ(a) in (42). 1D, 2D and 3D sampling domains have been considered, related to a line stretch of length a (dashed line), to a square (dotted line) and a cube of side a (solid line), respectively.

FIG. 2 .

 2 FIG.2. Reverberation times from Ref.7. Four configurations have been considered : a large room (solid line), a large damped room (dotted line), a small room lightly damped (dashed line) and a small damped room (dashed-dotted line).

FIG. 3 .

 3 FIG. 3. Comparison of the lower bound GP W as a function of frequency (lower x -axis) (with ten plane waves per mode) with the 2D-square spatial term GΩ as a function of a/λ (upper x -axis) for (a) the small room configurations and (b) the large room configurations of Ref. 7. A value of 3.3 m was considered for the parameter a.

  FIG.4. GP W (bold line) and GΩ (solid line with circles) as a function of the acoustical volume V λ of a room. GΩ is related to the 2D-section of the sampling volume for which two cases are considered : the sampling volume fits the entire room volume (α = 1) and only 10% of the room volume is used in the second case (α = 0.1). For GP W , frequencies below and above Schroeder's frequencies have been chosen (100Hz and 1kHz, respectively) and three values of reverberation times given in Fig.2were used, namely the minimum value of 1 s, the maximum value of about 6 s and an in-between value of 3 s.

FIG. 5 .

 5 FIG.5. The relative margin of error related to (48) for a 95% confidence interval and a reverberation time of about 3 s at 1 kHz. Sampling volumes in 1D (triangles), 2D (circles) and 3D (squares) are considered. The case for which the sampling volume Ω fits the entire room volume (α=1) is plotted in solid line -the case for which the sampling volume is 10% of the room volume is plotted in dashed line.

FIG. 6 .

 6 FIG.6. Variance of the spatial-average estimator in the 1D (triangles), 2D (circles) and 3D (squares) cases with respect to the ensemble-mean estimator (solid line with no markers), as a function of the discrete number of samples.

s

  in the 1D case N -0.9 s in the 2D case N -0.666 s in the 3D case (58) Spatial and Ensemble estimators 10

Spatial and Ensemble estimators 12

APPENDIX A: IDENTIFICATION OF THE VARIANCE

In order to identify the relative-ensemble variance σ 2 r,p 2 of the mean-square pressure defined as,

we must consider the error expression in (13) in a single point. The matrix P is then made of ones. Consequently, the matrix being of unitary rank, we have only a single non-zero eigenvalue λ 1 such as,

Accordingly, (32) can be expressed with the following simpler form,

Recalling that ν 2 2 = 1, the final variance is obtained computing the conditional average on the { k} set. However, ν 4 is independent of the { k} set, finally yielding σ 2 r,p 2 = ν 4 -ν 2 2 .

(A4)


Parameterize a territorial risk evaluation scale using multiple experts knowledge through risk assessment examples

Olivier Cailloux & Vincent Mousseau

Laboratoire Génie Industriel, École Centrale Paris, Châtenay-Malabry, France ABSTRACT: Evaluating and comparing the threats and vulnerabilities associated with territorial zones according to multiple criteria (industrial activity, population, etc.) can be a time-consuming task and often requires the participation of several experts and decision makers. Rather than a direct evaluation of these zones, building a risk evaluation scale and using it in a formal procedure permits to automate the assessment and therefore to apply it in a repeated way and in large-scale contexts and, provided the chosen procedure and scale are accepted, to make it objective. One of the main difficulty of building such a formal evaluation procedure is to account for the multiple experts knowledge and decision makers preferences. The procedure used in this article, ELECTRE TRI, uses the performances of each territorial zone on multiple criteria, together with preferential parameters, to qualitatively assess their associated risk level. The preferential parameters to be determined are the category limits, i.e. the limits on each criterion of the performance range associated with a given risk level, and the weights associated with the different criteria. To obtain these parameters with no direct questioning of the stakeholders, tools based on mathematical programming have been developed to deduce these preferential parameters from assessment examples. In this article, several such tools are applied together in order to build a complete example leading to the definition of a territorial risk evaluation scale taking into account the preferences of multiple stakeholders.

RISK ASSESSMENT MODELLED AS MULTICRITERIA SORTING MODELS

Assessing the risk levels associated with geographical zones can be a difficult task. It involves multiple, and often conflicting, point of views: a zone may have a low risk according to one criterion while being exposed to a critical risk according to an other one.

Associating one risk level to a zone involves aggregating these point of views. This article suggests to use the tools developed in the domain of multicriteria decision aid, which enable a formal approach to that aggregation problem when assessing risk. The paper is structured as follows. This section describes the relation between a risk assessment problem and a multicriteria sorting problem. Section 2 describes a specific multicriteria sorting method, ELEC-TRE TRI [START_REF] Figueira | ELEC-TRE methods[END_REF][START_REF] Mousseau | A user-oriented implementation of the ELECTRE TRI method integrating preference elicitation support[END_REF][START_REF] Roy | The outranking approach and the foundations of ELECTRE methods[END_REF], and discusses existing methodologies to elicit an ELECTRE TRI based risk model with experts. Section 3 suggests a strategy to build a risk evaluation model as an ELECTRE TRI model using indirect parameters elicitation. An illustrative example is given in Section 4. The algorithm used for the inference is detailed in Section 5.

MCDA problem statements

Real-world decision problems such as territorial risk evaluation usually involve multiple points of view on evaluation of zones they concern. These points of view are relevant for a Decision Maker (DM), or expert (who is either a single person or a collegial body). They are formally represented by real-valued functions called criteria. Multicriteria (MC) decision aiding aims at recommending a decision which is consistent with a value system of the DM. Various methodologies have been proposed to support DMs facing a MC decision problem [START_REF] Keeney | Decisions with multiple objectives: Preferences and value tradeoffs[END_REF][START_REF] Roy | Multicriteria Methodology for Decision Aiding[END_REF][START_REF] Bouyssou | Evaluation and decision models with multiple criteria: Stepping stones for the analyst[END_REF]. Three main problematics structure the MC decision aiding field: choice, ranking or sorting [START_REF] Tsoukiàs | On the concept of decision aiding process: an operational perspective[END_REF]. Given a finite set A of alternatives, the choice problematic consists in selecting a subset of alternatives (as small as possible) being judged as the most satisfactory. The ranking problematic consists in establishing a preference pre-order (either partial or complete) on the set of alternatives. These two problematics are said to be comparative, as they require to compare alternatives one to another in order to set up the choice set or the preference pre-order. The sorting problem-atic concerns ordinal classification of alternatives and consists in assigning each alternative to one of some pre-defined and ordered categories. The assignment of an alternative to an appropriate category relies on the alternative's intrinsic value (and not on the comparison to others). In this paper, we consider the multicriteria (MC) sorting problematic to represent qualitative risk assessment models.

Defining a qualitative risk assessment scale

Most concepts used in a risk assessment problem easily map to concepts used in the MC sorting problematic. What is called a category in such a problematic represents the set of possible risk levels that the territorial zones must be mapped into. They can be e.g.: {Critical risk, Medium risk, Low risk}. Ideally these risk levels should be defined according to some precautionary measures associated with them, to give these categories a precise meaning. The point of views involved in a risk evaluation problem correspond to the criteria in a MC sorting problem. Examples include the presence of a school or the percentage of vulnerable persons.

Building an evaluation scale thus amounts to build a MC sorting model.

Building a sorting model

Each zone being described by a vector of risk factors associated with the point of views involved in the problem, the task at hand consists in assigning these zones to a set of risk categories. This is done using a MC sorting preference model. This model contains a set of subjective data representing the preferences of the considered DM with regard to, e.g., the relative importance of each of the criteria. These objective and subjective data together with a sorting method allow to aggregate the different point of views to assess the risk level of each considered zone.

These preferential parameters may be elicited in a direct way, but this is often difficult as it requires the DM to undertand the fine details of their use in the considered MC sorting method. That is why it has been suggested to deduce the preferential parameters in an inverse way, by asking the DM examples of alternatives, or zones, and the category, or risk level, they would consider appropriate for these.

A supplementary difficulty arises when the evaluation method to be defined involve multiple DMs, as different stakeholders may favor different subjective parameter values. Applying inverse elicitation in a multiple DMs context amounts to ask each DM for a set of examples, which may be conflicting, and deduce preferential parameter values that may be either entirely shared by the DMs, or shared for a part of the parameters, and individual for other values. This is the approach we use in this article.

ELECTRE TRI FOR RISK EVALUATION

The MC sorting method used here is a simplified version of ELECTRE TRI. It is appropriate for a risk evaluation setting as it only requires ordinally evaluated performances on the different criteria. The version considered here is very close to the version studied by Bouyssou andMarchant (2007a, 2007b).

2.1 Sorting procedure ELECTRE TRI requires, as a definition of the preferences of a DM, criteria importance parameters and category limits separating the categories. The criteria importance parameters include a weight for each of the criteria and a majority threshold that defines when a set of criteria is good enough to be decisive. The category limits separate, for each criterion, two consecutive risk levels.

Consider a finite set of territorial zones A, a set of category limits B = {b 0 , . . . b k }, and a finite set of criteria indices J. A criterion g j (j ∈ J) is a function from A ∪ B to R where g j (a) denotes the performance of the zone a on the criterion g j . The zones have to be sorted in k risk levels, c 1 , . . . , c k , ordered by their desirability. c 1 is the worst (i.e. highest) risk level, and c k is the best (the lowest) one. Each risk level c h is defined by the performances of its lower frontier, or category limit, b h-1 and its upper frontier b h of B. The performances are here supposed to be such that a higher value denotes a better performance (i.e. associated with less risks) and the performances on the frontiers are non-decreasing, i.e. ∀j ∈ J,

1 ≤ h ≤ k : g j (b h-1 ) ≤ g j (b h ).
To sort the zones, ELECTRE TRI uses the concept of outranking relation. The assignment rule used here, known as the pessimistic rule, assigns a zone a to the highest possible risk level c h such that the zone outranks the category's lower frontier b h-1 .A zone a outranks a frontier b h-1 if and only if there is a sufficient coalition of criteria supporting the assertion "a is at least as good as b h-1 ", and no criterion strongly opposes (vetoes) that assertion. To compute this, preferential parameters given by a DM are used.The coalition of criteria in favor of the outranking, ∀a ∈ A, 1 ≤ h ≤ k, is defined as

j∈J w j C j (a, b h-1 ), (1)
where w j is the weight of the criterion g j , and C j (a, b h-1 ) ∈ {0, 1} measures whether a is at least as good as b h-1 from the point of view of the criterion j:

C j (a, b h-1 ) = 1 ⇔ g j (a) ≥ g j (b h-1
), 0 otherwise. The weights are defined so that they sum to one ( j∈J w j = 1). The coalition is compared to a majority threshold λ ∈ [0.5, 1] defined by the decision maker along with the weights. If j∈J w j C j (a, b h-1 ) < λ, the coalition is not a sufficient coalition and the zone does not outrank the frontier b h-1 and will therefore be assigned in a risk level below c h .

Even when the coalition is strong enough, a criterion may veto the outranking situation. It happens when g j (a) > v h-1 j . The veto threshold v h-1 j is a value that the DM may define and represents the performance that, if not reached by some zone a, forbids the zone to have a risk label of c h . To summarize, the zone a outranks the frontier b h-1 (and therefore is assigned to at least the category c h ) if and only if

j∈J w j C j (a, b h-1 ) ≥ λ and ∀j ∈ J : g j (a) > v h-1 j .
In a case involving a single DM, the weights and majority thresholds (defining the sufficient coalitions) and the category limits may be given directly by him. However, this requires that the DM understands how these values will be used. It is moreover a difficult process to directly ask the DM for these parameters. The approach used here supposes that he provides assignment examples which are used to infer the preferential parameters.

The situation is even more complex when several DMs are involved: there is no reason to suppose that all DMs a priori agree on the importance of the criteria or on the frontiers parameters. Hence, in addition to the difficulty related to the inference of an ELEC-TRE TRI model, an additional level of complexity is introduced by a possible divergence in the preferences expressed by the different DMs.

Inferrence of preferential parameters in a multiple DM context

The ELECTRE TRI preferential parameters to elicit are the category limits, the weights, and the vetoes. Previous works aiming to infer preferential parameters for the ELECTRE TRI procedure on the basis of assignment examples exist, but they mostly involve a single DM. Existing approaches suggest to find the entire ELECTRE TRI preference model parameters [START_REF] Mousseau | Inferring an ELEC-TRE TRI model from assignment examples[END_REF] from assignment examples, or find the importance coefficients only [START_REF] Mousseau | Using assignment examples to infer weights for ELEC-TRE TRI method: Some experimental results[END_REF], or the categories limits (Ngo The and Mousseau 2002), the other parameters being supposedly known. Robust approaches are suggested which compute for each alternative a range of possible categories to which alternatives can be assigned under incomplete determination of the parameters [START_REF] Dias | On computing ELEC-TRE's credibility indices under partial information[END_REF][START_REF] Dias | ELECTRE TRI for groups with imprecise information on parameter values[END_REF][START_REF] Dias | An aggregation/disaggregation approach to obtain robust conclusions with ELECTRE TRI[END_REF]. Some tools deal with the problem of non existing preference model solutions which may arise because of an inconsistent set of assignment examples (i.e. assignment examples that do not match ELECTRE TRI) [START_REF] Mousseau | Resolving inconsistencies among constraints on the parameters of an MCDA model[END_REF][START_REF] Mousseau | Dealing with inconsistent judgments in multiple criteria sorting models[END_REF]). While the above approaches target a unique DM, [START_REF] Damart | Supporting groups in sorting decisions: Methodology and use of a multi-criteria aggregation/disaggregation DSS[END_REF] propose a method involving a group of DMs that iteratively build, in parallel, in- dividual preference models and a collective preference model representing the group consensus. Table 1 presents a summary of the available tools for indirect preference elicitation related to ELECTRE TRI.

BUILDING A RISK EVALUATION MODEL WITH MULTIPLE EXPERTS

Contrary to the above approaches, the method we propose consists in eliciting category limits shared by all DMs while leaving weights of criteria possibly different for each DM. This enables a consensus to be approached by first validating a part of the preference model (the category limits) with all DMs before eliciting the weights of criteria. The process is the following.

• Obtain from each DM typical zones that correspond to each specific risk levels. These zones are defined by their evaluations on the criteria and will be used as indirect information to elicit the risk evaluation model.

• Search for an ELECTRE TRI model without vetoes representing the zone examples. A mathematical program is thus built in order to compute category limits compatible with the typical zones provided by the DMs. The program is defined in Section 5.

• If no ELECTRE TRI model without vetoes is able to represent all zone examples, it is possible to allow vetoes to be used and define the corresponding mathematical program.

• If still no satisfying ELECTRE TRI model is found, it is possible to check which maximal subsets of zone examples can be represented [START_REF] Mousseau | Resolving inconsistencies among constraints on the parameters of an MCDA model[END_REF][START_REF] Mousseau | Dealing with inconsistent judgments in multiple criteria sorting models[END_REF]).

The DMs may then be asked individually if they agree to remove or change specific zone examples to restore consistency.

• When shared category limits are found, they should be presented to the DMs for validation. If they disagree, they may propose additional zone examples to further constrain the model and iteratively converge towards satisfactory category limits.

• At this stage, the DMs agree on a set of category limits but weights of criteria are still possibly distinct for each DM. The approach suggested by [START_REF] Damart | Supporting groups in sorting decisions: Methodology and use of a multi-criteria aggregation/disaggregation DSS[END_REF] may then be used to build a consensus on the weights.

One of the important features of the proposed approach is that, as opposed to existing methods, the proposed method applies to a group DM context and does not suppose that part of the preference model is known beforehand.

ILLUSTRATIVE EXAMPLE

Let us illustrate the method on the following hypothetical scenario. This simplistic example role is to illustrate the proposed method rather than provide an in-depth case study.

A group of experts would like to develop a scale permitting to evaluate the level of risk of each territorial zone around a given industrial installation related to a possible hazard. Each zone is to be determined as belonging to one of the three categories {High risk ≺ Medium risk ≺ Low risk}. Each of these categories are associated with specific precautionary measures (e.g. evacuate the population). Building a risk evaluation scale relative to hazards related to industrial installations is typically important in decision processes such as the PPRT (Plan de Prévention des Risques Technologiques) used in France (PPRT 2011) .

The four members of the expert group consider that the following six criteria should be used to evaluate the risk associated to each zone. d The distance to the hazardous industrial installation, evaluated on a 3 points ordinal scale. 0: less than 500 meters; 1: between 500 and 2000 meters; 2: more than 2 km.

sc Presence or absence school in the zone, a binary assessment. No school is encoded with a 1, presence of school with a 0.

bv Building vulnerability, evaluated on a 5 points ordinal scale.

pe Presence of public environment and technical assets, evaluated on a 5 points ordinal scale.

i Impact on other industrial installations, which could lead to cascaded effects. Evaluated on a 5 points ordinal scale. We suppose that each expert provides 30 zone examples and their associated risk evaluation (that could correspond to real zones the experts have previously evaluated or fictitious zones defined by their evaluation vectors). The data is given in Annex A.

Using the assignment examples, the program Infer Limits (described in Section 5) is used to find category limits shared by the experts which match their individual assignment examples. The results are provided in Table 2.

The numerical values found as frontiers by the mathematical program should be interpreted in terms of the scales used for the evaluations of the related point of views. For example, the evaluation of the frontier b 2 on the point of view of the building vulnerability criterion, as displayed in Table 2, is 3.2. As that criterion scale is discrete, a value of at least 3.2 corresponds to a value of at least 4. This means that an evaluation of 4 on the point of view of the building vulnerability for a given zone counts as an argument in favor of that zone to be assigned to a category better than b 2 , i.e. the category "Low risk". If the zone has an evaluation of at least 2.8, the evaluation of b 1 on the building vulnerability criterion, but not better than 3.2, hence, an evaluation of 3, this point of view argues in favor of that zone being assigned to the risk level "Medium risk". If the evaluation is less than 2.8, thus is 0 to 2, the risk level recommanded by this point of view is "High risk".

The table shows a possible set of frontiers evaluations that may be shared by all DMs such that when used in an ELECTRE TRI model it is possible to reproduce their zone examples, provided adequate weights are used. At this stage the DMs do not share the weight values yet. Table 3 shows for each DM a set of weights matching the assignment examples with the common category limits. Note that these weights are not the only ones that reproduce all assignment examples with the frontier values shown in Table 2.

Choosing these category limits impose some constraints, for each expert, on the criteria importance that may be used to match the zone examples. It is possible to compute the restrictions imposed to each expert on the relative importance of the criteria. In Figure 1, for a given expert, an arrow is drawn from a criterion j to a criterion j ′ if all ELECTRE TRI models matching the zone examples for that expert, considering the category limits in Table 2, are such that the Figure 1: The restrictions on the weights imposed by the choice of the common frontiers, for each expert. An arrow from a criterion j to a criterion j ′ states that the choice of these frontiers implies a greater weight for j than for j ′ for this expert. weight on the criterion j is more important than the weight on the criterion j ′ . These importance relations are to be discussed and accepted by all experts for these shared frontiers to be accepted. If a DM disagree on some importance comparisons, the Infer Limits program should be run again with some additional constraints concerning unacceptable weights comparisons.

It may be seen that a consensus on the weights cannot be readily found using these category limits as e.g. expert 3 has sc necessarily more important than bv while experts 1 and 4 have the inverse preference. To get closer to a consensus, the third decision expert might be asked if he could consider bv to be more important than sc and the Infer Limits program be run again with that supplementary constraint included.

In our example, we suppose that when presented the category limits, one of the expert disagree with one of the values, saying that in his opinion the category limit b 1 on the criterion pe must have a value of at least 2, thus, that a zone should have an evaluation of at least 2 according to the point of view of the public environment and technical assets for this point of view to argue in favor of the zone to be in the "Medium risk" or "Low risk" category. This constraint is thus included in the program and a new so- 4.

A set of weights compatible with the zone examples using these frontiers is presented in Table 5.

Supposing that all DMs agree with the shared frontier values found, the approach suggested by [START_REF] Damart | Supporting groups in sorting decisions: Methodology and use of a multi-criteria aggregation/disaggregation DSS[END_REF] may then be used to iteratively build consensual weights among the group of experts.

ALGORITHMIC DESCRIPTION

In this section we detail the mathematical program which allow to elicit the category limits for a set of DMs. This requires the following supplementary notation. For each DM l ∈ L, the set of examples E l is the set of pairs (a, h) ∈ A × [1, k] specifying that the zone a is assigned to the risk level c h by l. We write a -→ l h to denote such a zone example.

Infer Limits is a mixed integer program that finds, if possible, a set of category limits such that it is possible to satisfy the zone examples of each DMs using individual criteria importance parameters without using veto thresholds.

Stating the problem

Having a set of zones A used as examples, a set of criteria indices J, the evaluations of the zones g j (a)∀a ∈ A, j ∈ J, the number of categories k, a set of DMs L, zone examples E l , the goal of Infer Limits is to determine the performances of category limits g j (b h ), ∀j ∈ J, 1 ≤ h ≤ k -1 shared among the DMs, together with individual weights w l j and majority thresholds λ l , matching all zone examples. Note that the lowest category limit, b 0 , must be outranked by every zone (because they must go into at least category c 1 , the worst one) and hence the performances g j (b 0 ) on every criterion necessarily consist in performances lower than every other performances on the same criterion. The upper frontier of the best category, that would be b k , is used only when an other assignment procedure known as the optimistic procedure is used. Therefore, the lowest and best frontiers are not considered in our inference programs. Satisfying these assignment examples amounts to satisfy, ∀l ∈ L, ∀a -→ l h ∈ E l : j∈J:g j (a)≥g j (b h-1 )

w l j ≥ λ l , and (2)
j∈J:g j (a)≥g j (b h ) w l j < λ l .

(3) Equation (2) ensures that the example zone is assigned to a risk level at least as good as c h , and (3) makes sure that it is assigned to a risk level not better than c h .

Constraints

∀j ∈ J, a ∈ A, 1 ≤ h ≤ k -1 :            1 M j + δ j ((g j (a) -g j (b h )) + δ j ) ≤ C j (a, b h ); C j (a, b h ) ≤ 1 M j (g j (a) -g j (b h )) + 1. (4)
The constraints in (4) define C j (a, b h ), binary variables indicating the agreement of the criterion j to say that a is at least as good as b h . See Figure 2. δ j , ∀j ∈ J, is a value smaller than the minimal difference between any two different performances on j (∀j ∈ J : min j = min a,b∈A |g j (a)g j (b)|) divided by k. If there is no two different performances for a given criterion, δ j is defined arbitrarily to 1. It is defined so that adding this value to the constraints does not restrict the space of possibilities about g j (b h ) and may thus be used to change the strict inequality to a loose inequality. ∀j ∈ J, M j = g jg j + min j is a scaling factor representing the maximal performance difference on the criterion j including one extreme frontier. g j and g j are the best and worst performances on the criterion j.

The sum of the weights are modelled using the following constraints, ∀l ∈ L, a -→

l h ∈ E l , j ∈ J :                      σ l j (a, b h ) ≤ w l j σ l j (a, b h ) ≥ 0 σ l j (a, b h ) ≤ C j (a, b h ) σ l j (a, b h ) ≥ C j (a, b h ) + w l j -1. (5)
The variables σ l j (a, b h ) represent the sum of the support for saying that a is at least as good as b h

-M j 0 M j g j (a) -g j (b h ) C j C j = 1 M j (g j (a) -g j (b h )) + 1 C j = 1 M j (g j (a) -g j (b h )) | ■ × Figure 3: constraining σ l j to the appropriate value 0 w j 1 0 1 C j σ j σ l j = C j σ l j = C j + w j -1 ⑦ ⑥ ✻ ❄ ❡
while avoiding the non-linear expression σ l j (a, b h ) = w l j C j (a, b h ) [START_REF] Meyer | September). Disaggregation of bipolar-valued outranking relations[END_REF]). See also Figure 3.

Then we need to ensure that the sum of supports for each examples are greater than the related majority thresholds, i.e., ∀l ∈ L, ∀a -→

l h ∈ E l , h < k : j∈J σ l j (a, b h ) + s ≤ λ l -δ λ , (6)
and ∀l ∈ L, ∀a -→

l h ∈ E l , h > 1 : j∈J σ l j (a, b h-1 ) ≥ λ l + s. ( 7
)
δ λ is an arbitrary small positive value, s is a slack variable used in the objective function.

∀l ∈ L :

j∈J w l j = 1. (8)
Constraints ( 8) ensure that the weights sum to one.

∀j ∈ J, 2 ≤ h ≤ k -1 : g j (b h-1 ) ≤ g j (b h ). (9)
Ensuring that the category limits are ordered correctly is achieved through Constraints (9).

The last constraints ensures that the best and worst frontiers evaluations are not more than min j apart from the closest zone evaluation. ∀j ∈ J :    g j (b k -1) ≤ g j + min j , g j (b 1 ) ≥ g jmin j .

(10)

Objective function

In order to maximize the separation between the sum of support and the majority thresholds, the slack variable s used in Constraints ( 6), ( 7) must be maximized.

It is also possible to simply search for a solution to the program without specifying any objective.

CONCLUSION

In this article we show how to express the elaboration of an ordinal multicriteria risk scale as a multicriteria sorting problem. This highlights the usefulness for risk analysis of some tools and approaches developed in the field of multicriteria decision analysis. Moreover, we propose a method to elicit an appropriate multicriteria sorting model with multiple experts on the basis of the ordinal evaluation of typical zones.

The method emphasizes and fosters the necessary interaction between experts during the elicitation of the model.

dm3 Zone21 5 5 3 1 1 50.6 Medium dm3 Zone31 2 3 3 3 1 0.5 Medium dm3 Zone46 1 2 3 1 1 22.2 High dm3 Zone50 3 5 3 2 0 90.6 High dm3 Zone51 4 4 2 1 0 32.1 High dm3 Zone55 2 5 1 4 1 16.7 High dm3 Zone58 2 1 2 3 1 74.2 Medium dm3 Zone59 2 1 2 2 0 82.0 High dm3 Zone64 2 2 3 2 1 40.7 Medium dm3 Zone69 3 3 1 2 1 56.9 Medium dm3 Zone70 3 1 1 2 0 51.4 High dm3 Zone71 4 3 1 5 0 83.0 Low dm3 Zone72 4 3 3 4 0 99.7 Low dm3 Zone74 1 5 2 4 0 17.0 High dm3 Zone75 1 5 2 4 0 85.7 Medium dm3 Zone78 1 2 2 5 1 84.1 Low dm3 Zone81 3 2 1 5 0

Figure 2 :

 2 Figure 2: constraining C j to the appropriate value

Table 1 :

 1 A summary of the main features proposed by other articles and by this one (last row). For each article, the second column indicates the expected input of the main tool proposed in the article, the last one shows its output. i designates assignment examples from a single DM, i * designates possibly inconsistent assignment examples from a single DM, g is a group of DMs' assignment examples, P is a set of category limits, W is a set of weights. The computations are based on linear (or mixed integer and linear) programming, except for the first one.

	Article	input output
	MS98	i	P, W (non-linear)
	MFN01	i, P	W
	NM00	i, W P
	DMFC02	i	robust model (P, W)
	MDFGC03 i *	how to restore consistency
	MDF06	i *	how to restore consistency
	DDM07	g, P	progressive collective model (W)
	This article g	collective model (P)

Table 2 :

 2 Inferred category limits. b 1 separate the categories "High risk" and "Medium risk" and b 2 separate the categories "Medium risk" and "Low risk".

	Frontier bv pe	d	i	sc	pr
	b 1	2.8 1.2 1.8 2.8 0.8 40.7
	b 2	3.2 3.2 2.2 3.2 0.8 81.8
	pr Proportion of the population who is vulnerable
	(children under 15 and elderly persons), evalu-
	ated in percentage.

Table 3 :

 3 A set of weights, as found by the Infer Limits program, matching zone examples of each expert when used together with the inferred category limits.

Table 4 :

 4 Inferred category limits on the second run, with the added constraint that b 1 on pe must be at least 2. b 1 separate the categories "High risk" and "Medium risk" and b 2 separate the categories "Medium risk" and "Low risk".

	Frontier bv pe	d	i	sc	pr
	b 1	2.8 2.8 1.8 2.2 0.8 40.8
	b 2	3.2 3.2 2.2 3.8 0.8 74.2

Table 5 :

 5 A set of weights, as found by the Infer Limits program, matching zone examples of each expert when used together with the inferred category limits.