François Septier
email: francois.septier@univ-ubs.fr

Patrick Armand

Christophe Duchenne

Sequential Monte Carlo Sampler Applied to Source Term Estimation in Complex Atmospheric Environments

Keywords: Source term estimation, Bayesian approach, dispersion model, adaptivity, sequential importance sampling, synthetic example

The accurate and rapid reconstruction of a pollution source represents an important but challenging problem. Several strategies have been proposed to tackle this issue among which we find the Bayesian solutions that have the interesting ability to provide a complete characterization of the source parameters through their posterior probability density function. However, these existing techniques have certain limitations such as their computational complexity, the required model assumptions, their difficulty to converge, the sensitive choice of model/algorithm parameters which clearly limit their easy use in practical scenarios. In this paper, to overcome these limitations, we propose a novel Bayesian solution based on a general and flexible population-based Monte Carlo algorithm, namely the sequential Monte Carlo sampler. Owing to its full adaptivity through the learning process, the main advantage of such an algorithm lies in its capability to be used without requiring any specific assumptions on the underlying statistical model and also without requiring from the user any difficult choices of certain parameter values. The performance of the proposed inference strategy is assessed using twin experiments in complex built-up environments.

Introduction

The threat of chemical, radiological, biological, and nuclear (CRBN) releases raises some complex and challenging scientific issues. In the event of a CBRN incident, it is of great importance to have as soon as possible an accurate assessment of the damage likely to be caused by the release which is usually undertaken 5 using an atmospheric dispersion model of the contaminant. However, such spatial and temporal forecast, especially in complex built-up environments, is only possible by providing the source term parameters to the dispersion model. As a consequence, the main objective consists in rapidly obtaining an accurate estimation of the source parameters from noisy measurements of the concentration 10 levels observed by several sensors deployed in a specific surveillance area.

Due to its ill-posed nature, this challenging inverse problem has been tackled using the proposition of different algorithms which can be grouped into two categories. The first one aims at obtaining a single point estimate of the unknown parameters by solving an optimization problem where a cost function has to be minimized using least squares or genetic algorithms, e.g [START_REF] Winiarek | Estimation of errors in the inverse modeling of accidental release of atmospheric pollutant: Application to the reconstruction of the cesium-137 and iodine-131 source terms from the Fukushima Daiichi power plant[END_REF].

Unfortunately, such approaches do not allow us to quantify the uncertainty associated to the unknowns, which could be really problematic in such an important context. To overcome this limitation, Bayesian algorithms have been designed to solve an inference problem by providing the complete probability density function of the parameters of interest given the measurements observed by a network of sensors. Owing to the complex and nonlinear nature of the source term estimation (STE) model, the exact computation of such a distribution is not feasible analytically, and one has to resort to some approximation techniques such as stochastic simulation algorithms.

More specifically, these techniques consists in obtaining samples from the distribution of interest by using either some Markov Chain Monte Carlo (MCMC) kernel [START_REF] Delle Monache | Bayesian Inference and Markov Chain Monte Carlo Sampling to Reconstruct a Contaminant Source on a Continental Scale[END_REF][START_REF] Chow | Source inversion for contaminant plume dispersion in urban environments using building-resolving simulations[END_REF]Keats et al., 2007a;[START_REF] Yee | Bayesian Inference for Source Reconstruction: A Real-World Application[END_REF] or using the principle of importance sampling [START_REF] Septier | Report No. 2 -Source Term estimation[END_REF][START_REF] Septier | A Bayesian inference procedure based on inverse dispersion modelling for source term estimation in built-up environments[END_REF]. As recently discussed in [START_REF] Septier | A Bayesian inference procedure based on inverse dispersion modelling for source term estimation in built-up environments[END_REF], the application of an Adaptive Multiple Importance Sampling (AMIS) technique on the challenging STE problem allows us to obtain significant gain compared to state-of-the-art MCMC algorithms in both synthetic and real data experiments.

Unfortunately, this AMIS algorithm has been designed for a specific observation model and in particular, it is assumed that both the measurement errors and the emission rate levels follow a normal distribution which is not fully appropriate to real situations since negative values can appear in the concentration levels and have to be removed to keep the physical meaning of the data. However, relaxing these model assumptions will result in an increasing dimension of the state to estimate, thus leading generally to poorer estimation results as the performance of this AMIS algorithm is quite sensitive to the dimension related to the unknown parameters. Moreover, the parameters of these normal distributions have to be set a priori by the user which could clearly be difficult in practice.

