Alexandre Lepoutre

Olivier Rabaste

Franc ¸ois Le

Multitarget likelihood for Track-Before-Detect applications with amplitude fluctuations

Track-Before-Detect methods jointly detect and track one or several targets from raw sensor measurements. They often require the computation of the measurement likelihood conditionally to the hidden state that depends on the complex amplitudes of the targets. Since these amplitudes are unknown and fluctuate over time this likelihood must be marginalized over the complex amplitude (i.e. phase and modulus). It has been demonstrated in [1] that this marginalization can be done analytically over the phase in the monotarget case. In this article, we first propose to extend the marginalization to the modulus in a monotarget setting, and we show that closed-forms can be obtained for fluctuations of type Swerling 1 and 3. Second, we demonstrate that, in a multitarget setting, a closed-form can be obtained for the Swerling 1 case. For Swerling 0 and 3 models, we propose some approximation to alleviate the computation. Since many articles consider the case of squared modulus measurements, we also consider this specific case in mono and multitarget settings with Swerling 0, 1 and 3 fluctuations. Finally, we compare the performance in estimation and detection for the different cases studied and we show the gain, both in detection and estimation, of the complex measurement method over the squared modulus method, for any fluctuation model.

I. INTRODUCTION

Filtering consists in recursively finding, at each time step k, the best estimate of an hidden state X k from noisy observations. This is a very general problem encountered in many fields such as econometrics, speech recognition and many others [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF]. In this article we focus more specifically on the surveillance tracking issue that aims at solving two problems:

• detecting over time the appearance and disappearance of one or several targets in the surveillance area covered by a sensor, for instance a radar, a sonar or an infrared (IR) sensor. This is a detection problem.

• estimating over time as precisely as possible target parameters, for instance position, velocity, bearing and so on, from all the observations provided by the sensor until time step k. This is an estimation problem.

In classic radar processing, the measurement data considered at each time step is a set of detection "hits" obtained by thresholding the raw sensor data after matched filtering. Each "hit" either corresponds to a target or to a false alarm due to clutter or sensor noise. Then, several strategies proposed in the literature can be applied to solve the detection and estimation problem. The MHT (Multiple Hypothesis Tracker) algorithm [START_REF] Blackman | Multiple hypothesis tracking for multiple target tracking[END_REF], [START_REF] Reid | An algorithm for tracking multiple targets[END_REF] tries to solve the two problems jointly by considering all possible association hypotheses between hits and tracks. Some other solutions solve the two problems separately: one specific algorithm is devoted to track initiation and termination [START_REF] Blackman | Multiple-target tracking with radar applications[END_REF] December 2, 2014 DRAFT while another algorithm performs the tracking itself, assuming that the number of targets is known and solving the hits/tracks association problem if necessary. Among all algorithms proposed to solve this tracking and association problem stands the J/PDAF (Joint/Probability Data Association Filter) algorithms [START_REF] Fortmann | Sonar tracking of multiple targets using joint probabilistic data association[END_REF], [START_REF] Bar-Shalom | The probabilistic data association filter[END_REF].

Whereas these two strategies are different, they share common features. First, the link between the detection "hit" and the target state is often quite simple, so that it is possible to use Kalman or Extended Kalman filter to estimate the target states. Second, both solutions try to solve the association problem that requires to consider all possible associations between "hits" and tracks. When the number of tracks and the number of "hits" are important, the number of possible associations becomes large and often prohibitive to manage. A possible solution to limit the computational cost consists in reducing the number of "hits" by choosing a high detection threshold. However, miss detections will then often occur for low-SNR (Signal to Noise Ratio) targets, making such targets more difficult to track as they will not provide a detection "hit" at each iteration. Thus, a trade-off must be done between the reduction of the number of associations and the detection and tracking capabilities of the system for low-SNR targets.