In this work, we propose a novel Bayesian STE strategy based on a more general and flexible population-based Monte Carlo framework. Originally proposed in [START_REF] Moral | Sequential Monte Carlo samplers[END_REF][START_REF] Peters | Topics in Sequential Monte Carlo Samplers[END_REF], this technique called Sequential Monte-Carlo (SMC) sampler is a class of sampling algorithms which combine importance sampling and resampling. They have been primarily used as "particle filter" to solve optimal filtering problems [START_REF] Doucet | On sequential Monte-Carlo sampling methods for Bayesian filtering[END_REF]. In this context, SMC methods/particle filters have benefited from wide-spread use in various applications (tracking, computer vision, digital communications) due to the fact that they provide a simple way of approximating complex filtering distribution sequentially in time. But in (Del [START_REF] Moral | Sequential Monte Carlo samplers[END_REF], the authors developed a general framework that allows SMC to be used to simulate from a single and static target distribution, thus becoming a promising alternative to standard MCMC methods. The SMC sampler framework involves the construction of a sequence of artificial distributions on spaces of increasing dimensions which admit the distributions of interests as particular marginals. The mechanism is similar to sequential importance sampling (resampling) [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF][START_REF] Liu | Monte Carlo Strategies in Scientific Computing[END_REF], with one of the crucial differences being the framework under which the random samples, also called particles, are allowed to move, resulting in differences in the calculation of the weights of the particles. These methods have several advantages over traditional and populationbased MCMC methods. Firstly, unlike MCMC, SMC methods do not require any burn-in period and do not face the sometimes contentious issue of diagnosing convergence of a Markov chain. Secondly, as discussed in [START_REF] Jasra | On population-based simulation for static inference[END_REF], compared to population-based MCMC, SMC sampler is a richer method since there is substantially more freedom in specifying the mutation kernels in SMC: kernels do not need to be reversible or even Markov. As a consequence, adaptive proposal distributions can be easily used, thus giving a lot more of opportunities to improve its efficiency. Moreover, unlike MCMC, SMC samplers provide an unbiased estimate of the normalizing constant of the posterior distribution which can be one quantity of interest in the inference problem to deal with. Indeed, this normalizing constant is the marginal likelihood and therefore could be used for selecting some model assumptions (e.g. choice of distribution for the measurement noise or for the prior of parameters). Let us finally denote that many other inference technique based on importance sampling such as Annealed Importance sampling [START_REF] Neal | Annealed Importance Sampling[END_REF], population Monte Carlo [START_REF] Cappé | Population monte carlo[END_REF] and its more advanced variants like the Adaptive Multiple Importance Sampling (AMIS) [START_REF] Cornuet | Adaptive Multiple Importance Sampling[END_REF] can all be considered as special cases of the SMC sampler.

Although this approach presents many advantages over traditional MCMC methods, the potential of these emergent techniques is however largely underexploited in practice. In [START_REF] Nguyen | Efficient Sequential Monte-Carlo Samplers for Bayesian Inference[END_REF], some strategies have been described and proposed in order to make easier its efficient implementation in practical problems. In this work, we propose to design an STE algorithm based on this sequential Monte Carlo sampler. The rest of this paper is organized as follows. Section 2 describes the statistical model and describes the Bayesian framework used for this STE problem. In Section 3, the general principle of the SMC sampler is firstly presented, then we describe the proposed SMC sampler applied to the STE problem. Numerical simulations are performed in Section 4 to assess the performances of the proposed approach. Conclusions are finally given in Section 5.

Statistical Model of the STE Problem

This section firstly describes the statistical model of the source term estimation problem, and then the Bayesian framework for estimating the characteristics of the source is discussed.

Atmospheric dispersion model

In this paper, a point-wise and static source fully characterized by the parameter β = [x s , q] is considered where x s = [x s , y s ] stands for the spatial position of the source and q is the release rate vector resulting from the discretization of the plausible emission time window into T s time intervals. A network of N c sensors is deployed over a 2-dimensional monitoring region to measure the concentration levels. By using an atmospheric dispersion model, the output simulated concentration at the location of the i-th sensor at time t j is defined as

y i,j = Ts n=1 q n C i,j (x s , ∆t n ), (1)
where j = 1, . . . , T c with T c the number of time samples collected by each sensor. Each value results from the superposition of the T s releases at the different time steps {∆t n } Ts n=1 weighted by their associated emission rates {q n } Ts n=1 of the source. C i,j (x s , ∆t n ) corresponds therefore to the simulated concentration obtained at the i-th sensor at time t j if a unitary release is made during the time step ∆t n from a source that is located at x s . Let us note that the proposed method can be used for any specific choice of atmospheric dispersion model as long as we are able to obtain C i,j (x s , ∆t n ). All these simulated concentration values obtained at the different time samples of all sensors can be written in the following matrix form:

y = C(x s )q, ( 2
)
where

y = y 1,1 • • • y 1,Tc • • • y Nc,1 • • • y Nc,Tc
T is the vector of simulated concentration values from the used atmospheric dispersion model and C(x s ), generally called source-receptor matrix [START_REF] Seibert | Source-receptor matrix calculation with a Lagrangian particle dispersion model in backward mode[END_REF], takes the following matrix form

C(x s ) =              C 1,1 (x s , ∆t 1 ) • • • C 1,1 (x s , ∆t Ts ) . . . . . . . . . C 1,Tc (x s , ∆t 1 ) • • • C 1,Tc (x s , ∆t Ts ) . . . . . . . . . C Nc,1 (x s , ∆t 1 ) • • • C Nc,1 (x s , ∆t Ts ) . . . . . . . . . C Nc,Tc (x s , ∆t 1 ) • • • C Nc,Tc (x s , ∆t Ts )              . (3)
The computation of this source-receptor matrix is an important part in an STE procedure as it links the source characteristics with the measurements and quantifies the predicted concentration value at some location and time from a dispersion model for a given source. As a consequence, in stochastic simulation based inference techniques, this matrix has to be computed for each generated sample which are quite numerous (at least several thousands) for a satisfactory estimation accuracy level. The computation of this matrix with a Lagrangian particle dispersion model (LPDM) in a forward mode constitutes the most timeconsuming step of the algorithm proposed in [START_REF] Rajaona | An adaptive Bayesian inference algorithm to estimate the parameters of a hazardous atmospheric release[END_REF]. In this study as proposed in (Keats et al., 2007a;[START_REF] Yee | Bayesian inversion of concentration data: Source reconstruction in the adjoint representation of atmospheric diffusion[END_REF][START_REF] Septier | A Bayesian inference procedure based on inverse dispersion modelling for source term estimation in built-up environments[END_REF], we use an alternative strategy which consists in using instead the backward mode of a LPDM. Using this backward mode is computationally advantageous if the number of receptors is less than the number of sources, which is generally the case in practice.

Definition of the observation model

The likelihood function expresses the probability of observing a specific set of concentrations from the N c sensors sensors given all the parameters of the source β. In other words, this likelihood function provides a probabilistic information about the discrepancy between the measured concentration values, denoted by the vector y, and the simulated concentration values, y, obtained from the dispersion model and defined in Eq. ( 2). The chosen parametric distribution for this likelihood should therefore characterize the three classical types of error: the dispersion modeling error, the observation error and the representativeness error due to the interpolation in both time and space of the dispersion model [START_REF] Koohkan | Accounting for representativeness errors in the inversion of atmospheric constituent emissions: application to the retrieval of regional carbon monoxide fluxes[END_REF]. Since the y is a complex but deterministic mapping of the source parameters β, let us denote the likelihood function by p(y| y;

Ψ l ) = p(y|β; Ψ l ), (4)
where Ψ l represents all the parameters of the specific parametric likelihood chosen in the study. The parameters Ψ l are assumed either known by a calibration from historical data or unknown and therefore need to be jointly estimated with the source parameters β. For the rest of the paper, we will keep this general formulation of the likelihood since the estimation method we propose in this work does not depend on the choice of this likelihood -the only assumption we need is the ability of evaluating this function point-wise.

As an illustration of a possible plausible choice of likelihood function which has been used in [START_REF] Lewellen | Analysis of concentration fluctuations from lidar observations of atmospheric plumes[END_REF][START_REF] Septier | Tracking of Multiple Contaminant Clouds[END_REF] in order to ensure the positivity of measurements, the noisy concentration of pollutant at a given location and time can be considered as conditionally independent and can be modeled by a clipped normal distribution, i.e.

p(y| y;

Ψ l ) = Nc i=1 Tc n=1 CN y i,n ; y i,n , σ 2 ǫ , (5)
where

CN y; µ, σ 2 ǫ =    0 if y < 0 Φ(-µ σǫ ) if y = 0 CN y i,n ; y i,n , σ 2 ǫ otherwise (6)
where Φ(y) is the cumulative distribution function (CDF) of the standard normal distribution. The delta function at zero corresponds to the intermittency (periods of zero concentration) in observed concentrations. With such a choice, the (possibly unknown) parameter of the distribution is Ψ l = σ 2 ǫ which corresponds to the variance (level) of the noise in the measurements. This likelihood is clearly more adapted to the problem under study than the traditional normal distribution used for example in [START_REF] Yee | Bayesian inversion of concentration data: Source reconstruction in the adjoint representation of atmospheric diffusion[END_REF][START_REF] Rajaona | An adaptive Bayesian inference algorithm to estimate the parameters of a hazardous atmospheric release[END_REF][START_REF] Septier | A Bayesian inference procedure based on inverse dispersion modelling for source term estimation in built-up environments[END_REF].

A priori knowledge about the source parameters

Our belief regarding the characteristics of the unknown state of interest, β, is encapsulated within the prior probability distributions of the proposed Bayesian model. As pointed out above with the likelihood, the method we propose in this work could be applied for any choice of prior probability distribution for the source parameters β as soon as sampling from such distribution and its point-wise evaluation are possible.

As an illustration, if we just know that the release could appear anywhere uniformly in the region of surveillance denoted here by Ω ⊆ R 2 , the following uniform prior distribution could be chosen for the position of the source:

p(x s ) = U Ω (x s ). ( 7
)
Of course, in some scenarios of interest, it could be more appropriate to incorporate a more informative distribution to represent our initial guess about this source location (nuclear plants, industrial sites, etc).

Regarding now the emission rate vector, q, as in (Keats et al., 2007b), we consider a point source which releases material at a steady rate of α and whose turn-on and turn-off times are t on and t off , i.e.:

q n = α1 [ton;t off ] (n), (8)
with 1(•) denoting the indicator function. As a consequence, the emission rate vector q is fully characterized by the three parameters: α ∈ R + , (t on , t off ) ∈ {1, . . . , T s } with t on ≤ t off . However, in practice, the emission rate α is assumed to be bounded such as

Q min ≤ α ≤ Q max .
If nothing else is known a priori about these unknown parameters, the following uniform distribution for these parameters can be chosen as prior distribution:

p(α) = U [Qmin;Qmax] (α), p(t on , t off ) ∝ 1 {1,...,Ts} (t on ) × 1 {ton,...,Ts} (t off ). ( 9
)
By using such simple prior distributions, the unknown parameters that fully characterize the source is given by β = [x s , α, t on , t off ] and moreover there does not exist any unknown hyperparameters, so Ψ p = ∅.

Source term estimation in a Bayesian framework

In this work, a Bayesian solution is considered in order to solve efficiently the challenging STE problem. Instead of just a point-wise estimation of the source characteristics, β, we are therefore interested in obtaining the full posterior distribution of the unknown parameters, p(β|y), which completely characterizes the available information on β given the measurements y obtained from all the sensors deployed in the field. With such a quantity, one can obtain all possible quantities of interest about the parameters such as, for example, point estimates and more importantly confidence intervals. Additionally in this problem, we want to jointly estimate the unknown parameters of both the prior and the likelihood probability functions, denoted respectively by Ψ p and Ψ l . As a consequence, the complete Bayesian solution is to consider the posterior distribution on the extended space (β, Ψ l , Ψ p ) which can be decomposed as follows

p(β, Ψ l , Ψ p |y) = p(y|β, Ψ l )p(β|Ψ p )p(Ψ p , Ψ l ) p(y) , (10)
where p(Ψ p , Ψ l ) represents the prior distribution on the unknown parameters (Ψ p , Ψ l ).

Unfortunately, even with some particular choice for the likelihood and prior distributions, this joint posterior distribution of interest in ( 10) is analytically intractable. Indeed, the dependence of the position of the source in the measurements is highly nonlinear due to the complex structure of the source-receptor matrix C(x s ) as discussed in Section 2.1. In this work, we will develop an efficient stochastic simulation based algorithm to approximate this complex posterior distribution p(β, Ψ l , Ψ p |y).

Proposed SMC Sampler to STE

In this section, we first present the general principle of the SMC sampler then we describe how to apply the SMC sampler to efficiently solve the STE problem.

General Idea of the SMC Sampler

The SMC Sampler methodology is a generic approach to approximate a sequence of probability distributions {π t } T t=1 regarding some parameters of interest θ and which are defined upon a common measurable space E (Del [START_REF] Moral | Sequential Monte Carlo samplers[END_REF], where the final distribution π T is the distribution of interest. The main ideas of the SMC sampler can be summarized as: a) Instead of sampling directly from the complex target posterior distribution, a sequence of intermediate {π t } T t=1 is designed such that the transition from a simpler distribution, π 0 to the one of interest, π T , is smooth. b) The problem is then solved by transforming this problem in the standard SMC filtering framework, where the sequence of target distributions on the path-space denoted by {π t } T t=1 , which admits {π t } T t=1 as marginals, is defined on the product space, i.e. supp(π

t ) = E × • • • × E = E t . This novel sequence of target distributions πt is defined as follows: πt (θ 1:t ) = γt (θ 1:t ) Z t (11)
where γt (θ

1:t ) = γ t (θ t ) t-1 k=1 L k (θ k+1 , θ k ) (12)
in which the artificial kernels {L k } t-1 k=1 are called backward Markov kernels since L t (θ t+1 , θ t ) denotes the probability density of moving back from θ t+1 to θ t . γ t and γt represent the unnormalized version of the target distribution π t and πt , respectively. Z t thus corresponds to the normalizing constant of the target distribution πt . By using such a sequence of extended target distributions {π t } T t=1 based on the introduction of backward kernels {L k } t-1 k=1 , sequential importance sampling can be used in the same manner as standard SMC filtering algorithms.

Within this framework, one may then work with the constructed sequence of distributions, π t , under the standard SMC algorithm [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF]. In summary, the SMC sampler algorithm involves three stages:

1. mutation, where the particles are moved from θ t-1 to θ t via a mutation kernel K t (θ t-1 , θ t ); 2. correction, where the particles are reweighted with respect to π t via the incremental importance weight (Equation ( 16)); and 3. selection, where according to some measure of particle diversity, such as effective sample size, the weighted particles may be resampled in order to reduce the variability of the importance weights.

In more detail, suppose that at time t-1, we have a set of weighted particles

θ (m) 1:t-1 , W (m) t-1 N m=1 that approximates πt-1 via the empirical measure πN t-1 (dθ 1:t-1 ) = N m=1 W (m) t-1 δ θ (m) 1:t-1 (dθ 1:t-1 ) (13)
where δ(•) is the Dirac delta function. These particles are first propagated to the next distribution πt using a Markov kernel K t (θ t-1 , θ t ) to obtain the set of particles θ

(m) 1:t N m=1
. IS is then used to correct for the discrepancy between the sampling distribution η t (θ 1:t ) defined as

η t (θ (m) 1:t ) = η 1 (θ (m) 1 ) t k=2 K k (θ (m) t-1 , θ (m) t ) ( 14
)
and πt (θ 1:t ). In this case the new expression for the unnormalized importance weights is given by

W (m) t ∝ πt (θ (m) 1:t ) η t (θ (m) 1:t ) = π t (θ (m) t ) t-1 s=1 L s (θ (m) s+1 , θ (m) s ) η 1 (θ (m) 1 ) t k=2 K k (θ (m) t-1 , θ (m) t ) ∝ w t (θ (m) t-1 , θ (m) t )W (m) t-1
(15) where w t , termed the (unnormalized) incremental weights, are calculated as,

190 w t (θ (m) t-1 , θ (m) t ) = γ t (θ (m) t )L t-1 (θ (m) t , θ (m) t-1 ) γ t-1 (θ (m) t-1 )K t (θ (m) t-1 , θ (m) t ) (16)
However, as in the particle filter, since the discrepancy between the target distribution πt and the proposal η t increases with t, the variance of the unnormalized importance weights tends therefore to increase as well, leading to a degeneracy of the particle approximation. A common criterion used in practice to check this problem is the effective sample size ESS, which is given by:

ESS t = N m=1 ( W (m) t ) 2 -1 = N m=1 W (m) t-1 w t (θ (m) t-1 , θ (m) t ) 2 N j=1 W (j) t-1 2 w t (θ (j) t-1 , θ (j) t ) 2 (17)
If the degeneracy is too high, i.e., the ESS t is below a prespecified threshold, ESS, then a resampling step is performed. The particles with low weights are discarded whereas particles with high weights are duplicated. After resampling, the particles are equally weighted.

The final weighted particles at distribution π T are considered weighted samples from the target distribution π of interest. As a consequence, the SMC sampler provides an estimate of this distribution

π N T (dθ) = N m=1 W (m) T δ θ (m) T (dθ) (18)
and the estimation of expectation of some function ϕ(•) with respect to the target distribution of interest is given by

E π N T [ϕ(θ)] = N m=1 W (m) T ϕ(θ (m) T ) (19)

Algorithm settings 195

As described in Section 2.4 -Eq. ( 10), the target distribution of interest in our case is defined as the posterior distribution over the unknown parameters, i.e.

π(θ) = p(θ|y), (20)
where all the parameters of interest that have to be jointly estimated is denoted by θ = {β, Ψ l , Ψ p }. The algorithm presented in the previous subsection is very general. There is a wide range of possible choices to consider when designing an SMC sampler algorithm, the appropriate sequence of distributions {π t } 1≤t≤T , the choice of both the mutation kernel {K t } 2≤t≤T and the backward mutation

200 kernel {L t-1 } T t=2
(for a given mutation kernels), see details in [START_REF] Moral | Sequential Monte Carlo samplers[END_REF]. In this subsection, we provide a discussion on how to choose these parameters of the algorithm for our estimation problem.

Sequence of distributions π t

There are many potential choices for {π t } leading to various integration and optimization algorithms. As a special case, we can set π t = π for all t ∈ N . Alternatively, to maximize π(θ), we could consider π t (θ t ) = [π(θ t )] ξt for an increasing schedule {ξ t } t∈N to ensure π T (θ) is concentrated around the set of global maxima of π(θ). In the context of Bayesian inference for static parameters which is the main focus of this paper, we consider the likelihood tempered target sequence [START_REF] Neal | Annealed Importance Sampling[END_REF] π t (θ) ∝ p(θ)p(y|θ) φt (21

)
where {φ t } is a non-decreasing temperature schedule with φ 0 = 0 and φ T = 1. We thus sample initially from the prior distribution π 0 = p(θ) directly and introduce the effect of the likelihood gradually in order to obtain at the end t = T an approximation of the posterior distribution p(θ|y). Tempering the likelihood could significantly improve the exploration of the state space in complex multimodal posterior distribution. The approximation of an expectation with respect to the posterior distribution of interest is therefore given by:

E π N [ϕ(θ)] = N i=1 W (i) T ϕ(θ (i) T ) (22)
As discussed in [START_REF] Nguyen | Efficient Sequential Monte-Carlo Samplers for Bayesian Inference[END_REF], several statistical approaches have been proposed in order to automatically obtain such a tempering schedule via the optimization of some criteria, which are known as on-line schemes [START_REF] Jasra | Inference for lévy-driven stochastic volatility models via adaptive sequential monte carlo[END_REF][START_REF] Zhou | Toward Automatic Model Comparison: An Adaptive Sequential Monte Carlo Approach[END_REF]. In particular, [START_REF] Zhou | Toward Automatic Model Comparison: An Adaptive Sequential Monte Carlo Approach[END_REF] propose to set the temperature at the t-th iteration of the algorithm such that the conditional ESS (CESS), decays in a regular predefined way. This CESS given by

CESS t = N N i=1 W (i) t-1 w t (θ (i) t-1 , θ (i) t ) 2 N j=1 W (j) t-1 w t (θ (j) t-1 , θ (j) t ) 2 . ( 23
)
is a slight modification of the ESS defined in Eq. ( 17). As shown in [START_REF] Zhou | Toward Automatic Model Comparison: An Adaptive Sequential Monte Carlo Approach[END_REF], by characterizing how good an importance sampling proposal π t-1 would be for the estimation of expectation under π t , the use of this quantity rather than the ESS leads to a reduction in estimator variance.

Sequence of mutation kernels K t

The performance of SMC samplers depends heavily upon the selection of the transition kernels {K t } T t=2 and the auxiliary backward kernels {L t-1 } T t=2 . There are many possible choices for K t which have been discussed in [START_REF] Moral | Sequential Monte Carlo samplers[END_REF]. In this study, we propose to employ MCMC kernels of invariant distribution π t for K t . This is an attractive strategy since we can use the vast literature on the design of efficient MCMC algorithms to build a good importance distributions (See [START_REF] Robert | Monte Carlo statistical methods[END_REF]).

More precisely, since we are interested in complex models with potentially multimodal posterior distribution, a series Metropolis-within-Gibbs kernels with local moves [START_REF] Robert | Monte Carlo statistical methods[END_REF] will be employed in order to successively move:

• the position of the source x s ,

• the level of emission α,

• the time of emission (t on , t off ),

• the hyperparameters of the likelihood distribution Ψ l = σ 2 ǫ .

• and the hyperparameters of the prior distributions Ψ p , if any.

More specifically, after transforming the parameters σ 2 ǫ and α to take value in the real line by using the log(•) transform, an adaptive random walk proposal is used for each parameter except for the time of emission which values are discrete. This proposal consists, at the t-iteration of the SMC sampler, in adding a random perturbation, which follows generally a centered Normal distribution with covariance Σ rw •,t which is computed adaptively from the past simulations, to the current value of the parameter. As an example, for the position of the source, the proposed sample is obtained at the i-th move of the MCMC kernel of the t-th iteration of the SMC sampler as:

x * s,t = x i s,t + b xs,t (24)
in which b xs,t is a Gaussian random variable with zero mean and covariance matrix Σ rw xs,t . As with any sampling algorithm, faster mixing does not harm performance and in some cases will considerably improve it. In the particular case of Metropolis-Hastings kernels, the mixing speed relies on adequate proposal scales. As a consequence, we adopt the strategy proposed in [START_REF] Jasra | Inference for lévy-driven stochastic volatility models via adaptive sequential monte carlo[END_REF]. The authors applied an idea used within adaptive MCMC methods [START_REF] Andrieu | On the ergodicity properties of some adaptive mcmc algorithms[END_REF] to SMC samplers by using variance of parameters estimated from its particle system approximation as the proposal scale for the next iteration, i.e., the covariance matrix of the random-walk move for the position of the source is thus given by:

Σ rw xs,t = N m=1 W (m) t-1 x (m) s,t-1 -µ xs,t-1 x (m) s,t-1 -µ xs,t-1 T (25) with µ xs,t-1 = N m=1 W (m) t-1 x (m) s,t-1
The motivation is that if π t-1 is close to π t (which is recommended for an efficient SMC sampler algorithm), then the variance estimated at iteration t -1 will provide a sensible scaling at time t. In difficult problems such as the one addressed in this paper, other approaches could be added in order to have appropriate scaling adaptation; one approach demonstrated in [START_REF] Jasra | Inference for lévy-driven stochastic volatility models via adaptive sequential monte carlo[END_REF] is to simply employ a pair of acceptance rate thresholds and to alter the proposal scale from the simply estimated value whenever the acceptance rate falls outside those threshold values. This scheme is to ensure that the acceptance rates in the Metropolis-Hastings steps did not get too large or small. Through all this work, we use this procedure which consists for example to multiply the covariance matrix by 5 (resp. 1/5) if the rate exceeded 0.7 (resp. fell below 0.2). The same procedure is used respectively to sample respectively the log transform of the level of emission log(α) and of the likelihood variance log(σ 2 ǫ ). Finally, for the time of emission, the proposal distribution used is defined as: q(t on , t off |t i on,t , t i off,t ) = q(t on |t i on,t )q(t off |t on , t i off,t ),

where q(t on |t i on,t ) ∝ 1 {max(1,t i on,t -∆ rw t,on ),...,min(Ts,t i on,t +∆ rw t,on )} (t on ), q(t off |t on , t i off,t ) ∝ 1 {max(ton,t i off,t -∆ rw t,off ),...,min(Ts,t i off,t +∆ rw t,off )} (t off ).

(

) 27
where ∆ rw t,on and ∆ rw t,off are adapted using the same strategy as previously described for the other proposals, i.e. for ∆ rw t,on :

∆ rw t,on = N m=1 W (m) t-1 t (m) on,t-1 -µ ton,t-1 2 1 2 (28) with µ ton,t-1 = N m=1 W (m) t-1 t (m) on,t-1 ∆ rw t,off is obtained by just replacing t (m) on,t-1 by t (m)
off,t-1 . This adaptive Metropolis within Gibbs used in the implementation of the SMC sampler through this work is summarized in Algorithm 2.

Sequence of backward kernels L t

The backward kernel L t is arbitrary, however as discussed in (Del [START_REF] Moral | Sequential Monte Carlo samplers[END_REF], it should be optimized with respect to mutation kernel K t to obtain good performance. [START_REF] Moral | Sequential Monte Carlo samplers[END_REF]) establish that the backward kernel which minimizes the variance of the unnormalized importance weights, W t , are given by

L opt t (θ t+1 , θ t ) = η t (θ t )K t+1 (θ t , θ t+1 ) η t+1 (θ t+1 ) (29)
However, it is typically impossible to use these optimal kernels as they rely on the marginalization of the joint distribution defined in Eq. ( 14) which do not admit any closed form expression, especially if an MCMC kernel is used as K t which is π t -invariant distribution. Thus we can either choose to approximate L opt t or choose kernels L t so that the importance weights are easily calculated or have a familiar form. As discussed in (Del [START_REF] Moral | Sequential Monte Carlo samplers[END_REF], if an MCMC kernel is used as forward kernel, the following L t is employed

L t-1 (θ t , θ t-1 ) = π t (θ t-1 )K t (θ t-1 , θ t ) π t (θ t ) (30)
which is a good approximation of the optimal backward if the discrepancy between π t and π t-1 is small; note that (30) is the reversal Markov kernel associated with K t . In this case, the unnormalized incremental weights becomes

w (m) t (θ (m) t-1 , θ (m) t ) = γ t (θ (m) t-1 ) γ t-1 (θ (m) t-1 ) = p(y|θ (m) t-1 ) (φt-φt-1) (31)
This expression ( 31) is remarkably easy to compute and valid regardless of the MCMC kernel adopted. Note that φ tφ t-1 is the step length of the cooling schedule of the likelihood at time t. As we choose this step larger, the discrepancy between π t and π t-1 increases, leading to increase the variance of the importance approximation when it deteriorates. Thus, it is important to construct a smooth sequence of distributions {π t } 0≤t≤T by judicious choice of an associated real sequence {φ t } T t=0 . Let us remark that when such backward kernel is used, the unnormalized incremental weights in Eq. ( 31) at time t does not depend on the particle value at time t but just on the previous particle set. As suggested in (Del [START_REF] Moral | Sequential Monte Carlo samplers[END_REF], in such case, the particles θ Based on this discussion regarding the different choices, the SMC sampler that will be used for Bayesian inference in the STE model is summarized in Algorithm 1.

Algorithm 1 SMC Sampler Algorithm for STE 1: Initialize particle system

2: Sample θ (m) 1 N m=1 ∼ η1(•) and compute W (m) 1 = γ 1 (θ (m) 1 ) η 1 (θ (m) 1 ) N j=1 γ 1 (θ (j) 1 ) η 1 (θ (j) 1 ) -1
and do resampling if ESS < ESS 3: for t = 2, . . . , T do 4:

Set the temperature schedule such that the CESS defined in Eq. ( 23) decays in a regular predefined way.

5:

Computation of the weights:

for each m = 1, . . . , N W (m) t = W (m) t-1 p(y|θ (m) t-1 ) (φ t -φ t-1 )
Normalization of the weights :

W (m) t = W (m) t N j=1 W (j) t -1 6: Selection: if ESS < ESS then Resample 7: Mutation: for each m = 1, . . . , N : Sample θ (m) t ∼ Kt(θ (m) t-1 ; •) where Kt(•; •) is a πt(•) invariant Markov kernel described in more details in Algo. 2. 8: end for
Algorithm 2 Adaptive Metropolis-within-Gibbs Kernel K t (•; •) for the m-th particle

1: Input: Set θ 0 = x 0 s , log(α) 0 , (t 0 on , t 0 off ), log(σ 2 ǫ ) 0 = θ (m) t-1
2: for i = 1, . . . , NMCMC do 3:

MHwGibbs for the position of the source:

4: Sample x * s ∼ N x i-1 s
, Σ rw xs,t with Σ rw xs,t defined in Eq. 25

5:

Compute the Acceptance ratio:

ρx = min 1, p(y|x * s , θ i-1 -xs ) φ t p(x * s ) p(y|x i-1 s , θ i-1 -xs ) φ t p(x i-1 s ) with θ i-1 -xs = {log(α) i-1 , (t i-1 on , t i-1 off ), log(σ 2 ǫ ) i-1
} corresponding to the set containing all the unknown parameters of interest excepting xs.

6:

Set x i s = x * s with probability ρx, otherwise set x i s = x i-1 s 7:
MHwGibbs for the level of emission:

8:
Sample log(α) * ∼ N log(α) i-1 , Σ rw α,t by following the same procedure as for xs 9:

Compute the Acceptance ratio:

ρα = min 1, p(y| log(α) * , θ i-1 -α ) φ t p(log(α) * ) p(y| log(α) i-1 , θ i-1 -α ) φ t p(log(α) i-1 ) with θ i-1 -α = {x i s , (t i-1 on , t i-1 off ), log(σ 2 ǫ ) i-1 } 10:
Set log(α) i = log(α) * with probability ρα, otherwise set log(α) i = log(α) i-1 11:

MHwGibbs for the time of emission:

12:

Sample (t * on , t * off ) ∼ q(ton, t off |t i-1 on , t i-1 off ) defined in Eq. 26

13:

Compute the Acceptance ratio:

ρt = min 1, p(y|t * on , t * off , θ i-1 -t ) φ t p(t * on , t * off )q(t i-1 on , t i-1 off |t * on , t * off ) q(t * on , t * off |t i-1 on , t i-1 off )p(y|x i-1 s , θ i-1 -t ) φ t p(t i-1 on , t i-1 off ) with θ i-1 -t = {x i s , log(α) i , log(σ 2 ǫ ) i-1 } 14:
Set (t i on , t i off ) = (t * on , t * off ) with probability ρt, otherwise set (t i on , t i off ) = (t i-1 on , t i-1 off )

15:

MHwGibbs for the variance of measurements' noise:

16: Sample log(σ 2 ǫ ) * ∼ N log(σ 2 ǫ ) i-1 , Σ rw σ 2 ǫ ,t
by following the same procedure as for xs 17:

Compute the Acceptance ratio:

ρ σ 2 ǫ = min 1, p(y| log(σ 2 ǫ ) * , θ i-1 -σ 2 ǫ ) φ t p(log(σ 2 ǫ ) * ) p(y| log(σ 2 ǫ ) i-1 , θ i-1 -σ 2 ǫ ) φ t p(log(σ 2 ǫ ) i-1 ) with θ i-1 -σ 2 ǫ = {x i s , log(α) i , (t i on , t i off )} 18: Set log(σ 2 ǫ ) i = log(σ 2 ǫ ) * with probability ρ σ 2 ǫ , otherwise set log(σ 2 ǫ ) i = log(σ 2
ǫ ) i-1 19: end for 20: Set the new particle value at time t as θ

(m) t = x N MCMC s , log(α) N MCMC , (t N MCMC on , t N MCMC off ), log(σ 2 ǫ ) N MCMC

Numerical Simulations

The performances of the source term estimate method presented in this paper are assessed using twin experiments assuming a fictitious malevolent release in a complex atmospheric environment. The exact nature of the release is not at stakes for this hypothetical event. It could be constituted of chemical products, biological pathogenic agents, radioactive-nuclear materials or explosives (CBRN-E) In order to extend the study done by [START_REF] Septier | A Bayesian inference procedure based on inverse dispersion modelling for source term estimation in built-up environments[END_REF], we consider the same actual urban area corresponding to the district of the Opera in Paris (France), where all buildings are explicitely accounted for. The 3D simulation domain has dimensions of 1.1km × 0.9km × 1.6km. It is meshed at a horizontal regular resolution of 2 meters and a vertical resolution of 2 meters near the ground decreasing with the elevation above the ground. Figure 1a depicts the simulated source and 20 virtual sensors set up in the urban district, which are represented respectively by a green asterisk and red crosses. The location of the source is chosen here as an example and it could be anywhere in the simulation domain. The sensors are displayed quite regularly in the street network. For the sake of simplicity, both the simulated source and sensors are assumed to be located at the same level close to the ground, immersed in the urban canopy. Even if the source term localization is a two-dimensional problem in these twin experiments, the flow and dispersion simulations are performed in three-dimensional space across the time dimension. The meteorological conditions considered during and after the fake noxious release correspond to a real weather sequence with the wind initially blowing from the west-northwest, then gradually rotating and, finally, blowing from north-northeast. Dispersion simulations are carried out with the model defined in Eq. ( 2) in backward mode from each virtual sensor emitting unit releases each minute over a 45-minute period. Owing to the reciprocity between direct and retrograde atmospheric transport, the backward computations allow us to obtain the source-receptor matrix C(•) in Eq. (3). As discussed in Section 2.1, backward atmospheric transport has the advantage of drastically reducing the number of computations necessary to obtain this matrix as there are many less virtual sensors than potential (N × N MCMC × T ) positions of the emitting source generated by the SMC sampler. More precisely, using the dispersion model in backward mode, N c × T c computations are needed where N c is the number of sensors and T c the number of time samples collected by each sensor. In the test-case presented in the paper, N c and T c are equal respectively to 20 and 45. One can notice that in a practical situation, these numbers would even be lower. For example, results are given for the source represented by the green asterisk in Figure 1a. The simulated concentration values on all sensors affected by the plume are shown as a function of time in Figure 1b with a different color for each sensor. The concentration histories with additional noise are assimilated to pseudo-measurements, which in turn are used by the source term estimate method. For example, realizations of the noisy concentration pseudo-measurements are presented in Figure (1c-1d), where the clipped normal distribution, defined in Eq. ( 5), is used with two different values of the variance for all types of errors (measurements, model, etc) denoted by σ 2 ǫ . The different colors represent different sensors. In this study, the 3D flow illustrated and dispersion simulations have been performed with the Parallel-Micro-SWIFT-SPRAY (PMSS) system presented in [START_REF] Oldrini | Description and preliminary validation of the PMSS fast response parallel atmospheric flow and dispersion solver in complex built-up areas[END_REF] and (Oldrini et al., 2019). PMSS is the efficient parallel version in time, space, and Lagrangian particles of Micro-SWIFT-SPRAY (MSS), which was originally developed to provide simplified CFD solutions for the flow and dispersion, respectively with SWIFT and SPRAY, in built-up environments in a limited amount of time [START_REF] Tinarelli | Description and preliminary validation of the PMSS fast response parallel atmospheric flow and dispersion solver in complex built-up areas[END_REF]. On the one hand, SWIFT is a 3D terrain-following mass-consistent diagnostic model taking account of the buildings and providing the 3D fields of wind, turbulence, temperature, and humidity. The resulting velocity fields of the wind used in all the simulations are depicted in Fig. 2. On the other hand, SPRAY is a 3D Lagrangian Particle Dispersion Model able to account for the presence of buildings. PMSS has been systematically validated against numerous wind tunnel and in-field experimental campaigns for short and prolonged releases [START_REF] Castelli | Validation of a Lagrangian particle dispersion model with wind tunnel and field experiments in urban environment[END_REF]Oldrini and Armand, 2019). Besides the compliance of PMSS results with the statistical acceptance criteria of dispersion models in built-up environments, SPRAY dispersion model can be run in both direct and retrograde mode what is a great advantage to evaluate the source-receptor matrix C(•) in Eq. ( 3).

The performances of the proposed SMC sampler are assessed with N = 200 particles and N MCMC = 10 MCMC moves within the Markov kernel K t (•; •). Figures 3 and 4 give the results obtained by the proposed SMC sampler for a source (represented respectively by a green asterisk in Fig. 1a) located at (300;320) with the variance of the model and measurements errors set to σ 2 ǫ = 10 -6 and σ 2 ǫ = 10 -5 , respectively. Firstly from Figs. 3a and 4a, we can remark that in both cases, the position of the source is correctly estimated since the mean of the posterior distribution of the source position p(x s , y s |y) approximated by the algorithm is close to the ground truth (≤ 10 meters for the noisiest scenario in Fig. 4a). As expected, the uncertainty given by the empirical measure of the posterior distribution for both p(x s |y) and p(y s |y) is higher when the errors are larger (i.e. σ 2 ǫ = 10 -5 ). The remarks holds regarding the estimation of the posterior distributions of σ 2 ǫ and q depicted in Figs ( 3b,4b) and (3c,4c), respectively. Let us remark that the SMC sampler learns only from the the measurements that the errors are larger in Fig. 4 than in Fig. 3 since the variance of the noise is unknown and estimated by the algorithm quite accurately. Having an estimation of such posterior distributions for all the unknown parameters is of high interest for users in order to be able to quantify the associated uncertainty which is not possible when an optimization technique is used to solve this STE problem. Moreover, the resulting posterior distributions can be easily obtained without requiring a complex tuning of parameters by a non-expert user thanks to the proposed complete automatic and adaptive Bayesian solution.

In Figure 5, the proposed SMC sampler is compared to a more traditional MCMC algorithm. The MCMC approach is using the same Adaptive Metropoliswithin-Gibbs Kernel described in Algo. 2 as in the proposed SMC sampler. The MCMC estimator obtained to compute the source position is obtained by taking the last 60 % of the simulated Markov chain values (the first samples are considered to be the burn-in period). Concerning the SMC sampler, the estimator is based only on the last population of N = 200 particles generated at a given iteration. The median of the Squared Error (SE) has been obtained by running 50 independent runs of both algorithms. From the results, the SMC sampler significantly outperforms the corresponding MCMC algorithm based on the same MCMC kernel. Firstly, the SMC sampler has the benefit of using a population of N samples compared to the MCMC in which the samples are generated successively. Additionally, the use of likelihood tempered target sequence clearly facilitates the exploration of the state-space and thus limits the problem of having the inference algorithm to get stuck in one region of the space, which could easily happen in our problem due to the presence of many buildings in the surveillance area. Over the 50 trials, the squared errors of the estimates are lower than 100 (which means that the estimate of the position of the source lies in a circle of 10 meters around the ground truth) only 10 times for the MCMC algorithm but 50 times for the proposed SMC sampler. These results clearly highlight the benefit of using the SMC sampler in this complex environment. This fact also confirms in our context the remark mentioned in [START_REF] Jasra | Inference for lévy-driven stochastic volatility models via adaptive sequential monte carlo[END_REF] stating that the SMC samplers, when using a given MCMC kernel, often out-perform the corresponding MCMC algorithm.

Figures 6 and7 show the results obtained using the SMC sampler for a source located at (300;320) but with two different emission profile, a continuous and a prompt release, which may be due to an insidious long emission of a deleterious substance and an explosion, respectively. In both scenarios, the algorithm is able to provide accurate estimates for all the parameters of interest. By comparing the two scenarios, the results empirically show that the level of posterior uncertainties in the parameters decreases with the length of duration of the release. The larger amount of pollutant in the atmosphere in the case of Fig. 6 due to a longer release allows us to have more useful information from the sensor measurements and therefore more confidence about the parameters of interest.

x-axis in units 8) shows the output of a single run of the proposed SMC sampler in the case of source located at (410;80). As can be seen in Fig ( 8a), it is difficult to distinguish any smooth useful output concentration levels as 395 in Fig. 1b from the measurements due to the relatively large level of noise. Nevertheless, even in this challenging the STE problem, the algorithm is still able to accurately estimate the source location as well as both the release rate and the level of observation noise.

Time in minutes

Conclusion

In this paper, a complete automatic Bayesian solution is proposed to solve the source term estimation problem. In a previous study [START_REF] Septier | A Bayesian inference procedure based on inverse dispersion modelling for source term estimation in built-up environments[END_REF], an adaptive Bayesian solution based on importance sampling was proposed. Unfortunately despite its performance, the use of such algorithm was conditioned on quite severe and not realistic assumptions regarding the model. Indeed in the previous work, both the likelihood and the prior function of the emission rate vector had to be assumed to be normal distributions, thus leading to the possibility of having negative values for some parameters such as the emission rate. Moreover, the normal prior distribution necessitated to be tuned quite finely by the user which could clearly be difficult in practice.

In this paper, we provide a novel framework that overcomes these limitations. More realistic and more general distributions of all physical quantities and parameters can therefore be used and moreover the model parameters are now assumed to be completely unknown and are included in the estimation pro-cedure. Compared to previous algorithms, the proposed strategy provides a full probabilistic solution without having to choose, sometimes arbitrarily, values for certain model parameters, which therefore makes its use easier for practitioners.

The proposed solution is based on recent advances in the field of stochastic simulation algorithms and more precisely on the Sequential Monte Carlo sampler. Although this approach presents many advantages over traditional Monte-Carlo methods, the potential of this emergent technique is still however largely underexploited in practice. Moreover, the SMC sampler has been efficiently designed for this STE problem by integrating some adaptive schemes in order to use a sequence of intermediate distributions using appropriate cooling temperatures and also to explore the state in the MCMC moves using adaptive covariance for each parameter. Numerical twin experiments clearly demonstrate the ability of the proposed solution to accurately estimate the posterior distribution of the parameters of interest in complex atmospheric environments such as industrial and urban ones. As future research directions, we plan to study the impact of the scenario's characteristics such as the environmental variables, source's characteristics, placement and noise of the sensors on the proposed Bayesian solution. Some real datasets, such as the FUSION Fields Trials 2007 (FFT07) [START_REF] Platt | Comparative investigation of source term estimation algorithms using FFT07 data[END_REF] or the MUST campaign Biltoft (2001), will also be used to assess the ability of the proposed SMC sampler to solve the STE problem in a large number of different and various scenarios.

Figure 1 :

 1 Figure 1: Scenario under study -(a) 20 sensors (red) & 1 source (green). (b) Output concentration levels, y i,n at the sensor location from the dispersion model obtained every minute by the 20 sensors from 9:00 to 9:45. A realization of the noisy measurements, y i,n obtained using the clipped normal likelihood defined in Eq. (5) with σ 2 ǫ = 10 -6 (c) and σ 2 ǫ = 10 -5 (d).

  (a) 2:00 p.m. (b) 2:09 p.m. (c) 2:18 p.m. (d) 2:27 p.m. (e) 2:36 p.m. (f ) 2:45 p.m.

Figure 2 :

 2 Figure 2: Wind Velocity fields of the considered scenario at different time instants and 2 meters from the ground.

Figure 3 :

 3 Figure 3: Results of the SMC sampler for STE with a release at location (300;320) and noise variance σ 2 ǫ = 10 -6 in Eq. (5). (a) Representation of the N = 200 particles regarding the source position (blue) along with the ground truth (green) and the estimation of p(xs|y)/p(ys|y) on the bottom/right. (b) Estimation in red of p(σ 2 ǫ |y) and the true value in dashed black. (c): Mean of p(q|y) (black) and the 90% confidence interval (grey) compared to the ground truth (red).

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Results of the SMC sampler for STE with a release at location (300;320) and noise variance σ 2 ǫ = 10 -5 in Eq. (5). (a) Representation of the N = 200 particles regarding the source position (blue) along with the ground truth (green) and the estimation of p(xs|y)/p(ys|y) on the bottom/right. (b) Estimation in red of p(σ 2 ǫ |y) and the true value in dashed black. (c): Mean of p(q|y) (black) and the 90% confidence interval (grey) compared to the ground truth (red).

Figure 7 :

 7 Figure 7: Results of the SMC sampler for STE with a prompt release [variance of the observation uncertainty σ 2 ǫ = 10 -6 in Eq. (5)]. (a) Representation of the N = 200 particles regarding the source position (blue) along with the ground truth (green) and the estimation of p(xs|y)/p(ys|y) on the bottom/right (b) Estimation in red of (a) p(σ 2 ǫ |y) and the true value in dashed black. (c): Mean of p(q|y) (black) and the 90% confidence interval (grey) compared to the ground truth (red).

Finally

Figure 8 :

 8 Figure 8: Results of the SMC sampler for STE with a release at another location (410;80). (a) Measurements obtained with a variance of σ 2 ǫ = 5.10 -5 in Eq. (5)]. (b) Representation of the N = 200 particles regarding the source position (blue) along with the ground truth (green) and the estimation of p(xs|y)/p(ys|y) on the bottom/right. (c) Estimation in red of p(σ 2 ǫ |y) and the true value in dashed black. (d): Mean of p(q|y) (black) and the 90% confidence interval (grey) compared to the ground truth (red).