To overcome these limitations, a different strategy, known as Track-Before-Detect (TBD, has been proposed in the past fifteen years [START_REF] Salmond | A particle filter for track-before-detect[END_REF], [START_REF] Kreucher | Multitarget Tracking using the Joint Multitarget Probability Density[END_REF], [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF]. Contrary to classic techniques that work on detection "hits", TBD methods directly work on raw sensor data. This allows preserving all the information provided by the data and managing jointly target detection and tracking. As all the available information is kept, one can expect detecting and tracking low SNR targets. Nevertheless, the nature of the measurement data is totally different from classic methods. Indeed, considering data after reception processing and matched filtering, the raw sensor measurement consists of a large multidimensional array where each cell value is provided by the corresponding matched filter output. Thus, two problems arise. The first issue is the nature of the measurement. Indeed the highly non-linear link between the observation and the hidden state as well as the large size of the data array do not allow to use techniques based on Kalman Filter. Therefore, the Bayesian filter cannot be solved exactly and one must resort to some approximation.

Popular solutions are particle filters [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF], [START_REF] Ristic | Beyond the Kalman filter. Particle filters for tracking applications[END_REF] and grid-based methods [1]. Particle filter techniques were first developed for TBD applications in the mono-target case in [START_REF] Salmond | A particle filter for track-before-detect[END_REF] and in the multi-target case in [START_REF] Kreucher | Multitarget Tracking using the Joint Multitarget Probability Density[END_REF], [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF]. The second issue is the computation of the measurement likelihood conditionally to the hidden state. In radar applications considered in this article, the measurement equation depends on the complex amplitudes of the targets which are temporally incoherent, i.e. they may fluctuate independently from measurement to measurement. The amplitude fluctuation is modelled by a uniform distribution for the phase and a Swerling model for the modulus [START_REF] Skolnik | Introduction to Radar Systems /2nd Edition[END_REF]. Note that in another specific framework such as optics, the signal would be real and only the magnitude of the target would be considered. Because of the temporal fluctuation, no information can be gathered over time to estimate the amplitude sequentially contrary to other state parameters. Thus, it is not relevant to compute the likelihood conditionally to these parameters and techniques must be found to compute the likelihood without their knowledge.

The objective of this article is therefore to compute the measurement likelihood in a general multitarget TBD context with unknown fluctuating amplitude parameters. Several solutions have been provided in the literature, mainly in a monotarget setting. The first solution that deals with the unknown complex amplitude considers a December 2, 2014 DRAFT mono-target setting and consists in working on the squared modulus of the complex signal [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF], [START_REF] Salmond | A particle filter for track-before-detect[END_REF], [START_REF] Ristic | Beyond the Kalman filter. Particle filters for tracking applications[END_REF], [START_REF] Rutten | Recursive track-before-detect with target amplitude fluctuations[END_REF], [START_REF] Davey | A comparison of detection performance for several track-before-detect algorithms[END_REF]. For such a radical solution that completely discards the phase dependency, two strategies can be considered to deal with the modulus fluctuation. The first one consists in marginalizing the whole likelihood with respect to the density of the modulus fluctuation [START_REF] Davey | A comparison of detection performance for several track-before-detect algorithms[END_REF]. In practice, this leads to intractable integrals that must be approximated numerically. The second strategy consists in marginalizing independently the likelihood in each cell [START_REF] Rutten | Recursive track-before-detect with target amplitude fluctuations[END_REF]. The advantage of this heuristic second solution is that a closed form can be obtained for fluctuations of type Swerling 0, 1 and 3 [16]. On the other hand the spatial coherence of the modulus, i.e. the fact that the modulus of the target amplitude takes the same value in all cells, is then lost, inducing a possible degradation of performance. Note also that the spatial coherence of the phase is lost for both strategies. This loss was shown in [1] to severely degrade the performance.

To avoid this last drawback, Davey et al. [1] have proposed a new strategy that allows preserving the spatial coherence of the phase. Their solution consists in directly working on complex measurement and marginalizing the complex likelihood of the whole data over the phase. It provides better performance than solutions based on squared modulus. However, they mainly investigated the case where the modulus does not fluctuate (i.e. Swerling 0 case); for modulus fluctuations, they only provide a general marginalization formula. One of the contributions of this article is an extension of their work with complex measurements to fluctuations of type Swerling 1 and 3; we show that closed-forms can be obtained for the monotarget likelihood in both cases.

Furthermore, all the previously discussed strategies with squared modulus or complex measurements were proposed in a monotarget setting. In fact, to our knowledge, the multitarget case has not been investigated in the literature, except for the Swerling 1 case with squared modulus of the measurement [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF]. Thus, another contribution of this article consists in considering the multitarget case both with squared modulus and complex measurements. In the complex measurement case, we provide a closed-form expression for the likelihood in the Swerling 1 case, and we propose in the other fluctuation cases some approximations to alleviate the computational cost. In the squared modulus case, we show that, as soon as at least two targets are present, all phase dependencies cannot be removed from the likelihood; in fact taking the squared modulus permits to remove only one phase, so that other phases must be marginalized. In that latter case, we also propose some approximations in order to reduce the computational complexity.

Overall, this paper intends to provide a generic framework for computing the measurement likelihood in TBD applications in the presence of multiple targets presenting Swerling amplitude fluctuations of type 0, 1 and 3.

Both complex and squared modulus measurements are considered. Closed-form expressions are provided whenever they are obtainable, and approximations are proposed otherwise. We believe that this article thus represents an exhaustive overview for the computation of the measurement likelihood in TBD, summarized in This paper is organized as follows. In section II we present the state and measurement models. Then in section III we present solutions for the likelihood computation from complex and squared modulus measurements. In section IV we derive, when possible, closed forms for the likelihood with Swerling fluctuations of type 0, 1 and 3 in the monotarget and multitarget cases; when not possible, we propose approximations to alleviate the computational time. Finally in section V we present simulation results that show the gain both in detection and in estimation of the complex measurement method over the squared modulus method.

II. PROBLEM FORMULATION

In this section we present the framework for tracking multiple targets from complex measurements in a Track-Before-Detect context.

A. Multi-target tracking

Let us denote by X k the multitarget state at a given time step k. Assuming that N k targets are present at this time step k, X k can be decomposed as

X k = x T k,1 , . . . , x T k,N k T ,

B. Track-Before-Detect Measurement model

In a Track-Before-Detect setting, raw data are directly provided to the tracking filter. Radar raw data for instance consists of range, Doppler and/or 1D or 2D angle measurements. Since matched filtering operations for these parameters are linear, measurements can be considered indifferently before or after matched filtering steps. As the latter approach generally leads to reduced vector sizes and reduced computational cost [START_REF] Rabaste | Long-time coherent integration for low SNR target via particle filter in track-before-detect[END_REF], we will consider here observations after matched filtering along all measurement dimensions. At step k, the measurement z k consists of a vector of size N c corresponding to all measurement cells. It can be expressed by the following non-linear equation that depends on the multi-target state X k :

z k = N k i=1 ρ k,i e jϕ k,i h(x k,i ) + n k , (1)
where:

• h(x k,i ) represents the possibly multidimensional ambiguity function of the i th target centered on the target state x k,i . For the sake of simplicity, h(x k,i ) will be denoted h k,i in the rest of the article.

December 2, 2014 DRAFT • n k is a zero mean circular complex Gaussian vector with covariance matrix Γ.

• ϕ k,i and ρ k,i are respectively the phase and the modulus of the i th target complex amplitude. All variables ϕ k,1:N k and ρ k,1:N k are supposed mutually independent, and independent from n k .

Each phase ϕ k,i is supposed to be unknown and uniformly distributed over the interval [0, 2π) at each time step k.

Concerning the modulus ρ k,i , we consider in this article Swerling fluctuations [START_REF] Skolnik | Introduction to Radar Systems /2nd Edition[END_REF] of type 0, 1 and 3, i.e. constant modulus for the Swerling 0 model, or slow fluctuations for the Swerling 1 and 3 models. These models will be detailed in section IV.

An important point to be stressed here is that variables ρ k,1:N k and ϕ k,1:N k are spatially coherent: this means that the complex target amplitude ρ k,i e jϕ k,j is identical over all cells where the signal ambiguity function spreads.

Taking into account this information can really increase the performance of the Track-Before-Detect algorithms [1].

On the contrary, these variables ρ k,1:N k and ϕ k,1:N k are not assumed coherent over time, i.e. from time sample k to k + 1, amplitude parameters fluctuate independently. As a consequence, no information can be gathered over time on these parameters. These dependencies are illustrated in Figure 1.

III. LIKELIHOOD COMPUTATION

A. Likelihood computation with complex measurements

Bayesian Track-Before-Detect algorithms require the computation of the likelihood p (z k | X k ), i.e. the likelihood of the observation conditionally to the target states. However the measurement equation (

p (z k | X k , ρ k,1:N k , ϕ k,1:N k ) = 1 π N det (Γ) exp -(z k -µ k ) H Γ -1 (z k -µ k ) . ( 2
)
Then by developing (2) we obtain the following expression:

p (z k | X k , ρ k,1:N k , ϕ k,1:N k ) = exp -z k Γ -1 z k π N det (Γ) ×exp - N k i=1 ρ 2 k,i h H k,i Γ -1 h k,i + N k i=1 2ρ k,i |h H k,i Γ -1 z k | cos (ϕ k,i -ξ k,i ) - N k i=1 N k l=i+1 2ρ k,i ρ k,l h H k,i Γ -1 h k,l cos (ϕ k,i -ϕ k,l -φ k,il ) , (3)
where

ξ k,i = arg h H k,i Γ -1 z k and φ k,il = arg h H k,i Γ -1 h k,l . Finally the likelihood p(z k | X k ) can be obtained by marginalization over parameters ρ k,1:N k and ϕ k,1:N k : p(z k | X k ) = • • • R N k ≥0 ×[0,2π) N k p (z k | X k , ρ k,1:N k , ϕ k,1:N k ) p (ϕ k,1:N k ) p (ρ k,1:N k ) dϕ k,1:N k dρ k,1:N k . ( 4
)
First, notice that the spatial coherence is preserved in this formulation thanks to the marginalization. Then, note that the likelihood in Bayesian TBD algorithms such as particle filters can be computed up to a constant, so the constant term

1 π N det (Γ) exp -z k Γ -1 z k
will be discarded in the rest of the paper, and the computation of the likelihood will be always provided up to this constant.

At last, an important point is that Eq. ( 4) is often intractable, even for two targets, and must then be computed numerically. We will see in the following that a closed-form can be obtained only for the Swerling 1 fluctuations.

For other fluctuations type, the numerical implementation implies the evaluation of multiple integrals over several parameters and the computational cost may be rapidly prohibitive in the multitarget case. Fortunately, target contributions can in many cases be separated so that the multitarget likelihood becomes equal to the product of monotarget likelihoods that can be computed in closed-form. This separation can be performed when targets do not interact in the likelihood computation, which can be translated mathematically by the following condition:

h H k,u Γ -1 h k,v ≈ 0, ∀(u, v), u = v, (5)
that permits to remove all cross terms in Eq.( 3). In practice, this hypothesis may arise for instance when Γ = I Nc and targets are far away from each other. Indeed, for each target the ambiguity vector h k,i has only significant values in a few number of cells around the target location and can be assumed equal to zero elsewhere. Therefore, the scalar product between ambiguity function h k,u and h k,v is approximately equal to zero for sufficiently distant targets. Note however that when Γ = I Nc , condition (5) cannot be verified as straightforwardly and should thus be carefully checked, even for distant targets; in particular the inner product induced by matrix Γ -1 may mix the target contributions, even when they are located far apart from each other.

December 2, 2014 DRAFT Finally under condition (5) the expression of likelihood p

(z k | X k , ρ k,1:N k , ϕ k,1:N k ) becomes p (z k | X k , ρ k,1:N k , ϕ k,1:N k ) ∝ exp - N k i=1 ρ 2 k,i h H k,i Γ -1 h k,i + N k i=1 2ρ k,i |h H k,i Γ -1 z k | cos (ϕ k,i -ξ k,i ) ∝ N k i=1 exp -ρ 2 k,i h H k,i Γ -1 h k,i + 2ρ k,i |h H k,i Γ -1 z k | cos (ϕ k,i -ξ k,i ) , (6)
where the i th term of the product, denoted by

Ξ z k ,x k,i (ρ k,i , ϕ k,i ) = exp -ρ 2 k,i h H k,i Γ -1 h k,i + 2ρ k,i |h H k,i Γ -1 z k | cos (ϕ k,i -ξ k,i ) , (7)
only depends on parameters ρ k,i and ϕ k,i . As variables ρ k,1:N k and ϕ k,1:N k are independent, the joint density (4) then simply becomes

p (z k |X k ) ∝ N k i=1 +∞ 0 2π 0 Ξ z k ,x k,i (ρ k,i , ϕ k,i )p(ϕ k,i )p(ρ k,i )dϕ k,i dρ k,i . (8)
Thus, everything happens as if each target is processed separately. This drastically alleviates the computational complexity of integral ( 4) and allows processing distant targets with parallel filters or by using the filter developed by Vo et al. [START_REF] Vo | Joint detection and estimation of multiple objects from image observations[END_REF] that made the hypothesis that likelihood can be factorized as in Eq. ( 8). Of course, when condition ( 5) is not verified, this simplification can be done only for separated targets, while targets that cannot be separated must be processed by the same filter.

In the monotarget case, integral (4) becomes

p (z k |X k ) ∝ +∞ 0 2π 0 p (z k |X k , ϕ k , ρ k )p(ϕ k )p(ρ k )dϕ k dρ k . ( 9
)
Davey et al. [1] have shown in this particular monotarget case that the marginalization can be done over the phase

ϕ k , providing p (z k | X k , ρ k ) ∝ 2π 0 p(z k |x k , ϕ k , ρ k )p(ϕ k )dϕ k , ∝ exp -ρ 2 k h H k Γ -1 h k I 0 2ρ k h H k Γ -1 z k . ( 10
)
Then, the likelihood is obtained by integrating [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF] over the density of parameter ρ k that depends on the fluctuation model considered.

B. Likelihood computation with squared modulus

In the previous section, we have presented the exact computation of the likelihood from complex measurements.

In this section, we now expose a different approach often considered in the literature, that consists in working only with the squared modulus of the complex data [START_REF] Rutten | Recursive track-before-detect with target amplitude fluctuations[END_REF], [START_REF] Davey | A comparison of detection performance for several track-before-detect algorithms[END_REF], [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF]. This approach is interesting in applications where only the squared modulus of the data is available but also because it permits to remove the phase dependency in a monotarget setting. This simplifies in some extent the computations, at the cost of losing the spatial coherence of the phase. Squared modulus were also considered in an application involving two targets with Swerling 1 amplitude fluctuations [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF]. In this specific application, the spatial coherence of the target amplitude was not considered, December 2, 2014 DRAFT thus simplifying the computation at the cost of some information loss. We will derive here the general multitarget likelihood in the squared modulus framework. It differs from expressions obtained in the literature since it does not make any approximation and thus properly takes into account the spatial coherence of the complex amplitude.

Moreover we show that the squared modulus approach does not permit in the multitarget setting to remove all phase dependencies. Thus, as with complex measurements, these phase variables must be taken into account, for instance by marginalization.

First, let us assume, as in the literature [1], [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF], [START_REF] Salmond | A particle filter for track-before-detect[END_REF], that the covariance matrix is given by Γ = 2σ

|z k | 2 = [|z 1 k | 2 , . . . , |z Nc k | 2 ] T .
As a consequence, samples from |z k | 2 are also independent conditionally to variables ρ k,1:N k and ϕ k,1:N k and we can write,

p(|z k | 2 | X k , ρ k,1:N k , ϕ k,1:N k ) = Nc l=1 p(|z l k | 2 | X k , ρ k,1:N k , ϕ k,1:N k ). (11)
The desired density

p(|z k | 2 | X k ) can then be obtained from p(|z k | 2 | X k , ρ k,1:N k , ϕ k,1:N k )
exactly in the same way as with complex measurements, by marginalizing over all variables ρ k,1:N k and ϕ k,1:N k . Remark that the hypothesis of independence is absolutely necessary here to establish Eq.( 11). The condition Γ = 2σ 2 I Nc can be generalized to diagonal covariance matrices, but the case where Γ is not a diagonal matrix is much more complicated even for two coupled variables: in that case, squared modulus samples are correlated, thus leading to distributions with no closed-form, for instance multivariate Rayleigh distribution in the Swerling 1 case [START_REF] Mallik | On multivariate Rayleigh and exponential distributions[END_REF]. Note also that in practice, this hypothesis is verified with classic matched filtering in presence of white Gaussian noise and an appropriate sampling rate, but it may not be verified anymore when modifying the reception processing, for instance by applying classic weighting windows such as Hamming, Bartlett, Hann, etc. [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF]. Before going further into the computation, we would like to highlight here an interesting property that arises when considering squared modulus of complex data, and that has never been discussed to our knowledge in the literature: although N k targets are present, providing N k different and independent random phases ϕ k,1:N k , it is possible to show, by changing the set of parameters,

that density p(|z l k | 2 | X k , ρ k,1:N k , ϕ k,1:N k ) effectively depends only on N k -1 phase variables. Indeed the variable |z l k | 2 can be defined up to an arbitrary phase ϕ since |z l k | 2 = |e ϕ z l k | 2
, and we can write for instance

z l k 2 = ρ k,1 h l k,1 + N k i=2 ρ k,i e jϕ k,i h l k,i + n l k 2 ( 12
)
where all n l k = n l k e -jϕ k,1 are still independent circular symmetric complex Gaussian noise samples, and phases ϕ k,i = ϕ k,i -ϕ k,1 are still uniform variables distributed over the interval [0, 2π). Thus, |z l k | 2 only depends on N k -1 phase variables. Therefore, taking the squared modulus of the complex signal leads to drop out the dependence of one and only one phase. As a consequence in a monotarget setting the density of |z l k | 2 does not depend any longer on December 2, 2014 DRAFT the phase ϕ k but only on the modulus; this is one of the main reasons to use such a technique in a TBD monotarget algorithm. On the contrary, in the multitarget setting, taking the squared modulus does not remove all dependencies on the phases! This dependency remains present through coherent summations of the target contributions in each cell. Discarding it may lead to losing all the information provided by the spatial coherence of the phase variables.

Conditionally to variables X k , ρ k,1:N k and ϕ k,2:N k , each sample

|z l k | 2
σ 2 follows a non central chi-square distribution with two degrees of freedom; indeed it corresponds to the sum of the square of two non-centered Gaussian variables.

The density p(|z

l k | 2 | X k , ρ k,1:N k , ϕ k,2:N k )
is thus given by:

p(|z l k | 2 | X k , ρ k,1:N k , ϕ k,2:N k ) = 1 2σ 2 exp    - |z l k | 2 2σ 2 - γ l ϕ k,2:N k , ρ k,1:N k 2    I 0     γ l ϕ k,2:N k , ρ k,1:N k |z l k | 2 σ 2     , ( 13
)
where I 0 is the modified Bessel function of the first kind and γ l ϕ k,2:N k , ρ k,1:N k is the non centrality parameter equal to

γ l ϕ k,2:N k , ρ k,1:N k = ρ k,i h l k,i + N k i=2 ρ k,1 e jϕ k,1 h l k,i 2 σ 2 . ( 14
)
At this step, mono and multitarget cases are different, and we will consider them separately in the following.

Finally, note that, as with complex measurements, the likelihood can be computed up to a constant. Therefore terms

1 2σ 2 exp - |z l k | 2 2σ 2
will be discarded in the rest of the paper.

1) The monotarget case: In a monotarget setting, the non-centrality parameter in each cell becomes

γ l (ρ k ) = ρ 2 k |h l k | 2 σ 2 (15)
and does not depend on ϕ k . The joint likelihood can then be obtained by marginalizing Eq.( 11) over the parameter

ρ k : p(|z k | 2 | X k ) = ∞ 0 Nc l=1 p(|z l k | 2 | X k , ρ k )p(ρ k )dρ k , ( 16
)
where p(ρ k ) is the density for the parameter ρ k . As for complex measurements, this marginalization allows preserving the spatial coherence of the parameter ρ k . Since integral ( 16) is, to our knowledge, intractable for Swerling fluctuations models of type 1 and 3 since it consists in integrating N c Bessel functions, it must be in that case approximated numerically. An heuristic solution was proposed by Rutten et al. [START_REF] Rutten | Recursive track-before-detect with target amplitude fluctuations[END_REF] that consists in assuming

that samples |z 1 k | 2 , . . . , |z Nc k | 2 are independent.
Recall that this is not true in general because of the spatial coherence of random variable ρ k . Under this assumption, it comes that

p(|z k | 2 | X k ) = Nc l=1 p(|z l k | 2 | X k ), (17)
where the likelihood in each cell is obtained by marginalization over ρ k :

p(|z l k | 2 | X k ) = ∞ 0 p(|z l k | 2 | ρ k , X k )p(ρ k )dρ k . ( 18
)
December 2, 2014 DRAFT

Clearly the spatial coherence of ρ k is lost as it is integrated independently in each cell and not over the whole likelihood. On the other hand, the calculation of integral ( 18) can be done analytically for Swerling fluctuation models of type 1 and 3, leading to simple closed-forms.

2) The multitarget case: As previously discussed, in the multitarget case the parameter γ l ϕ k,2:N k , ρ k,1:N k still depends on the phase variables ϕ k,2:N k . The likelihood must thus be obtained by marginalization over modulus