Nathalie Lassau
email: nathalie.lassau@gustaveroussy.fr

Imad Bousaid

Emilie Chouzenoux

Antoine Verdon

Corinne Balleyguier

François Bidault

Elie Mousseaux

Sana Harguem-Zayani

Loic Gaillandre

Zoubir Bensalah

Isabelle Doutriaux-Dumoulin

Michèle Monroc

Audrey Haquin

Luc Ceugnart

Florence Bachelle

Mathilde Charlot

Isabelle Thomassin-Naggara

Tiphaine Fourquet

Héloise Dapvril

Joseph Orabona

Foucauld Chamming's Q Michael

El Haik

Jules Zhang-Yin

Marc-Samir Guillot

Mickaël Ohana

Thomas Caramella

Yann Diascorn

Jean-Yves Airaud

Philippe Cuingnet

Umit Gencer

Littisha Lawrance

Alain Luciani

Anne Cotten

Jean-François Meder

Three artificial intelligence data challenges based on CT and ultrasound

Keywords: Artificial intelligence, Ultrasonography, Computed tomography, Radiology, Data management Abbreviations 2D: Two dimensional, 3D: Three dimensional, AI: Artificial intelligence, AUC: area under the curve, AUROC: area under the receiver operating characteristic curve, BI-RADS: Breast imaging-reporting and data system, C-index: Concordance index, CT: Computed tomography, DECT: Dual-energy computed tomography, ENT: Ear nose and throat, GDPR: General Data Protection and Regulation, Magnetic resonance imaging, ROC: Receiver operating characteristic, SFR:

The 2020 edition of these Data Challenges was organized by the French Society of Radiology (SFR), from September 28 to September 30, 2020. The goals were to propose innovative artificial intelligence solutions for the current relevant problems in radiology and to build a large database of multimodal medical images of ultrasound and computed tomography (CT) on these subjects from several French radiology centers. This year the attempt was to create data challenge objectives relative to the clinical routine of radiologists, with less preprocessing of data and annotation, leaving a large part of the preprocessing task to the participating teams. The objectives were proposed by the different organizations depending on their core areas of expertise. A dedicated platform was used to upload the medical image data, to automatically anonymize the uploaded data. Three challenges were proposed including classification of benign or malignant breast nodules on ultrasound examinations, detection and contouring of pathological neck lymph nodes from cervical CT examinations and classification of calcium score on coronary calcifications from thoracic CT examinations.

A total of 2076 medical examinations were included in the database for the three challenges, in three months, by 18 different centers, of which 12% were excluded. The 39 participants were divided into six multidisciplinary teams among which the coronary calcification score challenge was solved with a concordance index > 95%, and the other two with scores of 67% (breast nodule classification) and 63% (neck lymph node calcifications).

Introduction

In recent years, artificial intelligence (AI) has become increasingly important in medical research [START_REF] Syeda-Mahmood | Role of big data and machine learning in diagnostic decision support in radiology[END_REF]. Indeed, the continuous creation of medical data makes machine learning and deep learning more effective, helping radiologists in their work [START_REF] Recht | Artificial intelligence: threat or boon to radiologists?[END_REF][START_REF] Thrall | Artificial intelligence and machine learning in radiology: opportunities, challenges, pitfalls, and criteria for success[END_REF], either in complex tasks like classification and detection or repetitive segmentation tasks [START_REF] Litjens | A survey on deep learning in medical image analysis[END_REF][START_REF] Society | What the radiologist should know about artificial intelligence -an ESR white paper[END_REF]. These new technologies represent an important advancement not only in terms of diagnostic capabilities but also in terms of quality of work life for radiologists [START_REF] Thrall | Artificial intelligence and machine learning in radiology: opportunities, challenges, pitfalls, and criteria for success[END_REF]. This requires a change in work habits in the health community, as it is important to involve patients, radiologists and AI researchers to work together, while abiding by the new regulations of the General data Protection and Regulation (GDPR) [START_REF] Liew | The future of radiology augmented with artificial intelligence: a strategy for success[END_REF][START_REF] Pesapane | Artificial intelligence as a medical device in radiology: ethical and regulatory issues in Europe and the United States[END_REF]. The first two editions of this Data Challenge were held during the Journées Francophones de Radiologie in 2018 and 2019. Its main goals were: (i), to help radiologists answer difficult issues in their field; (ii), create large medical databases that can be reused in the future; and (iii), gather and challenge multidisciplinary teams to raise awareness of new challenges among all the persons involved in AI and radiology research.

Several publications have resulted from prior Data Challenges [START_REF] Lassau | Five simultaneous artificial intelligence data challenges on ultrasound, CT, and MRI[END_REF][START_REF] Lassau | Three artificial intelligence data challenges based on CT and MRI[END_REF], and also on related subjects for which the work was initiated from these challenges, such as the works on meniscus tear detection on two-dimensional (2D) magnetic resonance imaging (MRI) [START_REF] Roblot | Artificial intelligence to diagnose meniscus tears on MRI[END_REF][START_REF] Couteaux | Automatic knee meniscus tear detection and orientation classification with Mask-RCNN[END_REF], thyroid cartilage classification on 2D computed tomography (CT) [START_REF] Santin | Detecting abnormal thyroid cartilages on CT using deep learning[END_REF], breast lesion characterization on 2D MRI [START_REF] Herent | Detection and characterization of MRI breast lesions using deep learning[END_REF], liver lesion classification on 2D ultrasound [START_REF] Schmauch | Diagnosis of focal liver lesions from ultrasound using deep learning[END_REF], kidney cortex segmentation using 2D CT [START_REF] Couteaux | Kidney cortex segmentation in 2D CT with U-Nets ensemble aggregation[END_REF], sarcopenia assessment on 2D CT [START_REF] Blanc-Durand | Abdominal musculature segmentation and surface prediction from CT using deep learning for sarcopenia assessment[END_REF], prediction of clinical disability in patients with multiple sclerosis using fluid attenuated inversion recovery MRI [START_REF] Roca | Artificial intelligence to predict clinical disability in patients with multiple sclerosis using FLAIR MRI[END_REF] and classification of pulmonary CT [START_REF] Blanc | Artificial intelligence solution to classify pulmonary nodules on CT[END_REF] are a few examples. The third edition of these challenges took place in September 2020. The objective was to answer the issues faced in clinical routine with limited preprocessing performed on the data by the radiologists. This made the objectives more difficult to reach, but made the AI more efficient. The data protection was optimized to be more compliant with the GDPR standards, by automatically anonymizing while uploading the medical examinations.