ρ k,1:N k and phases ϕ k,2:N k : p(|z k | 2 | X k ) = • • • R N k ≥0 ×[0,2π) N k -1 Nc l=1 p(|z l k | 2 | X k , ρ k,1:N k , ϕ k,2:N k )p(ρ k,1:N k )p(ϕ k,2:N k )dρ k,1:N k dϕ k,2:N k . ( 19
)
As in the monotarget case, this expression is to our knowledge intractable. An heuristic solution consists then again in marginalizing independently each sample from ϕ k,2:N k and ρ k,1:N k as in [START_REF] Vo | Joint detection and estimation of multiple objects from image observations[END_REF], providing

p(|z l k | 2 | X k ) = • • • R N k ≥0 ×[0,2π) N k -1 p(|z l k | 2 | X k , ρ k,1:N k , ϕ k,2:N k )p(ρ k,1:N k )p(ϕ k,2:N k )dρ k,1:N k ϕ k,2:N k . ( 20
)
Note, however, that contrary to the monotarget case there is in general no closed-form for the integral [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF], so that numerical integration must still be performed.

Finally, as with complex measurements, target contributions can often be separated so that the multitarget likelihood then resorts to a product of monotarget likelihoods. This separation is obtained under the condition h l k,i h l k,j ≈ 0, ∀i, j, i = j that permits to eliminate all cross terms in Eq.( 14). It can be obtained for instance when targets are far apart enough so that their respective contributions do not overlap spatially.

IV. LIKELIHOOD COMPUTATION FOR SWERLING MODELS

In this section, we will derive the measurement likelihood with three different Swerling models: Swerling 0, Swerling 1 and Swerling 3. For each model, we will consider first the case of complex measurements and second the case of squared modulus measurements. Whenever closed-forms are not obtainable, we will propose approximations that permit to compute the likelihood at a lower cost.

A. Complex measurements 1) Swerling 0 case: The modulus ρ k,i of each target is assumed constant and equal to ρ i . Although these constants are unknown, they can be estimated over time and their estimates injected into the likelihood. For instance, in particle filter, amplitude parameters ρ 1:N k can be inserted in the state vector, with artificial Markovian dynamics, and sampled as the other state parameters [START_REF] Andrieu | Particle methods for change detection, system identification, and control[END_REF]. Then integral (4) that corresponds to the complex measurement likelihood must just be computed over parameters ϕ k,1:N k . In the general multi-target case, this integral is, according to our knowledge, intractable and must be approximated except for the particular single target case. A first solution consists in calculating numerically the integral over the domain [0, 2π) N k but this may become rapidly computationally demanding. Thus, we propose to replace the intractable likelihood by its Laplace approximation December 2, 2014 DRAFT that has been already successfully used in particle filter application [START_REF] Musso | Introducing the Laplace approximation in particle filtering[END_REF]. Let 1 , . . . , e jϕ k,N k ] T . Equation ( 4) can be rewritten as follows:

H k = [ρ 1 h k,1 , . . . , ρ N k h k,N k ] and let Ψ k = [e jϕ k,
p (z k | X k ) ∝ [0,2π) N k exp {Υ (ϕ k,1:N k )} dϕ k,1:N k . ( 21
)
where

Υ (ϕ k,1:N k ) = - N k i=1 ρ 2 i h H k,i Γ -1 h k,i + N k i=1 2ρ i |h H k,i Γ -1 z k | cos (ϕ k,i -ξ k,i )- N k i=1 N k l=i+1 2ρ i ρ l h H k,i Γ -1 h k,l cos (ϕ k,i -ϕ k,l -φ k,il ) .
The integral ( 21) can be approximated using the Laplace's method [START_REF] Musso | Introducing the Laplace approximation in particle filtering[END_REF] and is then given by

p SW 0 (z k | X k ) ≈ exp {Υ ( ϕ k,1:N k )} (2π) Nr 2 |det (-∇ 2 Υ ( ϕ k,1:N k ))| 1 2 , ( 22
)
where ϕ k,1:N k is the maximum of function Υ (.) and ∇ 2 Υ ( ϕ k,1:N k ) is the Jacobian matrix calculated at the maximum. This maximum cannot be obtained analytically even for two targets and an optimization method such as a gradient descent must be used. However, this optimization step can be avoided by using the classic least square estimator

Ψ k = H H k Γ -1 H k -1 H H k Γ -1 z k ( 23
)
of the complex amplitude and approximating each component by

ϕ k,i = arg Ψ k,i . ( 24
)
In practice, this estimator is in most of the situations close to the actual maximum. However, in some situations, for instance when components h k,1:N k are almost colinear, the difference can be greater. In that latter case, an optimization must be performed or the filter will loss in sensitivity. A compromise must then be done between the quality of the estimate and the computational time required to reach it.

2) Swerling 1 case: Each modulus ρ k,i follows a Rayleigh distribution:

p SW 1 (ρ k,i ) = ρ k σ 2 ρi exp - ρ 2 k,i 2σ 2 ρi ( 25
)
where σ ρi is the parameter of the Rayleigh distribution such that E[ρ 2 k,i ] = 2σ 2 ρi , while the phases ϕ k,1:N k are uniformly distributed over [0, 2π). As for the Swerling 0 case, this parameter is unknown but it can be added to state vector as proposed in [START_REF] Andrieu | Particle methods for change detection, system identification, and control[END_REF]. Although the integral (4) with respect to the Swerling 1 densities for parameters ρ k,1:N k seems to be intractable directly, in practice the density p(z k | X k ) can be obtained using other probabilistic considerations. Indeed, in the Swerling 1 model, since ρ k,i follows a Rayleigh distribution of parameter σ ρ,i and ϕ k,i is uniformly distributed over [0, 2π), each variable ρ k,i e jϕ k,i in the measurement equation ( 1) is a zero-mean circular symmetric complex Gaussian variable with variance 2σ 2 ρ,i . Therefore z k , which is then the sum of independent December 2, 2014 DRAFT Gaussian vectors with zero-mean, is a complex Gaussian vector with zero-mean and covariance matrix Σ N k given by

Σ N k = Γ + N k i=1 2σ 2 ρ,i h k,i h H k,i . (26)
Clearly, this matrix is definite positive, so that the multi-target likelihood is finally given in closed form by:

p SW 1 (z k | X k ) ∝ 1 det (Σ N k ) exp -z H k Σ -1 N k z k . ( 27
)
In practice, the computation of the likelihood requires the evaluation of det (Σ N k ) and Σ -1 N k that can be computationally demanding since matrix Σ N k is a square matrix of size equal to the length of considered vector h k,i . Fortunately, these quantities can be easily computed by using classic linear algebra formulas. Indeed, the matrix Σ N k can be written

Σ N k = Γ + UVU H , ( 28
)
where

U = [h k,1 , • • • , h k,N k ] is a matrix with N k columns and V = diag 2σ 2 ρ,1 , • • • , 2σ 2 ρ,N k .
Then using a classic matrix inversion lemma (see [START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF], p. 117), it comes

Σ -1 N k = Γ -1 -Γ -1 U V -1 + U H Γ -1 U -1 U H Γ -1 . ( 29
)
The inverse of matrix Γ can be pre-computed, while V is a diagonal matrix and matrix

V -1 + U H Γ -1 U is an
N k -by-N k matrix of much smaller size than Σ N k as long as the number of targets N k remains small compared to the number of considered cells. In that case its inversion implies a drastically reduced cost compared to the inversion of Σ N k . Furthermore, the computational cost of the determinant can also be reduced using the matrix determinant lemma (see [START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF], p. 117)

det (Σ N k ) = det V -1 + U H Γ -1 U det (V) det (Γ) . ( 30
)
Note that no hypothesis was made here about the closeness of the targets and therefore this closed-form expression is valid both for distant and close targets. Finally, for the particular monotarget case, the likelihood simply becomes

pSW 1 (z k | X k ) ∝ 1 1 + 2σ 2 ρ h H k Γ -1 h k exp 2σ 2 ρ |h H k Γ -1 z k | 2 1 + 2σ 2 ρ h H k Γ -1 h k . ( 31
)
3) Swerling 3 case: Each squared modulus ρ 2 k,i follows a chi-square distribution with four degrees of freedom so that the corresponding density for the modulus ρ k,i is provided by:

p SW 3 (ρ k,i ) = 8ρ 3 k,i ν 2 ρi exp - 2ρ 2 k,i ν ρi , ( 32
)
where the parameter ν ρi is such that E[ρ 2 k,i ] = ν ρi . As for the Swerling 0 and Swerling 1 case, the parameter ν ρi is unknown but can be added to the state vector [START_REF] Andrieu | Particle methods for change detection, system identification, and control[END_REF]. According to our knowledge, no closed form can be obtained for Eq (4) in the Swerling 3 case and a numerical approximation must be done, implying the numerical computation of N k integrals over modulus ρ k,1:N k and N k integrals over phases ϕ k,1:N k . However, it is possible to December 2, 2014 DRAFT avoid the numerical integration over the parameters ρ k,1:N k by approximating the chi-square distribution by a Rice distribution; note indeed that the Swerling 3 model can be viewed as an approximation of a Rice distribution [START_REF] Richards | Relationship between the gamma, Erlang, chi-square, and Swerling 3/4 probability density functions[END_REF].

Using a Rice distribution instead of the Swerling 3 model, the density of the modulus ρ k,i becomes

p Rice (ρ k,i ) = 2ρ k,i (1 + a 2 ) ν ρi exp -a 2 -ρ 2 k,i (1 + a 2 ) ν ρi I 0 2a ρ 2 k,i (1 + a 2 ) ν ρi , ( 33
)
where a is the ratio between the dominant scatterer and the weaker ones. By choosing a = 1 + √ 2, it can be easily checked that densities of the squared modulus ρ 2 k,i under Swerling 3 and Rice models provide the same means and variances [START_REF] Richards | Relationship between the gamma, Erlang, chi-square, and Swerling 3/4 probability density functions[END_REF]. Now consider the complex amplitude ρ k,i e jϕ k,i where ρ k,i is distributed according to the Rice distribution (33). Recall first that this Rice distribution is the distribution of the modulus of a complex Gaussian variable with mean µ SW 3,i = a νρ i (1+a 2 ) and variance 2σ 2 SW 3,i = νρ i (1+a 2 ) . Then we can replace each variable ρ k,i e jϕ k,i in (1) by a variable ξ k,i e jψ k,i where the variables ξ k,i and ψ k,i are respectively Gaussian and uniform, and such that ξ k,i e jψ k,i follows the same distribution as ρ k,i e jϕ k,i . Then, conditionally to phases ψ k,1:N k , the observation z k is a complex Gaussian vector with mean

µ k,SW 3 = N k i=1 µ SW 3,i e jψ k,i h k,i
and covariance matrix

Φ N k = Γ + N k i=1 2σ 2 SW 3,i h k,i h H k,i .
The density is then given by,

p Rice (z k | X k , ψ k,1:N k ) ∝ 1 det (Φ N k ) exp -z k -µ k,SW 3 H Φ -1 N k z k -µ k,SW 3 (34)
Clearly, the computational cost of Φ -1 N k and det (Φ N k ) can be reduced as in the Swerling 1 case. Then, it just remains to marginalize (34) over the phases ψ k,1:N k . This marginalization cannot be computed analytically and must then be calculated numerically, except in the monotarget case.

In the particular monotarget case, a closed-form can be obtained both for the chi-square distribution and the Rice distribution. For the chi-square distribution, the expression in Eq.( 10) must be integrated over density (32). The following result (see [START_REF] Gradshteyn | Table of integrals, series, and products[END_REF], p. 1097 Eq. 6.663)

+∞ 0 x 3 exp -αx 2 I 0 (βx) dx = 2 α 2 1 + β 4α exp β 2 4α , ( 35
)
where

α ∈ R * ≥0 and β ∈ R, is used with α = 2 νρ + h H k Γ -1 h k and β = 2 h H k Γ -1 z k .
Then, the likelihood for the chi-square Swerling 3 model in the monotarget case is given by

p SW 3 (z k | X k ) ∝ 4 2 + ν ρ h H k Γ -1 h k 2 1 + ν ρ h H k Γ -1 z k 2 + ν ρ h H k Γ -1 h k exp ν ρ h H k Γ -1 z k 2 2 + ν ρ h H k Γ -1 h k . ( 36
)
For the Rice distribution, it is possible to integrate Eq.(34) over the phase ψ, a computation similar to the one providing Eq.( 10), Then, the likelihood for the Rice Swerling 3 model in the monotarget setting is given by

p Rice (z k | X k ) ∝ 1 + a 2 exp -a 2 1 + a 2 + ν ρ h H k Γ -1 h k exp ν ρ h H k Γ -1 z k 2 + a 2 1 + a 2 1 + a 2 + ν ρ h H k Γ -1 h k I 0 2a h H k Γ -1 z k (1 + a 2 ) ν ρ 1 + a 2 + ν ρ h H k Γ -1 h k B.

Squared modulus measurements

As we have seen, the likelihood computation with the squared modulus can be done in two ways, either by taking into account the spatial coherence of the phases and modulus with Eq.( 19) or by marginalizing independently in each cell with Eq. [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF]. As this two cases are different we treat them separately in the following.

1) The coherent case: In the coherent case, the likelihood is obtained according to Eq.( 19) by replacing the generic density p (ρ k,i ) by the density of the fluctuation considered. However, according to our knowledge, it cannot be done analytically for the Swerling models and the integral must be approximated numerically. Moreover, Note that it can be really intensive in terms of computational resources especially when the number of target is large since the size of the integration domain increases exponentially with the number of targets. For this reason,

we propose an heuristic solution that consists in replacing the parameter γ l (ρ k,1:N k , ϕ k,2:N k ) by its expectation

E[γ l (ρ k,1:N k , ϕ k,2:N k )] = N k i=1 E ρ 2 i |h l k,i | 2 σ 2 , ( 38
)
where E ρ 2 i just depends on the parameter of the fluctuations density. Thus, integrals [START_REF] Mallik | On multivariate Rayleigh and exponential distributions[END_REF] are just the product of the densities in Eq.( 11) for all the cells. This is a strong approximation for the likelihood, but as you will see in V, it gives interesting performance and it is really faster than the numerical integration which is costly in terms of computational resources. In the monotarget case, the likelihood is given by Eq.( 16) that requires the integration only over parameter ρ k and therefore the numerical approximation can be done in a reasonable time.

2) The non coherent case: The non coherent case consists in calculating Eq. [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF]. In practice for the Swerling 0 case, it has no interest and Eq.( 19) must be used as it takes into account the spatial coherence of variables ϕ k,2:N k .

Nevertheless, for the Swerling 1 and 3 cases, probabilistic considerations can be used to calculate Eq. [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF]. Indeed, in the Swerling 1 case Boers et al. [START_REF] Boers | Multitarget particle filter track-before-detect application[END_REF] noticed that sample |z l k | 2 follows an exponential distribution with parameter

λ l k = 1 2σ 2 + N k i=1 2σ 2 ρ i |h l k | 2 , so that p(|z l k | 2 | X k ) = 1 λ l k exp - |z l k | 2 λ l k . ( 39
)
For the Swerling 3 case, the integration over parameters ρ k,1:N k can be avoided with the Rice fluctuations. Indeed, by replacing each variable ρ k,i e jϕ k,i by a variable ξ k,i e jψ k,i , each sample

|z l k | 2 σ 2 + N k i=1 σ 2 SW 3,i h l k 2
conditionally on variables ψ k,1:N k follows a non central chi-square distribution with two degrees of freedom; with non-centrality parameter

γ l Rice ψ k,2:N k = µ SW 3,1 h l k,1 + N k i=2 µ SW 3,i e jψ k,i h l k,i 2 σ 2 + N k i=1 σ 2 SW 3,i h l k 2 , ( 40
)
December 2, 2014 DRAFT that does not depend on parameters ρ k,1:N k anymore. The density of |z l k | 2 conditionally on ψ k,2:N k is given by Eq. ( 13) where σ 2 is substituted by

σ 2 + N k i=1 σ 2 SW 3,i h l k 2 and γ l ϕ k,2:N k , ρ k,1:N k by γ l Rice ψ k,2:N k . Finally the likelihood |z l k | 2
is obtained only by integrating over variables ψ k,2:N k . In the monotarget case integral (18) can be computed analytically both for the Rice distribution and the chi-square. For the Rice distribution, no integration over phase ψ k,1 is required and the likelihood is given by

pSW 3,Rice z l k 2 | X k ∝ 2σ 2 (1 + a 2 ) exp -a 2 2σ 2 (1 + a 2 ) + νρ h l k 2 exp   νρ |h l k | 2 |z l k | 2 2σ 2 + 2σ 2 a 2 1 + a 2 2σ 2 (1 + a 2 ) + νρ h l k 2   I0 2a h l k z l k (1 + a 2 ) νρ 1 + a 2 + νρ h l k 2 . ( 41
)
For the chi-square distribution, result ( 35) is used with α =

νρ|h l k | 2 +4σ 2 2νρσ 2 and β = |h l k ||z l k | σ 2
. Then, integral [START_REF] Vo | Joint detection and estimation of multiple objects from image observations[END_REF] becomes

p SW 3,χ 2 z l k 2 | X k ∝ 4σ 2 2 4σ 2 + ν ρ h l k 2 2 1 + z l k 2 2σ 2 ν ρ h l k 2 4σ 2 + ν ρ h l k 2 exp z l k 2 2σ 2 ν ρ h l k 2 4σ 2 + ν ρ h l k 2 (42)

C. Summary

In this section, we have given several way to compute the likelihood in a Track-Before-Detect context for complex amplitudes fluctuation of type Swerling 0, 1 and 3. For the computation of the Likelihood with the complex measurement, we have shown that a closed-form can be obtained for all the Swerling fluctuations considered in the monotarget case. In the multitarget case, a closed-form can be obtained only in the Swerling 1 case, in the other cases a numerical integration must be performed, however we propose several methods in order to alleviate the time calculation. For the Likelihood with the squared modulus of the complex measurement, we have derived the right expression in order to keep the spatial coherence information of complex amplitude parameters and we have shown that only the dependency of one phase can be removed, however this leads to an intractable integral for all the Swerling models. Then approximation must be performed and we give the principal way to do so. In table IV, we give a sum-up of the different techniques to calculate the likelihood with the existing methods on those proposed in this article.

V. SIMULATION AND RESULTS

In this section, we first study the performance in detection and estimation of a single target particle filter that considers either complex measurements or squared modulus. We show the improvement of using complex measurements both in detection and in estimation only for the Swerling 1 and 3 model as Davey et al. have already shown the benefits of doing so in the Swerling 0 case. Then, we study the behaviour of a simple multitarget particle filter for two close targets. Performance are evaluated in terms of estimation of the two target states and track loss for fluctuations of type Swerling 0, 1 and 3.

A. Target Model

We assume a discrete time model, with a fixed time step T, and we define the state vector of the i th target as

x k,i = [x k,i , ẋk,i , y k,i , ẏk,i ] T
, where (x k,i , y k,i ) and ( ẋk,i , ẏk,i ) are, respectively, the location and the velocity of December 2, 2014 DRAFT the target in Cartesian coordinates. If the i th target appears at step k, its position (x k,i , y k,i ) is uniformly drawn in the pavement defined in polar coordinate P = [r min , r max ] × [θ min , θ max ] where r min , r max , θ min and θ max are respectively the minimum and maximum target ranges and bearings. The velocity ( ẋk,i , ẏk,i ) is uniformly drawn

in the area C = {( ẋk , ẏk ) | v min ẋ2 k + ẏ2 k v max }
where v min and v max are respectively the minimum and maximum target velocity. On the contrary, if the i th target was already present at previous step k -1, we assume that the target state evolves according to the following linear equation:

x k,i = Fx k-1,i + v k,i , (43)
where v k,i is a white Gaussian noise with covariance matrix

Q =   Q s 0 0 Q s   , with Q s = q s   T 3 /3 T 2 /2 T 2 /2 T   ,
and F is the transition matrix defined by:

F =   F s 0 0 F s   with F s =   1 T 0 1   .
Concerning the modulus, we consider fluctuations of type Swerling 0, 1 and 3 with parameters ρ i , 2σ 2 ρi and ν ρi respectively.

B. Ambiguity function

In Eq.( 1), the ambiguity function h k (x k ) is possibly multidimensional depending on the radar application under consideration. We consider here a simple scenario with only range and bearing measurements.For the range, the transmitted pulse is assumed to be a linear frequency modulated signal ("chirp") with band B and duration T e , whose range ambiguity function is given by [START_REF] Levanon | Radar signals[END_REF]:

h l r (x k ) = sin πBτ l 1 -|τ l | Te πBτ l for |τ l | ≤ T e ,
with τ l = 2 (r k -r l ) /c, c the celerity of the electromagnetic wave,

r k = x 2 k + y 2 k and r l = r min +(l+ 1 2 )∆ r , l ∈ [0, N r -1] with ∆ r = c
2B the range resolution and N r = rmax-rmin ∆r . For the bearing, we consider at the reception a linear phased array of N a antennas spaced by λ 2 (where λ is the wavelength of the carrier frequency). Then the ambiguity in bearing is given by [START_REF] Van Trees | Optimum Array Processing, Part IV of Detection, Estimation, and Modulation Theory[END_REF],

h m θ (x k ) = sin NaΦ m 2 N a sin Φ m 2 , with Φ m = 2πda λ [cos(θ k ) -cos (θ m )], θ k = arctan( y k x k ) the target bearing and θ m = θ min + (m + 1 2 )∆ θ , m ∈ [0, N θ -1] with ∆ θ = 1.772 λ
Na the half-power beamwidth [START_REF] Van Trees | Optimum Array Processing, Part IV of Detection, Estimation, and Modulation Theory[END_REF] and

N θ = θmax-θmin ∆ θ .
The overall range-bearing ambiguity function in cell (l, m) is given by the product

h lm (x k ) = h l r (x k ) × h m θ (x k ). Finally h k (x k ) is a vector of size N c = N r × N θ defined by h k (x k ) = h 11 (x k ), h 12 (x k ), • • • , h 12 (x k ), • • • , h NrN θ (x k ) .
December 2, 2014 DRAFT (where P f a is a given false alarm probability). Concerning the velocity, we simply choose a uniform prior over the area C defined in V-A. Finally, parameters σ 2 0,p and ν 0,p are drawn uniformly over the interval corresponding to a target SNR between SNR min and SNR max .

The approximated posterior density is a mixture with two components, one for the birth particles and one for the continuing particles. Therefore unnormalized weights are calculated separately. For the birth particles it is given by

wp k,b ∝ N P f a N × 1 N b p (z k | x p k ) , (46)
where N P f a is the number of cells exceeding the threshold ν P f a . For the continuing particles it is given by

wp k,c ∝ 1 N c p (z k | x p k ) . ( 47
)
Finally, The probability of each component is given by

p c = Mc Mc + Mb , p b = 1 -p c , (48)
where

Mc = (1 -P d ) P e k-1 Nco p=1 wc k,p , Mb = P b (1 -P e k-1 ) N b q=1 wb k,q , (49)
with P e k-1 the estimated probability of target existence at step k-1. The estimated probability P e k of target existence is obtained by

P e k = Mc + Mb Mc + Mb + P d P e k-1 + (1 -P b )(1 -P e k-1 ) . ( 50
)
The weights are normalized separately for the birth and continuing particles. Finally, N c particles are resampled from the N co + N b particle cloud.

3) Detection scheme: The proposed particle filter allows to estimate the probability of existence P e k but does not take any decision on the presence or absence of the target. It appears that a common detection scheme consists in simply thresholding the probability of existence. However, this way of proceed can lead to miss detections. Indeed for an outlier measurement z k , the estimated probability of existence P e k may become pretty small in one iteration, even though particles are properly located around the target; in that case the target will not be detected at that iteration. To avoid this issue, we use an adaptive threshold that depends on the detection status of the target at previous step.

Let us thus call d k the decision at time step k that can take value 0 (no detection) or 1 (detection), the following detection scheme is used: For the simulation of the radar measurements, the parameters used are: r min = 100 km,r max = 120 km θ min = -10 • , θ max = +10 • , N r = 40, N θ = 14, σ 2 = 0.5, B = 150 kHz, T e = 6.67 × 10 -5 s, N a = 70, λ = 3 cm, c = 3 × 10 8 m.s -1 . Note that a small radar window is chosen here to avoid using an important number of particles and thus limit the computational cost.

d k =    1 if P e k > T h (d k-1 ), 0 otherwise.
Three filters are used to detect and estimate the hidden target state x k , based on different assumptions for the likelihood computation:

1) The first filter, denoted by CSM PF (Coherent Squared Modulus Particle Filter), considers squared modulus to compute the likelihood and takes into account the spatial coherence of the amplitude parameter ρ k : it corresponds to Eq. ( 16).

2) The second filter, denoted by NCSM PF (Non Coherent Squared Modulus Particle Filter), considers squared modulus but does not take into account the spatial coherence of the amplitude parameter ρ k : it corresponds to Eq. [START_REF] Vo | Joint detection and estimation of multiple objects from image observations[END_REF].

3) The third filter, denoted by CM PF (Complex Measurement Particle Filter), considers complex measurements and spatial coherence: it corresponds to Eq. [START_REF] Kreucher | Multitarget Tracking using the Joint Multitarget Probability Density[END_REF].

N M C = 1000 Monte-Carlo simulation were performed for performance measurement. Finally, the RMSE in position and velocity is measured when d k = 1 and the estimated target location is in the detection area of two range-bearing cells from the actual target position. a) Detection performance: In figures 2 and 3, we present the average of variable d k is measured at each step for the Swerling 1 and 3 models respectively. In both case, filters that use the complex measurement outperform those that use squared modulus. Furthermore, the difference between SMCPF and SMNCPF is quite small, therefore it seems that taking into account the spatial coherence of the phase is more important than taking into account the modulus information. Moreover, SMCPF requires numerical approximation that leads to increase the computational time for a very small gain in detection. b) Estimation performance: In figures 4 and 5, we present the result in terms of RMSE in position and velocity for the Swerling 1 and 3 models respectively. As for all the detection results, particle filters that used the complex measurement outperform filters that work on squared modulus. Moreover, note that the RMSE in position seems to be better at the beginning which is not expected since the tracking alogrithm should improve the RMSE. However, this can be explained by the fact that the RMSE is calculated only over the iteration where the target has been detected (i.e. d k = 1) and at the beginning only a few simulations have detected the target (recall the SNR of 5dB used here). These detections correspond to favorable cases After a few iteration, target has been detected in most of December 2, 2014 DRAFT the case and thus particularly cases where is located at the edge of the cell and thus it has lost of few dB, which leads to increase the RMSE. Besides with the filter used here, some birth particles are spread everywhere in the observation window, and used for the situation as well, thus leading to an increase in the RMSE. This phenomena could be managed by clustering the particle cloud before the estimate computation.

D. Multi-Target Simulation and Results

1) Multi-Target Scenario:

We now consider a scenario with two targets present during all the experiment. Both targets follow a uniform rectilinear trajectory. Target states x k,1 and x k,2 are uniformly initialized over P × C such that:

• the two velocity vectors ( ẋk,1 , ẏk,1 ), ( ẋk,2 , ẏk,2 ) form an angle of π 4 .

• the minimum distance between targets is reached at time step k c = 35 and is set to d min = 150 m, i.e. the minimum distance is smaller than the range resolution.

Targets are set to 10dB and we consider fluctuations of type Swerling 0, 1 and 3.

2) Multi-Target Particle Filter: For the simulation, we consider here the particle filter proposed by Kreucher et al. [START_REF] Kreucher | Multitarget Tracking using the Joint Multitarget Probability Density[END_REF]. We assume that the number of targets is known since the objective here is to measure the effect of the likelihood computation on the particle filter for two close targets. Therefore, the particle state is defined as

X p k = [x p k,1 , x p k,2
] T , where x p k,1 and x p k,2 are the single state vectors of the first and second targets respectively of particle p, p ∈ {1, . . . , N p }. In the following, we detail the instrumental density used in the particle filter. At step k = 0, each particle target state x p 0,i is initialized from the actual target state according to the following procedure:

• for the position, a Gaussian noise with variance σ 2 r is added to the actual target range r 0,i = x 2 0,i + y 2 0,i and a Gaussian noise with variance σ 2 θ is added to the true target bearing θ 0,i = arctan For the particle propagation, we consider two cases:

• Either for each particle, state X p k verifies [START_REF] Blackman | Multiple-target tracking with radar applications[END_REF]. Then, the likelihood for each target state x p k,i can be computed separately and we propose to use the Independent Partition instrumental density (IP) [START_REF] Kreucher | Multitarget Tracking using the Joint Multitarget Probability Density[END_REF], i.e. sample the state of the particles according to the distributions defined by the likelihood of each target.

• Or hypothesis (5) is not verified for all the particles and (IP) cannot be used any longer. In that latter case, we just propagate particles according to the prior distribution Eq. (43).

3) Calculation of probability of track loss:

The probability of track loss is evaluated from N M C Monte-Carlo simulation with the following procedure: at each time step k and for each target, we compute the binary loss variable

l k,i =          1 if   rk,i -r k,i θk,i -θ k,i   P   rk,i -r k,i θk,i -θ k,i   > α, 0 otherwise, (52)
where

rk,i = x2 k,i + ŷ2 k,i , θk,i = arctan ŷk,i xk,i , P =   1 ∆r 0 0 1 ∆ θ 
 and α = 5.99 is the value of the quantile function of the chi-square distribution with two degrees of freedom evaluated at 0.95. In other words, at each iteration, we check if the position estimator for each target is located within the 0.95% confidence ellipse around the true target position. Finally, a track is declared to be lost if at least one of the variables l k,i equals 1 during at least five consecutive iterations. We define by f m the loss variable for the m th Monte-Carlo run that takes value 1 if the filter failed to track the two targets during all the experiment and 0, otherwise. Then, the probability of track loss is given by Ploss

= 1 N M C NMC m=1 f m .

4) Calculation of the Root Mean Square Error (RMSE):

The mean RMSE of the two targets is estimated from N M C Monte-Carlo runs with the following procedure: at each iteration, we obtain an estimator of the target state for each target provided by

xk,i = 1 N p Np p=1 x p k,i , i ∈ {1, 2},
and we associate each estimator to a target such that the sum of the Euclidean distances between the estimates and the actual state is minimum. Finally, the RMSE is computed at each iteration k for simulations where both targets have not been declared lost (i.e. l k,1 = 0 and l k,2 = 0) by taking the mean RMSE of the two targets over these simulations.

5) Simulations:

The particle filter is performed with the following parameters: T = 1 s, q s = 10 -3 , σ 2 r = 3.6 × 10 -3 , σ 2 θ = 1.022 × 10 -4 , σ 2 v = 0.01 and σ 2 n = 0.1. Parameters for the simulation of the radar measurements are the same as for the monotarget simulation, except for the radar window for which we take r min = 100 km, r max = 150 km, θ min = -20 • and θ min = +20 Then, as for the monotarget case, performance is evaluated for the three different filters CSM, NCSM and CM already defined, a fourth filter is also used and denoted by ESM (Expectation Squared Modulus) and corresponds to the case where the expectation of the non-centrality parameter is taken to compute the likelihood. Note that for the Swerling 0 case there is no interest of using NCSM since CSM requires integration only over N k -1 phases, therefore we replace the SMNC filter by the LC (for Laplace Coherent), where the likelihood is calculated via its Laplace approximation (see IV-A1).

When the particle states x p k,1 and x p k,2 are well separated, the likelihoods are calculated in closed-form according to the corresponding monotarget likelihood expression. When particle states are too close to each other to be assumed disjoint, the likelihoods are computed according to the multitarget likelihood expressions. When this computation requires a numerical integration, this integration is done over 10 points for each parameter. This small number of integration points is explained by the overall computational cost induced when several parameter dimensions are involved.

a) Estimation performance: the performance in terms of RMSE in position and velocity is presented in figures 6, 7 and 8 for the Swerling 0, 1 and 3 models respectively. First we observe that in all cases, CM provides the best performance. Then, the difference between the CSM and NCSM is quite small so that it does not seem relevant to take into account the spatial coherence of parameters ρ k,1:N k and ϕ k,1:N k with squared modulus. An other important point is to compare the computational time with respect to performance. Thus, in Swerling 0 the CL filter is approximatively six times faster than CM filter with almost the same performance. Likewise, in Swerling 1 and Swerling 3, the NCSM filters are approximatively 60 times faster than CSM. Finally, note that the RMSE in velocity increases when targets are close. This can be explained by the fact that the likelihood does not depend directly on the velocity. For all the Swerling models, the track-loss is minimum for the complex measurement, but filter CSM and NCSM are relatively close to it. The poorest performance is obtained with the ESM filter where the likelihood is computed with a rough approximation but has the advantage to be much more faster than CSM and NCSM.

VI. CONCLUSION

In this paper, we have investigated different methods for computing the likelihood in a radar Track-Before-Detect context. In practice, the likelihood of the complex measurement depends on the unknown complex amplitude parameters of the targets that must be marginalized. We have shown that closed-form expressions can be obtained in the monotarget case for all the Swerling models. In the multitarget case, a closed-form can be obtained only for the Swerling 1 case; for the others models, we propose some possible approximations to alleviate the computational time and it may be interesting to investigate other approximations that may lead to faster computational time while preserving acceptable performance. We have also considered the case where the data are the squared modulus of the complex measurements. In that case, no closed-form can be obtained and approximations must be performed.

Finally, we have demonstrated via Monte-Carlo simulation the benefits of taking into account the spatial coherence of the complex amplitudes both in detection and in estimation compared methods that works on the square modulus of the complex signal. The main conclusions that can be stated based on this work are the following:

• in a TBD context, complex measurements should be used whenever they are available since it appears that the

Fig. 1 .

 1 Fig. 1. Received signal (noise-free) corresponding to the l th target at two adjacent time steps k and k+1, where dots represent the corresponding measured samples. ρ k,l and ϕ k,l are the same for all cells of z k (we denote this feature spatial coherence) but their values change independently and randomly over time; there is no temporal coherence from step k to step k + 1.

  (51) At step k = 0, we assume that d -1 = 0. December 2, 2014 DRAFT 4) Simulations: For the simulation of the target scenario, the following parameters are used: SNR = 5dB, T=1 s, v min = 100 m/s, v max = 300 m/s, SNR min = 2 dB, SNR max = 10 dB, q s = 10 -3 , P f a = 0.1 and σ 2 n = 0.05. The transition probabilities for the particle filter are set to P b = P d = 0.05. The number of continuing particles is set to N c = 2000 and the number of newborn particles to N b = 1000. Concerning the detection strategy, we choose T h (d k-1 = 0) = 0.9 and T h (d k-1 = 1) = 0.2.

Fig. 2 .

 2 Fig. 2. Monte-Carlo simulation results for the single target case with the Swerling 1 model. average of the detection variable d k . SNR = 5 dB.

Fig. 3 .

 3 Fig. 3. Monte-Carlo simulation results for the single target case with the Swerling 3 model. average of the detection variable d k . SNR = 5 dB.

Fig. 4 .

 4 Fig. 4. Monte-Carlo simulation results for the single target case with the Swerling 1 model. Top: RMSE in position. Bottom: RMSE in velocity. SNR = 5 dB.

Fig. 5 .

 5 Fig. 5. Monte-Carlo simulation results for the single target case with the Swerling 3 model. Top: RMSE in position. Bottom: RMSE in velocity. SNR = 5 dB.

  y0,i x0,i . • The velocity is initialized around the true velocity in Cartesian coordinates by adding a Gaussian noise with December 2, 2014 DRAFT covariance matrix σ 2 v I 2 .

Fig. 6 .Fig. 7 .

 67 Fig. 6. Monte-Carlo simulation results in a multi-target setting with the Swerling 0 model. Top: RMSE in position. Bottom: RMSE in velocity.

Fig. 8 .

 8 Fig. 8. Monte-Carlo simulation results in a multi-target setting with the Swerling 3 model. Top: RMSE in position. Bottom: RMSE in velocity.

  1) depends on phase and amplitude parameters ϕ k,1:N k and ρ k,1:N k that are unknown and not temporally coherent. Therefore, the expression of the measurement likelihood p (z k | X k ) cannot be, in general, directly computed. However it is possible to obtain

	the conditional likelihood p (z k | X k , ρ k,1:N k , ϕ k,1:N k ) that corresponds to a complex Gaussian density with mean
	N k
	µ k =	ρ k,i e jϕ k,i h k,i and covariance matrix Γ:
	i=1
	December 2, 2014	DRAFT

  Nc . It means that complex noise samples n k are mutually independent. Note however that, since modulus ρ k,1:N k and phases ϕ k,1:N k are random variables and spatially coherent at time k, this hypothesis does not allow to establish that signal samples from z k are independent; these samples are independent only conditionally to variables ρ k,1:N k and ϕ k,1:N k . With a slight abuse of notation, we denote by |z k | 2 the vector of squared modulus of the complex signal :

2 I

TABLE I ESTIMATED

 I PROBABILITY OF TRACK LOSS FOR THE DIFFERENT MULTITARGET PARTICLE FILTERS WITH Swerling 0 FLUCTUATIONS.

		Probability of track loss
	CM SW1	1.6 × 10 -2
	CSM SW1	3.1 × 10 -2
	NCSM SW1	4 × 10 -2
	ESM SW1	6.9 × 10 -2

TABLE II ESTIMATED

 II PROBABILITY OF TRACK LOSS FOR THE DIFFERENT MULTITARGET PARTICLE FILTERS WITH Swerling 1 FLUCTUATIONS.

TABLE III ESTIMATED

 III PROBABILITY OF TRACK LOSS FOR THE DIFFERENT MULTITARGET PARTICLE FILTERS WITHSwerling 3 FLUCTUATIONS.

	December 2, 2014	DRAFT

(37) December 2, 2014 DRAFT

& Target Tracking Conference, 2012. December 2, 2014 DRAFT

The noise covariance matrix is assumed to be Γ = 2σ 2 I Nc . Throughout this section, the SNR will be defined by SNR = 10 log 10

. Note that this definition provides a simple relationship between the fluctuation parameters ρ, σ ρ and ν ρ and the target SNR and in the following we will always give the SNR instead of the value of the parameter. Note also that in practice, the ambiguity function have only significant values in a small set of cells around the target location. Therefore, to avoid unnecessary computation, it can be advantageous to truncate the function h(x k ) and compute the likelihood only in the set of cells V x k where the ambiguity function remains significant [START_REF] Salmond | A particle filter for track-before-detect[END_REF]. For a target located in cell (l, m), we define V x k as,

In the rest of the paper, h(x k ) will thus be computed only over the set of cells V x k . Furthermore, note that in that case and with Γ = 2σ 2 I Nc , condition (5) holds when

C. Single Target Simulation and Results

1) Scenario of the Simulation:

We consider a scenario with 100 time steps. The target appears at time step k b = 10 and disappears at step k d = 75. At time step k b , the target state is initialized with the prior distribution defined in section V-A and until time step k d the state is propagated according to (43) (with q s = 0). We also assume that the entire trajectory is contained within area P. The SNR of the target is fixed to 5dB and we consider fluctuations of type Swerling 1 and 3.

2) Single Target Particle Filter: For the simulations, we consider the TBD monotarget particle filter described in [START_REF] Rutten | Recursive track-before-detect with target amplitude fluctuations[END_REF]. In the following, we briefly recall the principle of this filter. To perform the detection stage, a Markovian variable s k is added to the state vector and takes the value 1 if the target is present and 0 otherwise. Thus, the process (s k ) k∈N is entirely defined by its transition probabilities

and

The probability of presence p (s k = 1 | z 1:k ) can be computed recursively and allows to make a decision about the target presence in the radar window. Moreover, as we have already mentioned, parameter 2σ 2 ρ for the Swerling 1 case is unknown and so it is added to the particle state as proposed in [START_REF] Andrieu | Particle methods for change detection, system identification, and control[END_REF]. Obviously this parameter is replaced by ν k,p in the Swerling 3 case. At step k -1, the particle filter consists of a set of N co particles (continuing particles) with uniform weights

that are propagated according to Eq.(43). For the modulus parameter, it is simply propagated according to

where v k is Gaussian noise, with variance σ n .

Then N b new particles (birth particles) are initialized in the radar window. Note that using the uniform prior density for the particle initialization leads to poor performance. Therefore we resort here to the method proposed in [START_REF] Lepoutre | Optimized instrumental density for particle filter in track-before-detect[END_REF] that initializes the particle location (x p k , y p k ) in cells exceeding a given threshold ν P f a = -2σ phases information is very important to improve the performance.

• Multitarget likelihood are not simple to compute except for the particular Swerling 1 case. Thus monotarget likelihood should be performed whenever it is possible to factorize the overall joint density.


domain: math nlin phys.meca.mefl

text:
Jean-Marc Hérard

Hippolyte Lochon

A simple turbulent two-fluid model

We present in this paper a simple turbulent two-phase flow model using the two-fluid approach. The model, that relies on the classical ensemble averaging, allows the computation of unsteady flows including shock waves, rarefaction waves and contact discontinuities. It requires the definition of adequate source terms and interfacial quantities. The hyperbolic turbulent two-fluid model is such that unique jump conditions hold within each field. Closure laws for the interfacial velocity and the interfacial pressure comply with a physically relevant entropy inequality. Moreover, source terms that account for mass, momentum and energy interfacial transfer are in agreement with the entropy inequality. Particular attention is also given to the jump conditions when assuming a perfect gas equation of state within each phase ; this enables to recover expected bounds on the mean density through shock waves.

Introduction