Material and methods

Clinical objectives

Radiologists from several organizations of the French Society of Radiology (SFR) proposed the objectives, in order to cover all the different domains, related to the current reality of medical diagnosis. These objectives were discussed with data scientists of the National Research Institute in Digital Sciences and Technologies (INRIA, Saclay, France) in order to verify the comprehensiveness of the objectives by AI, and the relevance of using the AI instead of another technology. This first step was essential in order to determine the important issues in both the medical and AI fields. A review of the scientific literature was done to verify if the proposed AI projects were not already answered in the publications, on Kaggle and on Grand-Challenge websites. Finally, a feasibility assessment was done to verify the collection of a sufficient number of medical examinations for the different challenges. The final dataset contained about 600 medical image series for each challenge.

Security and data protection

All the stages of the project, from data uploading to the challenge phase, were designed to be aligned with the GDPR. The Commission Nationale Informatique et Libertés was consulted, and the SFR assumed the role of Data Protection Officer. Automatic anonymization of the examinations was ensured by the platform developed by the startup Cleverdoc Ⓒ . This process ensured the protection of patients' data, before the participants had access to it. The DICOM metadata on the age, gender, image size, modality, and pixel spacing was retained post anonymization. In addition, every radiologist involved was asked to send an information letter to patients about the use of their medical examinations, with the option of refusal of consent.

A data chart was also sent to each radiologist to help them abide by the GDPR rules, as well as guidelines on the terms of use of the data for the participants to the challenge.

Communication and uploading

The data uploading phase began on June 19, 2020, however the first step was to contact as many radiologists as possible. The SFR contacted its network of radiologists. The interested radiologists had to register on the platform, and they were welcomed by email with a user guide, the technical and clinical specifications of the question and the format of the medical data to be uploaded was sent to them. Our aim was to have an all-in-one interface through the platform for the radiologists, to upload and annotate with a uniform format throughout the dataset. They were trained through webinars and provided with a support contact email and phone numbers. A follow-up of their progress was also done. In addition, a weekly update was sent to all radiologists that included an overall report on the project's progress. The purpose of this newsletter was to motivate radiologists to participate actively in the project and for them to add the required medical examinations. The uploaded data was monitored on a daily basis to check the conformity of the medical image examinations and an expert checked the whole datasets before sending it to the participants.

Team gathering and challenge phases

Each team had to have multidisciplinary member background: at least one radiologist, one engineer/data-scientist and an engineering/PhD student. The team could be a startup, a big company or a research laboratory. Each team would have to register on the platform with all the members' details. Three datasets were sent to them for each challenge: the first dataset batch was sent on September 3, the second on September 28, for the teams to train their algorithms. Finally, the validation dataset was sent during the data challenge on September 30 at 1PM, and teams had until 2 PM to submit their results. The jury calculated the scores of each team, with the help of an INRIA researcher, and the winners were announced on the next day, on October 1. A prize of €3000 was awarded to the winner of each challenge.

Results

There

Breast nodule classification on ultrasound

The radiologists use the Breast Imaging-Reporting and Data System (BI-RADS) classification for mammography, ultrasound and MRI examinations to define the abnormalities seen and to indicate clinical management, (i.e., return to screening, close follow-up or biopsy). The purpose of this question was to determine the risk of malignancy of the findings to establish appropriate diagnostic strategy for the patient. It was challenging for even the human eye to categorize some. Thus, there was a need to have a more powerful algorithm for the particular abnormalities in question. The test dataset composed of ultrasound images of the breast with abnormal findings, was annotated as "benin" or "malignant" depending on the corresponding histological proof.

Pathological lymph nodes classification on cervical CT examinations

The presence of abnormal lymph nodes is an important finding in the interpretation of ear nose and throat (ENT) CT examinations in patients with infection or cancer. The presence of cervical lymph nodes has an initial prognostic impact in ENT cancer. CT is widely used for a variety of clinical and therapeutic situations involving ENT adenopathy. The challenge was therefore to detect and assess the pathological lymph nodes on cervical dual energy (DE) CT examinations. The test dataset composed of cervical three-dimensional (3D) DECT examinations with normal and/or pathological lymph nodes. The locations of some specific lymph nodes were annotated on the images, with pathological or normal assignments. All nodes were not annotated as there were more than 70 lymph nodes present per patient, in which the amount of normal lymph nodes is much more than the pathological ones.

Coronary calcifications classification on thoracic CT examinations

Coronary calcification is responsible for many cardiovascular diseases and its early detection is an important prevention tool. When calcium score is low, the patient is followed up to limit his or her risk of disease. The challenge was to determine the calcification class, based on Agatston's calcium score for a score of 0 as Class A, score between 1 and 9 as Class B, between 10 and 99 as Class C, between 100 and 399 as Class D, for a score greater than 400, classified as class E. The test dataset was composed of thoracic 3D CT examinations with or without coronary calcification. The location of the calcifications was annotated on the images.

Score computation

Breast nodule classification

The final score was calculated using area under the receiver operating characteristic curve (AUROC) method. The final score ratio of the AUROC between benign / malignant was defined using the following equation:

Final score = AUROCC (Benign/Malignant)

Pathological lymph nodes

For each patient's examination, a different score was considered depending on the patient's condition. If the individual was sick the score was the Dice score of the area of the pathological nodes, else if the individual was healthy, the score was defined using the following equation:

Score = ,

where M was the mean of the surface area of pathological nodes on the examinations of sick patients and A was the surface area of pathological nodes erroneously detected by the team.

There was a need to define an alternative to the Dice score due to some examinations with no pathological nodes. For these images, the Dice sore was not representative. The teams were evaluated by penalizing them if they detected pathological nodes on these examinations. The bias arising due to this evaluation method was minimized by including the variable M in the calculations.

Coronary calcifications

The final score was defined as the concordance index (C-index) score that compared the classification of examinations based on the Agatston's calcification score. It was basically similar to the area under the curve (AUC) analysis, but allowed a score to be established for multi-class predictions by establishing an order among classes on a scale (class A to class E).

At the beginning, the mean AUC on the five classes was opted, but the C-index score allowed an order of the classes, which was more relevant for this challenge because classes A to E were defined as the successive stages of coronary calcification.

Communication and team gathering

Six among the eight registered teams participated in the data challenge (BrAIn, CS/EURECOM; GAMC, Owkin, Philips, Radioadvisor). Thirty-nine team members participated in this challenge, including six researchers, 13 students, 12 employees, and eight radiologists, in six different teams. The upload phase began on June 19, and the first images were uploaded by mid-July. A total of 2,085 images were uploaded from 18 radiology services (15 public hospitals, three profit private clinics, and three cancer centers) with 681 for the breast nodules, 737 for the pathological lymph nodes and 667 for the coronary calcifications (Table 1).

Data processing