We propose herein a turbulent two-fluid model for the prediction of two-phase flows. Actually, though we use the classical ensemble averaging [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flows[END_REF][START_REF] Drew | Theory of Multicomponent Fluids[END_REF], we follow here the approach used in a recent series of papers. Our main goal is to derive a two-fluid two-phase flow model that accounts for Reynolds stress tensors in a very simple way, so that the fundamental properties invoked in [START_REF] Coquel | Closure laws for a twofluid two-pressure model[END_REF] are preserved. This means that the following specifications are enforced:

• (i) the model should be such that one could retrieve the standard Baer-Nunziato model in the laminar case;

• (ii) an entropy inequality should hold for smooth solutions, and meanwhile it should provide some guidelines for closure laws associated with interfacial mass, momentum and energy transfer;

• (iii) unique jump conditions should be valid so that meaningful and unique shock waves might be predicted.

These requirements are mandatory if we intend to predict relevant shock solutions in two-phase flows, such as those that arise in vapour explosions or other similar situations, while using the standard verification and validation process (see [START_REF] Berry | Relap-7: Demonstrating seven-equation, two-phase flow simulation in a single-pipe, two-phase reactor core and steam Separator/Dryer[END_REF] for instance). A class of two-fluid models that are capable to predict unsteady situations has emerged from the recent literature, either for gas-particle flows (see [START_REF] Baer | A two-phase mixture theory for the deflagrationto-detonation transition (DDT) in reactive granular materials[END_REF][START_REF] Bdzil | Two-phase modeling of deflagration-to-detonation transition in granular materials: A critical examination of modeling issues[END_REF][START_REF] Kapila | Twophase modeling of DDT: structure of the velocity-relaxation zone[END_REF][START_REF] Glimm | Renormalization group solution of two-phase flow equations for rayleigh-taylor mixing[END_REF][START_REF] Jin | Compressible two-pressure two-phase flow models[END_REF] for instance), or for water-vapour flows (see [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Gavrilyuk | Mathematical and numerical modeling of twophase compressible flows with micro-inertia[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF] among others). These models essentially differ from one another through the definition of the interface pressure-velocity couple. We will thus focus here on the approach suggested in [START_REF] Coquel | Closure laws for a twofluid two-pressure model[END_REF], that gave some new enlightment on the admissible closure laws that enable to comply with both items (ii) and (iii) recalled above. As we will see, the constraint (ii) will in turn provide some closure laws for the interfacial pressure (see equation ( 7)), whereas (iii) will enable us to define relevant closure laws for the interface velocity (see equation ( 6)).

Moreover, in order to account for turbulent effects in a rather simple way, we will rely on the single-phase proposal introduced in [START_REF] Hérard | Problème de riemann pour un modèle simple de turbulence monophasique compressible[END_REF], that inherits from the earlier work described in [START_REF] Forestier | Solveur de type godunov pour simuler les écoulements turbulents compressibles[END_REF][START_REF] Berthon | Contribution à l'analyse numérique des équations de Navier-Stokes compressibles à deux entropies spécifiques[END_REF][START_REF] Audebert | Contribution à l'analyse des modèles aux tensions de Reynolds pour l'interaction choc turbulence[END_REF][START_REF] Gavrilyuk | Geometric evolution of the reynolds stress tensor[END_REF]. As we will see, this will minimize the number of unknowns, and inlet/outlet boundary conditions, and meanwhile will allow complying with the above mentionned requirements (i,ii,iii). We also emphasize that the present work should not be confused with the one of [START_REF] Hérard | Numerical modelling of turbulent two phase flows using the two fluid approach[END_REF], which was mainly inspired by [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] and [START_REF] Coquel | Closure laws for a twofluid two-pressure model[END_REF]. Actually, though it agrees with the former items (i) and (ii), the latter turbulent model [START_REF] Hérard | Numerical modelling of turbulent two phase flows using the two fluid approach[END_REF] is not suitable for shock solutions. Indeed, jump conditions are not defined in a unique way, due to the occurence of non-conservative products that are active in genuinely nonlinear fields ; as a result, different mesh-converged solutions issuing from distinct schemes may emerge in practical computations, as it is now quite well-known (see [START_REF] Berthon | Contribution à l'analyse numérique des équations de Navier-Stokes compressibles à deux entropies spécifiques[END_REF][START_REF] Guillemaud | Modélisation et simulation numérique des écoulements diphasiques par une approche bifluide à deux pressions[END_REF]).

The paper is organised as follows. We first provide the governing set of equations including source terms accounting for all possible transfers between phases. The choice of relevant interfacial velocity and pressure is discussed. Then we focus on the key property of the model which is the entropy inequality. It is compared with the laminar case and we underline how turbulent energy affects the different relaxation terms. A third section discusses the main properties of the convective system: hyperbolicity, structure of fields and jump conditions. For simplicity, we restrict ourselves to the Baer-Nunziato closure in this section but few remarks are given about other possible closures. Particular attention is given to the jump conditions when assuming perfect gas equation of state within each phase. In that case, we show that density ratios are bounded through shock waves and that they are in agreement with the laminar frame. Moreover turbulent energy ratios are bounded as well, unlike pressure ratios. The last section is devoted to a few remarks on the Riemann problem. Though the derivation of the model is quite different, we will also address in section 3 the differences and similarities with models proposed in [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF].

Governing equations

The two-phase flow model has been obtained by a statistical averaging of the single-phase Navier-Stokes equations. One additional topological equation on the statistical void fraction is also needed (see [START_REF] Drew | Theory of Multicomponent Fluids[END_REF]). For the sake of simplicity, we do not detail this averaging procedure here but we underline the fact that the tensor involving turbulent effects is modelled following the approach of [START_REF] Hérard | Problème de riemann pour un modèle simple de turbulence monophasique compressible[END_REF] for the single-phase Reynolds stress tensor. Thus the governing set of equations takes the form (when neglecting viscous terms):

∂ t W + ∂ x F (W ) + C(W )∂ x W = S(W ) (1)
with W , F (W ), S(W ) in R 7 , and C(W ) in R 7×7 . The state variable W is:

W =           α 1 α 1 ρ 1 α 2 ρ 2 α 1 ρ 1 u 1 α 2 ρ 2 u 2 α 1 E 1 α 1 E 2           (2)
where α k , ρ k , u k , p k and E k are respectively the statistical void fraction, the mean density, the mean velocity, the mean pressure and the mean total energy of phase k, k = 1, 2. The statistical void factions are such that:

α 2 = 1 -α 1
and the mean total energy E k is given by:

E k = 1 2 ρ k u 2 k + ρ k ε k + K k where ε k = ε k (ρ k , p k )
is the mean internal energy of phase k and K k is the turbulent kinetic energy:

K k = K k,0 ρ 5/3 k (3)
(with K k,0 > 0). It is important to underline the fact that the mean internal energy only depends on the mean density and the mean pressure, which is a crude assumption on the statistical thermodynamics. However, one can easily prove that it is verified in the case of simple Equation Of State (EOS) such as perfect gas or stiffened gas (see [START_REF] Hérard | Une approche bifluide statistique de modelisation des écoulements diphasiques à phases compressibles[END_REF]). We define the set of admissible states Ω by:

Ω = W ∈ R 7 ; α 1 ∈ ]0, 1[ , ρ k > 0, ε k > 0 (4)
We also introduce the celerity of density waves c k in the pure phase k and its temperature T k :

ρ k c 2 k = (∂ p k ε k ) -1 p k ρ k -ρ k (∂ ρ k ε k ) 1 T k = (∂ p k ε k ) -1 (∂ p k s k )
where

s k = s k (ρ k , p k )
is the specific entropy complying with the constraint:

c 2 k (∂ p k s k ) + (∂ ρ k s k ) = 0
The convective part of the system is defined by:

F (W ) =           0 α 1 ρ 1 u 1 α 2 ρ 2 u 2 α 1 (ρ 1 u 2 1 + p 1 + 2 3 K 1 ) α 2 (ρ 2 u 2 2 + p 2 + 2 3 K 2 ) α 1 u 1 (E 1 + p 1 + 2 3 K 1 ) α 2 u 2 (E 2 + p 2 + 2 3 K 2 )           , C(W )∂ x W =           u I ∂ x α 1 0 0 -p I ∂ x α 1 -p I ∂ x α 2 -u I p I ∂ x α 1 -u I p I ∂ x α 2           (5)
where we use the following closure laws on interfacial velocity u I and pressure p I :

u I = au 1 + (1 -a) u 2 , a ∈ 0, m 1 m 1 + m 2 , 1 (6)
p I = b p 1 + 2 3 K 1 + (1 -b) p 2 + 2 3 K 2 , b = 1-a T1 1-a T1 + a T2 (7)
In the case when a ∈ {0, 1}, it corresponds to the so-called Baer-Nunziato closure [START_REF] Baer | A two-phase mixture theory for the deflagrationto-detonation transition (DDT) in reactive granular materials[END_REF] in the laminar case. The third choice a = m1 m1+m2 has been proposed in [START_REF] Coquel | Closure laws for a twofluid two-pressure model[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF] in the laminar case too. Those 3 possibilities for (u I , p I ) have been motivated by two requirements: the enforcement of a relevant entropy inequality which will be discussed later on, and the structure of the field associated with λ = u I which is assumed to be linearly degenerate. Other closure laws could be found [START_REF] Ransom | Hyperbolic two-pressure models for twophase flow[END_REF][START_REF] Glimm | Renormalization group solution of two-phase flow equations for rayleigh-taylor mixing[END_REF][START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Berry | The discrete equation method (DEM) for fully compressible, two-phase flows in ducts of spatially varying cross-section[END_REF]] but they will not be considered here.

The source part of the systeam reads:

S(W ) =           Φ Γ -Γ D + UΓ -D -UΓ Q + UD + HΓ -p I Φ -Q -UD -HΓ + p I Φ           (8) where U = u 1 + u 2 2 , H = u 1 u 2 2
and the relaxation terms are given by:

Φ = α 1 α 2 δ p p 1 + 2 3 K 1 -p 2 + 2 3 K 2 Γ = m 1 m 2 δ µ µ2 T2 -µ1 T1 D = m 1 m 2 δ u (u 2 -u 1 ) Q = m 1 m 2 δ T (T 2 -T 1 ) (9)
with :

µ k = ε k + p k ρ k -T k s k + 5K k 3ρ k
which corresponds to the Gibbs free enthalpy in the laminar case when K k = 0. The scalar functions δ ϕ are not detailed here but are assumed to be positive. For practical purposes, the pressure relaxation time scale involved in δ p is provided by the closure law detailed in [START_REF] Gavrilyuk | The structure of pressure relaxation terms: the one-velocity case[END_REF] ; besides, other relaxation time scales δ µ , δ u , δ T embedded in mass, momentum and energy interfacial transfer terms, are taken from the classical two-fluid literature. Hence, we emphasize the fact that all possible transfers between phases (mass, momentum and energy) are accounted for in this model.

Entropy inequality

A key property of the governing set of equations ( 1) is the following:

Proposition 1. Define the entropy-entropy flux pair (η, f η ):

η = -α 1 ρ 1 s 1 -α 2 ρ 2 s 2 , f η = -α 1 ρ 1 s 1 u 1 -α 2 ρ 2 s 2 u 2
Then the following inequality holds for smooth solutions of (1):

∂ t η + ∂ x f η ≤ 0
Proof. Classical manipulations of the system give us the evolution law of the specific entropy for smooth solutions:

∂ t (α k ρ k s k ) + ∂ x (α k ρ k s k u k ) + 1 T k p I -p k + 2 3 K k (u k -u I ) ∂ x α k = (-1) k+1 T k Q + (U -u k ) D -µ k Γ + p k + 2 3 K k -p I Φ
Then we use the closure laws on interfacial velocity and pressure ( 6)-( 7) to get the evolution law of the global entropy:

∂ t η + ∂ x f η = 1 T 2 - 1 T 1 Q + 1 2T 1 + 1 2T 2 (u 1 -u 2 ) D + µ 1 T 1 - µ 2 T 2 Γ + 1 -b T 1 + b T 2 p 2 + 2 3 K 2 -p 1 + 2 3 K 1 Φ
The conclusion is now obvious using closure laws [START_REF] Forestier | Solveur de type godunov pour simuler les écoulements turbulents compressibles[END_REF]. This entropy inequality is really close to the one in the laminar case, therefore it leads to similar source terms. Turbulent energies are involved in pressure and free enthalpy relaxation terms but one recovers the laminar source terms when turbulent kinetic energies vanish. Temperature and velocity relaxation terms remain exactly the same as in the laminar case.

From now on, we will restrict to the case:

(u I , p I ) = u 1 , p 2 + 2 3 K 2
in the sequel.

Remark 1.

• The model discussed afterwards is in that case quite similar to the one detailed in [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF], pages 298-299, equations (18a-18l), when neglecting the "turbulent entropy" dissipation o χ k (equations 18k-18l), and thus retaining the obvious solution: χ k (x, y, z, t) = (χ k ) 0 , which is actually the straightforward counterpart of equation (3) in the current paper. Note that the contribution o χ k in [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF], page 309, vanishes when the asymptotic pressure-velocity mechanical equilibrium is reached. We also emphasize that the problem arising in the definition of jump relations through shock waves, which is due to the occurrence of non-conservative first-order contributions in the whole system, is not addressed in [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF] ; this will be part of the discussion in the next section (see proposition 4).

• On the other hand, the present model should not be confused with the Discrete Equation Method (DEM), since closure laws involved in the latter for the interface velocity and the interface pressure are indeed totally different (see [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF], equations (33), page 308, to be compared with equations (6), [START_REF] Coquel | Closure laws for a twofluid two-pressure model[END_REF] in the current paper).

Main properties of the convective system

We provide below some of the main properties of the convective system of equations. Propositions 2 and 4 arise as expected.

Proposition 2. The homogeneous convective subset :

∂ t W + ∂ x F (W ) + C(W )∂ x W = 0 ( 10
)
is hyperbolic. It admits seven real eigenvalues:

               λ 1,2 = u 1 λ 3 = u 1 -c1 λ 4 = u 1 + c1 λ 5 = u 2 λ 6 = u 2 -c2 λ 7 = u 2 + c2 with c2 k = c 2 k + 10K k 9ρ k ( 11
)
and associated vectors span the whole space R 7 , unless |u 2 -u 1 | /c 2 = 1. Fields associated with eigenvalues λ 1,2 and λ 5 are linearly degenerate (LD). Other fields are genuinely non linear (GNL).

We notice that the turbulent kinetic energy modifies the celerity of the GNL waves which is straightforward when focusing on single phase turbulent models (see [START_REF] Forestier | Solveur de type godunov pour simuler les écoulements turbulents compressibles[END_REF] for instance). Thus it will also affect resonant waves of the system. Proposition 3 (Riemann invariants). The five Riemann invariants of the 1-2 LD field associated with the void fraction coupling wave are the following:

I 1 1,2 (W ) = u 1 I 2 1,2 (W ) = s 2 I 3 1,2 (W ) = α 2 ρ 2 (u 2 -u 1 ) I 4 1,2 (W ) = α 1 p 1 + 2 3 K 1 + α 2 p 2 + 2 3 K 2 + α 2 ρ 2 (u 2 -u 1 ) 2 I 5 1,2 (W ) = ε 2 + p2 ρ2 + 5 3 K2 ρ2 + 1 2 (u 2 -u 1 ) 2 (12)
The Riemann invariants associated with the other waves read:

I 3 (W ) = α 1 , ρ 2 , u 2 , p 2 + 2 3 K 2 , s 1 , u 1 + φ 1 T I 4 (W ) = α 1 , ρ 2 , u 2 , p 2 + 2 3 K 2 , s 1 , u 1 -φ 1 T I 5 (W ) = α 1 , ρ 1 , u 1 , p 1 + 2 3 K 1 , p 2 + 2 3 K 2 , u 2 T I 6 (W ) = α 1 , ρ 1 , u 1 , p 1 + 2 3 K 1 , s 2 , u 2 + φ 2 T I 7 (W ) = α 1 , ρ 1 , u 1 , p 1 + 2 3 K 1 , s 2 , u 2 -φ 2 T (13)
where

φ k = p k p0 ck ρ k c 2 k (p k , s k )dp k .
Once again, we retrieve the classical results [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF] when the turbulent kinetic energies are not accounted for (thus setting K k = 0). Proposition 4 (Jump conditions). Within each isolated field, unique jump conditions hold. We denote [f ] = f r -f l the jump between the (l)eft and (r)ight states on each side of a discontinuity travelling at speed σ. Turning to the genuinely non linear fields, jump conditions may be written:

               [α k ] = 0 -σ [ρ k ] + [ρ k u k ] = 0 -σ [ρ k u k ] + ρ k u 2 k + p k + 2 3 K k = 0 -σ [E k ] + u k E k + p k + 2 3 K k = 0 (14)
Those jump conditions may be easily rewritten as follows:

         [α k ] = 0 σ = [ρ k u k ] / [ρ k ] (ρ k ) R (ρ k ) L [u k ] 2 = p k + 2 3 K k [ρ k ] 2 ε k + K k ρ k + p k + 2 3 K k r + p k + 2 3 K k l 1 ρ k = 0 (15)
More over, looking at the case of a perfect gas EOS, we have the following:

Proposition 5 (Jump conditions for perfect gas EOS). We assume that phase k complies with the perfect gas EOS:

p k = (γ k -1)ρ k ε k , γ k > 1.
We also assume that the left state is admissible W l ∈ Ω according to (4) and that no vacuum occurs in the solution. Then the jump conditions [START_REF] Glimm | Renormalization group solution of two-phase flow equations for rayleigh-taylor mixing[END_REF] ensure that the right state is admissible: W r ∈ Ω. It also provides bounds for the density ratio whereas the pressure ratio has no bound:

(β k ) -1 < (ρ k ) r (ρ k ) l < β k ( 16
)
with

β k = γ k +1 γ k -1 .
Proof. We note z k = (ρ k ) r /(ρ k ) l the density ratio and π k = (p k ) r /(p k ) l the pressure ratio. The third equation of [START_REF] Guillemaud | Modélisation et simulation numérique des écoulements diphasiques par une approche bifluide à deux pressions[END_REF] gives the following relation between states through genuinely non linear fields associated with λ 3 or λ 6 :

π k (β k -z k ) + 1 -β k z k -g k (z k ) = 0 where g k (z) = 2 3 (K k ) l (p k ) l z 8/3 -4z 5/3 + 4z -1 . Moreover z k > 1 and it implies g k (z k ) > 0.
Then we get the pressure ratio:

π k = β k z k -1 + g k (z k ) β k -z k
and the conclusion is thus straightforward. A similar proof holds through fields associated with λ 4 or λ 7 .

We first note that we recover the bounds arising in the pure laminar case (see for instance [START_REF] Smoller | Shock waves and reaction diffusion equations[END_REF]). This was expected since the instantaneous inequality:

β -1 φ L < φ R < βφ L (17)
(where β is constant with respect to the statistics) obviously provides the following relation between the mean quantities:

β -1 < φ L > ≤ < φ R > ≤ β < φ L > (18)
It is important to notice that, unlike in [START_REF] Forestier | Solveur de type godunov pour simuler les écoulements turbulents compressibles[END_REF][START_REF] Hérard | Numerical modelling of turbulent two phase flows using the two fluid approach[END_REF], we do not require any approximate jump condition. Eventually, we can easily check that the ratio of left and right turbulent kinetic energies is also bounded.

We recall that all the properties described here suppose that the Baer-Nunziato closure holds (u I , p I ) = u 1 , p 2 + 2 3 K 2 . In the case when a = m1 m1+m2 , proposition 2 is slightly different since there is another distinct linearly degenerate field associated with λ = u I , where (u I -u 1 )(u I -u 2 ) = 0. Therefore proposition 3 is also different since other Riemann invariants arise. This implies a more intricate coupling wave associated with the eigenvalue λ = u I , as it already happens in the laminar case (K k = 0).

Remarks

We now make a few remarks on the one-dimensional Riemann problem:

   ∂ t W + ∂ x F (W ) + C(W )∂ x W = 0 W (t = 0, x) = W L , if x < 0 W R , if x > 0 ( 19
)
where W L and W R are admissible states:

W L , W R ∈ Ω.
First of all, we remind that one cannot provide the solution of the general Riemann problem in the laminar case, even when a perfect gas EOS is assumed within each phase (see [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF]). Nonetheless, a simple result may be given, which is the straightforward counterpart of what happens in the laminar case: Remark 2. We assume that (α 1 ) L = (α 1 ) R and that a perfect gas EOS holds within each phase:

p k = (γ k -1)ρ k ε k , γ k > 1.
The Riemann problem (19) has a unique entropy consistent solution involving constant states separated by shocks, rarefaction waves and contact discontinuities if and only if: [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flows[END_REF] for k = 1, 2.

(u k ) R -(u k ) L < (p k ) L 0 ck ρ k c 2 k (p k , s k )dp k + (p k ) R 0 ck ρ k c 2 k (p k , s k )dp k
Proof. Since (α 1 ) L = (α 1 ) R , phases evolve independently. Therefore, we can use the property from [START_REF] Hérard | Problème de riemann pour un modèle simple de turbulence monophasique compressible[END_REF] on the monophasic Riemann problem.

In order to consider real void fraction coupling waves, we have to use the connection through the 1,2-wave given by the Riemann invariants I 1,2 in property 3. This leads to the following: Remark 3. We now suppose that (α 1 ) L = (α 1 ) R . Assume that a perfect gas EOS holds within each phase and that W (x, t) is the solution of the Riemann problem [START_REF] Hérard | Une approche bifluide statistique de modelisation des écoulements diphasiques à phases compressibles[END_REF]. Then the connection of two constant states separated by a simple wave (either a shock, a rarefaction wave or a contact discontinuity) guarantees that all states of W (x, t) are admissible: W (x, t) ∈ Ω , ∀ x, t.

Our last remark addresses the problem of defining K k,0 . At the initial time, we have:

K k,0 = K k ρ 5/3 k (x, t = 0)
Thus, a natural choice of K k,0 using the initial data is the mean value:

K k,0 = 1 |V | V K k ρ -5/3 k (x, t = 0)dx
The turbulent two-phase flow model introduced in this paper is in agreement with requirements (i,ii,iii), and thus allows the computation of shock solutions in turbulent two-phase flows. The extension of the current results to the framework of multiphase flows seems feasable, considering [START_REF] Muller | Closure conditions for nonequilibrium multi-component models[END_REF][START_REF] Hérard | A class of compressible multiphase flow models[END_REF], but this point has not been investigated yet.

Acknowledgments

The authors would like to thank the reviewer for his useful comments. The second author has benefited from financial support of EDF and ANRT through an EDF/CIFRE contract number 2013/0459.