Only a few preprocessing operations were performed on the data, such as automatic anonymization of the examinations and a pixel slope. The aim was to preserve the data as close to the ones used in reality by not filtering the image noise, as these AI models were intended to be used in routine in medical imaging without being pre-processed. Therefore, this strategy was adopted from the point of choosing the clinical objectives till the end of the challenge.

An important issue faced was in the gathering of data from multiple centers with the same characteristics, for the AI to train on a homogeneous format. Each challenge had its own characteristic format. For the Breast Nodules Challenge, data consisted of a B-mode ultrasound image in JPEG format. For the Pathological Lymph Nodes Challenge, data were extracted from 3D DECT examinations obtained after intravenous administration of iodinated contrast material with the patient in supine position in DICOM format with a soft part filter.

For Coronary Calcification Challenge, data were extracted from 3D CT examinations with a slice thickness of 2.5-or 3 mm without gaps and iteratively reconstructed without filters.

When imposing these characteristics on radiologists, some poor quality data were excluded such as cervical CT examinations containing metal artifacts. Some images were also not suitable due to features that would have disturbed the training of the AI models, such as calipers within the lesion. However, some noisy images were not removed if the noise was not in the area of interest. The teams had to avoid these details to focus on the area of interest.

Therefore, the medical images were sorted to have relevant datasets, with the same format with homogeneous distribution between classes and quality (Tables 2, 3 and4). On October 1st, 2020, the three winners of the data challenges were announced based on the results calculated by all the teams. Philips team obtained the best results for Breast Nodule Classification Challenge (AUROC = 0.666 compared to 0.624 for Owkin and 0.643 for Radioadvisor) and GAMC team for Pathological Lymph Node Challenge (score = 0.631) and

Coronary Calcification Challenge (C-index = 0.951 compared to 0.909 for Owkin).

Discussion

The Data Challenge 2020 was designed to resemble routine clinical examinations of the patients. The data challenge in itself was an innovation because the AI models produced by the teams were useful for the radiologists. The aim was, for the AI to be effective on real-time medical examinations. In addition, the objectives concerned the human capacities of the radiologists, to classify the benign/malignant breast nodules and also the tedious task of contouring pathological neck lymph nodes on 3D examinations represented as a very relevant task for radiologists and as a complex task for an AI. The technical difficulty was therefore present on both detection and classification objectives. The optimization of uploading also presented an improvement in the method to organize the challenge. Gathering data from different medical centers across France was a complex task. A dedicated platform for the recuperation of data was created the previous year, but this year it was outsourced it to a startup Cleverdoc Ⓒ , which specializes in medical data collection and AI models, to propose a more optimized solution to facilitate uploading and annotation. Thus, it was simpler and faster for radiologists, and made it possible to automate their registration, uploading and anonymization of examinations. Such a platform was a good solution for the future as it ensured simplicity and being completely compliant with the GDPR rules.

The results of the Data Challenge 2020 were quite disappointing for the Breast Nodules Nevertheless, the score used presented a certain bias, and it would have been better to have an alternative solution in further Data Challenges, to limit this bias and to have an objective score.

The BI-RADS classification of breast nodules gives a probability of malignancy of a nodule and thus helps radiologists diagnose benign or malignant nodules. In particular, nodule category helps decide to perform a biopsy and reduce the rate of useless biopsies of benign nodules. This classification was initially made by clinical observations on mammography, breast ultrasound, and breast MRI [START_REF] Lee | BIRADS 3: Current and future use of probably benign[END_REF], but no study has been able to determine precisely the malignancy of these nodules [START_REF] Melnikow | Supplemental screening for breast cancer in women with dense breasts: a systematic review for the U.S. Preventive Service Task Force[END_REF]. Through the creation of this classification, numerous shape parameters are used on different image modalities such as "skin", "vascular", "coarse or popcorn-like", "large rod-like" [START_REF] Spak | BI-RADS fifth edition: A summary of changes[END_REF]. The Breast Nodule Classification Challenge was focused on ultrasound, one single image, and some patient data (age, gender), which represented a major challenge and also a significant step forward for clinicians [START_REF] Ghosh | Artificial intelligence using open source BI-RADS data exemplifying potential future use[END_REF]. Three machine learning algorithms were used for this challenge and the results showed that AI remains significantly less effective alone than with the help of a radiologist. Therefore, as anticipated, results were quite disappointing. Nevertheless, the development of deep learning algorithms for breast nodule classification should be encouraged as computer-aided diagnosis tools [START_REF] Morgan | Applications of artificial intelligence in breast imaging[END_REF].

Presence of pathological lymph nodes is an important finding in patients with ENT infection or cancer. The presence of cervical lymph nodes has an initial prognostic impact in ENT cancers taking into account by the Tumor Node Metastasis classification [START_REF]TNM classification of malignant tumours[END_REF]. In addition, lymph node is a frequent form of recurrence of ENT cancers. A meta-analysis gives a precise definition of the radiological semiology of ENT adenopathy as well as the performances of the semiological signs and imaging techniques for the diagnosis [START_REF] Bondt | Detection of lymph node metastases in head and neck cancer: a meta-analysis comparing US, USgFNAC, CT and MR imaging[END_REF]. There are only a few articles on the detection or contouring of lymph nodes whereas studies have used texture analysis from DECT data to classify lymph nodes into pathological or benign ones [START_REF] Seidler | Dual-energy CT texture analysis with machine learning for the evaluation and characterization of cervical lymphadenopathy[END_REF]. For the Pathological Lymph Nodes challenge, machine learning algorithms were used on DECT data, and on a rather limited database (50 patients with 412 lymph nodes). The objective was to detect and contour lymph nodes from 3D DECT examinations, which was challenging.

However, GAMC succeeded to build a solution yielding a Dice score of 0.631.

In conclusion, three Data Challenges with over 2000 GDPR-compliant, multi-centric, 3D CT and ultrasound databases were organized for 6 multidisciplinary teams. For the future challenges, a single platform to manage the whole data challenge workflow is recommended.

It could be also stimulating to work on one or two objectives but to increase the number of data for each of them.

Figures & tables

Table 1. Participating centers that provided imaging data for the three Data Challenges.

  were three main challenges for the Data Challenge 2020: The classification of benign or malignant breast nodules present from ultrasound examinations proposed by the Société d'Imagerie de la Femme (SIFEM; French Women's Imaging Society), the detection and contouring of pathological neck lymph nodes from cervical CT examinations, proposed by the College d'Imagerie pour la Recherche et l'Enseignement en Otorhino Larygologie (CIREOL, Francophone Society of Head and Neck Imaging), and the classification of calcium score on coronary calcifications from thoracic CT examinations proposed by the Société Française d'Imagerie Cardiaque et Vasculaire (SFICV, French Society of Cardiac and Vascular Imaging).

  Challenge and the Pathological Lymph Nodes challenge. The scoring system was also a matter of debate. For the Breast Nodules Challenge, AUC is a commonly used score for binary prediction and was relevant. However, for the Coronary Calcification Challenge, which was a non-binary classification problem, an initial option was to take a mean AUC over the five classes. Since these classes are ordered by Agatston's score (AA < AB < AC < AD < AE), it seemed preferable to use a score based on the C-index, as the values on which they were ordered was more important than identifying the right class. Finally, for the Pathological Lymph Nodes challenge, which combined detection and segmentation, selection of a suitable score was difficult considering that the Dice score did not take into account clinical findings for patients with no pathological lymph nodes. Thus, an alternative score for those examinations had to be determined. Having two different scores for two different types of examinations implied a bias in the final score, which was mitigated by using a factor related to the mean of the pathological lymph node surfaces, such that the two scores were close.

Table 2 .

 2 Examinations received and included for each challenge.

Table 3 .

 3 Training, validation, and testing dataset for each challenge.

Table 4 .

 4 Classification of data for each challenge.

Acknowledgements

All authors thank the Société Française de Radiologie for the opportunity to organize this challenge and its support and Cleverdoc Ⓒ for the platform and for the data hosting.

Funding

This work did not receive any grant from funding agencies in the public, commercial, or notfor-profit sectors.

Human and animal rights

The authors declare that the work described has been carried out in accordance with the Declaration of Helsinki of the World Medical Association revised in 2013 for experiments involving humans as well as in accordance with the EU Directive 2010/63/EU for animal experiments.

Informed consent and patient details

The authors declare that they obtained a written informed consent from the patients and/or volunteers included in the article. The authors also confirm that the personal details of the patients and/or volunteers have been removed.

All the stages of the project, from data uploading to the challenge phase, was designed to be aligned with the GDPR. The Commission Nationale Informatique et Libertés was consulted, and the SFR assumed the role of Data Protection Officer. Automatic anonymization of the examinations was ensured by the platform developed by the startup Cleverdoc Ⓒ . Therefore, the image data was anonymized during the upload, if not anonymized before such that the participants received anonymized medical data. This process ensured the protection of patients' data. In addition, every radiologist involved was asked to send an information letter to patients about the use of their medical examinations, with the option of refusal of consent.

A data chart was also sent to each radiologist to abide by the GDPR laws, as well as guidelines on the terms of use of the data for the participants to the challenge.

Disclosure of interest

The authors declare that they have no competing interest.

Author contribution

All authors attest that they meet the current International Committee of Medical Journal Editors (ICMJE) criteria for Authorship.