Rong Yang

Shalom Claude

Jauffret Annie-Claude

Pérez Gee

Wah Ng
email: ngeewah@dso.org.sg.

Wah Gee

Ng

Yaakov Bar-Shalom

Claude Jauffret
email: jauffret@univ-tln.fr

Annie-Claude Pérez
email: annie-claude.perez@univ-tln.fr

Maneuvering Target Tracking Using Continuous Wave Bistatic Sonar with Propagation Delay

, which estimates the trajectory of a constant velocity target. The results showed that the estimate bias caused by the propagation delay is not negligible, especially for a bistatic system. This paper develops an interacting multiple model unscented Gauss-Helmert filter with numerical Jacobian (IMM-UGHF-NJ) to track a maneuvering target with propagation delay using a bistatic sonar system. The IMM-UGHF-NJ can overcome the two tracking challenges introduced by the delay, namely, implicit state transition model and lack of analytical expression of the Doppler shifted frequency in the measurement model.

Simulation tests have been conducted, and the results show that the IMM-UGHF-NJ can reduce the estimation error significantly, especially for long range or fast moving targets.

Maneuvering Target Tracking Using

Continuous Wave Bistatic Sonar with Propagation Delay I. INTRODUCTION Continuous active sonar (CAS), also known as high duty cycle (HDC) sonar, with multistatic setup has attracted the research interest recently. In such a system, the signal is transmitted almost in a full duty cycle. Compared to the commonly used pulse active sonar (PAS) system, which transmits only a short pulse in a cycle, the CAS system has continuous detection capability, and is less disturbing to underwater fauna through using a low intensity signal.

There are two main types of CAS systems according to the signal waveforms transmitted, namely, frequency modulated (FM) waveforms and continuous constant frequency waveforms (CW). The FM-CAS can provide good target bistatic range information, whereas the CW-CAS has good Doppler shifted frequency measurement (linked to target range rate). The FM-CAS needs to separate indirect path signal from strong direct path signal via methods, such as m-sequence modulation [START_REF] Deferrari | The Application of m-Sequences to Bi-static Active Sonar[END_REF] and Dopplergram [START_REF] Yang | Acoustic Dopplergram for Intruder Defense[END_REF]. The FM-CAS has a frequency bandwidth limitation issue in multistatic system, as broadband waveforms are transmitted repeatedly [START_REF] Grimmett | Multistatic Tracking for Continous Active Sonar using Doppler-Bearing Measurements[END_REF]. The CW-CAS transmits a single fixed frequency waveform, so that it has no frequency bandwidth limitation problem as the FM-CAS. However, due to lack of range information, the observability of a target trajectory in CW-CAS is not as good as FM-CAS, especially for bistatic (a single transmitter-receiver pair) system.

We focus on target tracking using CW-CAS in this paper. A few approaches on this problem have been proposed in literature before. A Gaussian mixture probability hypothesis density (GMPHD) filter was developed in [START_REF] Grimmett | Multistatic Tracking for Continous Active Sonar using Doppler-Bearing Measurements[END_REF]. It tracks multiple constant velocity (CV) targets using bearings and Doppler frequencies detected by multistatic CW-CAS. Results show that CV targets can be tracked using more than two transmitter-receiver pairs when target range is not available. This research does not take signal propagation delay into consideration. The effect of propagation delay of CW-CAS has been studied in [START_REF] Jauffret | Doppler-only Target Motion Analysis in a High Duty Cycle Sonar System[END_REF] recently. An exact Doppler frequency model with propagation delay was proposed. and a maximum likelihood (ML) estimator based on this model was developed to perform batch estimation for a CV target. The simulation results showed that the estimation bias induced by the propagation delay is not negligible, especially for a bistatic system.

In this paper, the propagation delay problem raised in [START_REF] Jauffret | Doppler-only Target Motion Analysis in a High Duty Cycle Sonar System[END_REF] is studied further. We extend the target CV trajectory estimation using a batch parameter estimation technique to the dynamic recursive estimation, which can handle not only CV motion but also maneuvering motion. This extension faces two challenges. Firstly, the "target time" t k and target position (x k ,y k ) in the state [defined later in [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation: Theory, Algorithms and Software[END_REF] and [START_REF] Bar-Shalom | Tracking and Data Fusion: A Handbook of Algorithms[END_REF]] are highly correlated after propagation delay is introduced. This leads to a state transition equation in an implicit form instead of the commonly used explicit form in [START_REF] Yang | Interacting multiple model unscented Gauss-Helmert filter for bearings-only tracking with state-dependent propagation delay[END_REF] [START_REF] Yang | UGHF for acoustic tracking with state-dependent propagation delay[END_REF]. Secondly, the Doppler shifted frequency (one of the measurements) does not have an analytical expression in terms of the target state. This is because the Doppler frequency is a function of the bistatic range rate which cannot be described analytically after propagation delay is introduced. Details will be given later in Section II-B. The two challenges mentioned above were overcome in [START_REF] Jauffret | Doppler-only Target Motion Analysis in a High Duty Cycle Sonar System[END_REF] through solving a 2nd order polynomial equation for CV target. However, the approach in [START_REF] Jauffret | Doppler-only Target Motion Analysis in a High Duty Cycle Sonar System[END_REF] cannot be applied to a maneuvering target with coordinated turn (CT) motion, and the new approach in this paper will be shown to handle this.

A dynamic estimation problem uses two basic models, namely, the state transition model and the measurement model. The state transition model describes the evolution of the target state with time, and it is (in most cases) an explicit expression of the state at the current time in terms of the state at the previous time. The measurement model relates the measurement to the state. The two challenges of the dynamic estimation problem considered in this paper are: (i) the implicit state transition model;

(ii) the lack of an analytical measurement model. These make this problem impossible to solve using existing filters.

We will develop an interacting multiple model unscented Gauss-Helmert filter with numerical Jacobian (IMM-UGHF-NJ) to cope with the challenges mentioned above. The IMM [START_REF] Bar-Shalom | Tracking and Data Fusion: A Handbook of Algorithms[END_REF] is a well known hybrid algorithm to handle motion model uncertainty in maneuvering target tracking. The UGHF [11][12][14] [START_REF] Yang | Bearings-Only Tracking with Fusion from Heterogenous Passive Sensors: ESM/EO and Acoustic[END_REF] is a recently developed algorithm for bearings-only tracking (BOT) with implicit state transition model introduced by the acoustic propagation delay. It can be applied to our problem. For the measurement model without analytical form, the NJ (numerical Jacobian) algorithm, which computes the Jacobian numerically, can be utilized [START_REF] Xiong | The linear fitting Kalman filter for nonlinear tracking[END_REF][13] [START_REF] Xiong | Linear fitting Kalman filter[END_REF]. The Doppler shifted frequency is a function of the bistatic range rate, ṙ, which has no analytical form due to the unknown time delay. We can compute ṙ (derivative of range r) using the NJ.

The structure of the rest of paper is as follows. Section II formulates the problem. Section III presents the IMM-UGHF-NJ. Simulation results and conclusions are in Sections IV and V, respectively.

II. PROBLEM FORMULATION

The problem is illustrated in Fig. 1. At dynamic estimation cycle k, the transmitter emits a CW signal with constant frequency f T at time t T k , and the receiver receives the Doppler shifted frequency f R at time t R k via the target reflection at time t k . We assume the transmitter and receiver are stationary and located at (x T ,y T ) and (x R ,y R ), respectively. The target is moving and its location is [x(t k ),y(t k )] at reflection time t k . The ranges between the target at t k to the transmitter and the receiver are r T k and r R k , respectively. The target states to be estimated for the CV and CT models at time t k are

x CV (t k ) = [x(t k ) y(t k ) ẋ(t k ) ẏ(t k ) t k ] (1)
x

CT (t k ) = [x(t k ) y(t k ) ẋ(t k ) ẏ(t k ) ω(t k ) t k ] (2)
where x, y, ẋ and ẏ are the target positions and velocities in the x and y coordinates, respectively, ω is the target turn rate, and t k is the target time (or reflection time) corresponding to the emission time t T k and the reception time t R k of the transmitter and receiver, respectively. The measurement vector at time

t R k is z(t R k ) = [b(t R k ) f R (t R k )] (3)
where b is the target bearing from the receiver at time t R k to the target at time t k , measured clockwise from True North, and f R is the Doppler shifted frequency at the receiver.

A. State transition models

The state transition model describes the evolution of the target state with time. For a generic discrete problem, it is an explicit form given by

x(t k ) = f [x(t k-1 )] + Γv(t k-1 ) (4)
where k is the discrete estimation cycle index, v(t k-1 ) is the process noise, and Γ is the process noise gain. However, there is no explicit state transition model for our problem. It can be seen from December 20, 2016 DRAFT Fig. 2 that the target time, t k , is unknown due to the unknown propagation delay τ R . There is an implicit constraint between the known t R k and unknown t k given by

t k = t R k -τ R k ( 5
)
where

τ R k = [x(t k ) -x R ] 2 + [y(t k ) -y R ] 2 c p (6)
and c p is the signal propagation speed in the medium. It can be seen that t k is on the both sides of the constraint equation ( 5), since x(t k ) and y(t k ) are functions of t k . It is difficult to obtain an explicit express of t k . This leads to use a Gauss-Helmert (GH) state transition model, which describes an implicit constraint systemically [11][12]. The GH model is given by

g[x(t k ), x(t k-1 )] + Γv(t k-1 ) = 0 (7)
The GH models for the CV motion 1 and CT motion are given next.

1) Constant velocity Gauss-Helmert model

The GH model for CV motion is given by

g CV x CV (t k ), x CV (t k-1 ) + Γ CV v CV (t k-1 ) = 0 5 (8)
where g CV [•] is the implicit GH state transition function, which combines the CV motion constraints and the delay constraint between x(t k ) and x(t k-1 ). It is given by

g CV (•) = g CV 1 (•) g CV 2 (•) g CV 3 (•) g CV 4 (•) g CV 5 (•) (9)
where

g CV 1 (•) = x(t k ) -[x(t k-1 ) + ẋ(t k-1 )∆ k ] (10)
g CV 2 (•) = y(t k ) -[y(t k-1 ) + ẏ(t k-1 )∆ k ] (11)
g CV 3 (•) = ẋ(t k ) -ẋ(t k-1 ) (12)
g CV 4 (•) = ẏ(t k ) -ẏ(t k-1 ) (13)
g CV 5 (•) = t k -(t R k -τ R k ) (14)
with τ R k given in [START_REF] Jauffret | Doppler-only Target Motion Analysis in a High Duty Cycle Sonar System[END_REF] and

∆ k = t k -t k-1 (15)
Based on the discrete white noise acceleration (WNA) model [START_REF] Bar-Shalom | Tracking and Data Fusion: A Handbook of Algorithms[END_REF], the gain matrix Γ CV and the zeromean white Gaussian process noise v CV in (8) compensate for small accelerations and the uncertainty of the sound speed. The noise gain matrix Γ CV is given by

Γ CV =             1 2 (∆ k ) 2 0 0 0 1 2 (∆ k ) 2 0 ∆ k 0 0 0 ∆ k 0 0 0 1             (16)
The covariance of v CV is

q CV = diag(σ 2 ẍ σ 2 ÿ σ 2 t ) (17)
where σ 2 ẍ and σ 2 ÿ are the variances on small target accelerations in the x and y coordinates respectively, and σ 2 t is the process noise variance on the target time. The covariance of the error in the model ( 8) is given by

Q CV (∆ k ) = Γ CV q CV (Γ CV ) (18)
2) Coordinated Turn Gauss-Helmert model

The GH state transition model for the CT motion is given by

g CT [x CT (t k ), x CT (t k-1 )] + Γ CT v CT (t k-1 ) = 0 6 (19)
where

g CT (•) = g CT 1 (•) g CT 2 (•) g CT 3 (•) g CT 4 (•) g CT 5 (•) g CT 6 (•) (20)
with

g CT 1 (•) = x(t k ) -x(t k-1 ) + sin[ω(t k-1 )∆ k ] ω(t k-1 ) ẋ(t k-1 ) - 1 -cos[ω(t k-1 )∆ k ] ω(t k-1 ) ẏ(t k-1 ) (21) g CT 2 (•) = y(t k ) -y(t k-1 ) + sin[ω(t k-1 )∆ k ] ω(t k-1 ) ẏ(t k-1 ) + 1 -cos[ω(t k-1 )∆ k ] ω(t k-1 ) ẋ(t k-1 ) (22)
g CT 3 (•) = ẋ(t k ) -{cos[ω(t k-1 )∆ k ] ẋ(t k-1 ) -sin[ω(t k-1 )∆ k ] ẏ(t k-1 )} (23) g CT 4 (•) = ẏ(t k ) -{sin[ω(t k-1 )∆ k ] ẋ(t k-1 ) + cos[ω(t k-1 )∆ k ] ẏ(t k-1 )} (24) g CT 5 (•) = ω(t k ) -ω(t k-1 ) (25)
g CT 6 (•) = t k -(t R k -τ R k ) (26)
The noise gain matrix Γ CT is given by

Γ CT =                 1 2 (∆ k ) 2 0 0 0 0 1 2 (∆ k ) 2 0 0 ∆ k 0 0 0 0 ∆ k 0 0 0 0 ∆ k 0 0 0 0 1                 (27)
q CT = diag(σ 2 ẍ σ 2 ÿ σ 2 ω σ 2 t ) (28)
where σ 2 ω is the variance of the Gaussian process noises of ω. The covariance of the error in (19) for the (nearly) CT motion, Q CT (∆ k ), is computed by

Q CT (∆ k ) = Γ CT q CT (Γ CT ) (29)

B. Measurement model

The measurement model relates the state at time t k to the measurement at time t R k , which is given by

z(t R k ) = h[x(t k )] + w(t R k ) (30)
where w(t R k ) is the measurement noise, and

h(•) = [h 1 (•) h 2 (•)] (31)
with

h 1 (•) = b(t R k ) = tan -1 x(t k ) -x R y(t k ) -y R (32) h 2 (•) = f R (t R k ) = f T (t T k ) 1 - ṙ(t R k ) c P ( 33
)
The challenge is how to obtain ṙ(t R k ) in (33). We know

r(t R k ) = r T k + r R k = [x(t k ) -x T ] 2 + [y(t k ) -y T ] 2 + [x(t k ) -x R ] 2 + [y(t k ) -y R ] 2 (34)
and

ṙ(t R k ) = d[r(t R k )] d(t R k ) = ẋ(t k )[x(t k ) -x T ] + ẏ(t k )[y(t k ) -y T ] r T k dt k d(t R k ) + ẋ(t k )[x(t k ) -x R ] + ẏ(t k )[y(t k ) -y R ] r R k dt k d(t R k ) (35)
When the signal propagation delay is negligible (for example, for a radar signal), one has

t k = t R k and dt k d(t R k ) = 1 (36)
The analytical form of

ṙ(t R k ) is then ṙ(t R k ) = ẋ(t k )[x(t k ) -x T ] + ẏ(t k )[y(t k ) -y T ] r T k + ẋ(t k )[x(t k ) -x R ] + ẏ(t k )[y(t k ) -y R ] r R k (37)
However, the acoustic signal in our problem has significant propagation delay and t k = t R k . The analytical function

t k = f (t R k ) (38)
is impossible to obtain for a target in CT motion. This causes a major challenge for mapping the state to the measurement. An appropriate filter to cope with this challenge will be developed next.

III. INTERACTING MULTIPLE MODEL UNSCENTED GAUSS-HELMERT FILTER WITH NUMERICAL

JACOBIAN

The IMM estimator [START_REF] Bar-Shalom | Tracking and Data Fusion: A Handbook of Algorithms[END_REF] is the most commonly used hybrid approach to handle model uncertainty in target tracking. This section describes an IMM-UGHF-NJ filter with the implicit CV and CT models described in Section II and lack of analytical expression for the measurement function.

Similarly to the original IMM estimator, the IMM-UGHF-NJ performs the state estimation in four steps: mixing, mode-matched filtering, mode probabilities updating and final state combination:

1) In the mixing step, the m hypotheses (where m is the number of models in the filter) at time k-1 expand to m 2 hypotheses using the mixing probabilities based on the mode Markov chain, which is governed by the m × m mode probability transition matrix Π consisting of the mode transition probabilities, p ij . The m 2 hypotheses are then merged into m hypotheses based on the mixture equations [START_REF] Bar-Shalom | Tracking and Data Fusion: A Handbook of Algorithms[END_REF].

2) In the mode-matched filtering step, the mixed state estimates are updated by UGHF-NJs (given later) in parallel.

3) The mixing probabilities are obtained, and the updated mode probabilities are computed based on the innovations in the mode-matched UGHF-NJs. The updated mode probabilities together with the mode-conditioned estimated states and covariances are brought to the next step.

4) The final state estimate and its covariance for the current time cycle are computed based on the mixture equations using the latest mode probabilities in the combination step.

Since the states in the CV and CT models described in Section II have different dimensions, the unbiased mixing approach [START_REF] Yuan | A multiple IMM estimation approach with unbiased mixing for thrusting projectiles[END_REF] is applied in the IMM filter to increase the CV state from 5 to 6.

Before the mixing step, the CV state estimate and its error covariance are augmented with the turn rate information from the CT model.

The IMM-UGHF-NJ differs from the standard IMM in the mode-matched filters, which are UGHF-NJ. The UGHF-NJ handles the implicit GH state transition model and evaluates f R (t R k ) in the measurement vector (3) numerically. The UGHF-NJ prediction, state-to-measurement mapping and update steps are given in Algorithms 1-3, respectively. In these algorithms, the model superscripts "CV" and "CT" for the states and GH functions are omitted for simplicity.

Algorithm 1 predicts the state x( tk-1 ) from time tk-1 to an unknown target time, t k , corresponding to the signal reception time t R k . The relationship between t k and t R k is given by the implicit constraint [START_REF] Grimmett | Multistatic Tracking for Continous Active Sonar using Doppler-Bearing Measurements[END_REF]. An unscented Gauss-Helmert approach is used for the state prediction with the implicit constraint. Firstly, 2n x + 1 sigma points of x( tk-1 ) are generated using SigPt(•) (given in the December 20, 2016 DRAFT Algorithm 1 UGHF-NJ prediction Generate (2n x + 1) sigma points for x(t k-1 ):

[{x i ( ti k-1 )}, {w i }] = SigPt[x( tk-1 ), P( tk-1 ), κ]
Predict sigma points using Gauss-Newton algo.:

for all xi ( ti k-1 ), i ∈ {1, . . . , 2n x + 1} do x 0 = xi ( ti k-1 ) xi ( ti k | ti k-1 ) = GaussN[g(x 1 , x 0 )] end for
Regen sigma points with process noise:

x( tk | tk-1 ) = 2nx+1 i=1 w i xi ( ti k | ti k-1 ) P( tk | tk-1 ) = 2nx+1 i=1 w i xi ( ti k | ti k-1 )(x i ( ti k | ti k-1 )) + Q(∆ k ) [{x i ( ti k | ti k-1 )}, {w i }] = SigPt[x( tk | tk-1 ), P( tk | tk-1 ), κ] where xi ( ti k | ti k-1 ) = xi ( ti k | ti k-1 ) -x( tk | tk-1 )
κ is a spread scalar of the sigma points. Appendix), where n x is the dimension of the state vector. Secondly, each sigma point is predicted to ti k using the Gauss-Newton algorithm GaussN(•) (also given in the Appendix) based on the Gauss-Helmert function g(x 1 , x 0 ), where i is the index of the sigma points. The 2n x + 1 GaussN(•) find

x 1 = xi ( ti k | ti k-1 ) from x 0 = xi ( ti k-1
) iteratively. Thirdly, the predicted sigma points are re-generated with considering also the process noise (with the approprate larger prediction covariance).

Algorithm 2 maps the predicted state to the measurement space. The challenge here is that we cannot obtain the Doppler shifted frequency f R (t R k ) in the measurement from the predicted state directly. The range rate ṙ(t R k ) in (33) cannot be derived from the bistatic range r(t R k ), which has no analytical form in terms of t R k . We use an numerical approach, called numerical Jacobian (NJ), to obtain ṙ(t R k ) from r(t R k ). It is known that the slope of the tangent line is the derivative of a nonlinear function at a point of interest. The principle of the NJ(•) (given in the Appendix) is to find the best linear fit to a nonlinear function based on a few weighted points around the point of interest. If we can provide these weighted points around [t R k ,r(t R k )], its derivative ṙ(t R k ) can then be computed using NJ(•). Firstly, we generate the reception time set around t R k using SigPt(•), i.e.,

{t R,j k } = {t R k , t R k -σ t R k , t R k + σ t R k } j = 1, 2, 3 (39)
Algorithm 2 UGHF-NJ mapping the predicted state to measurement

[{t R,j k }, {w j }] = SigPt[t R k , σ t R k , κ] for all xi ( ti k | ti k-1 ), i ∈ {1, . . . , 2n x + 1} do x 0 = xi ( ti k | ti k-1 ) for j = 1 : 3 do xi,j ( tj k | ti k-1 ) = GaussN g(x 1 , x 0 )| t R k =t R,j k ri,j (t R k ) ← xi,j ( tj k | ti k-1 ) end for ri (t R k ) = NJ[{t R,j k }, {r i,j (t R k )}, {w j }] f R,i (t R k ) ← using (33) bi (t R k ) ← using (32) z i (t R k ) = [ bi (t R k ) f R,i (t R k )] end for ẑ(t R k ) = 2nx+1 i=1 w i ẑi (t R k )
where σ t R k is a very small shift from t R k . Its weight set is {w j }. Secondly, we use GaussN(•) to obtain the predicted state set {x i,j ( tk | tk-1 )} corresponding to the reception time set {t R,j k } for the ith sigma point of the predicted state (obtained from Algorithm 1). The bistatic range can then be computed using (34). The set of bistatic ranges corresponding to {t R,j k } for the ith sigma point of the predicted state is

{r i,j (t R k )} = {r i (t R k ), ri (t R k -σ t R k ), ri (t R k + σ t R k )} j = 1, 2, 3 (40)
Thirdly, we use these two sets, {t R,j k } and {r i,j (t R k )}, which form three points around

[t R k , ri (t R k )] to evaluate the range rate ri (t R k ) using NJ(•). Once ri (t R k ) is obtained, f R,i (t R k
) can be computed using (33), and the predicted measurement z i (t R k ) follows. Algorithm 3 updates the predicted state based on the measurment z(t R k ). This step is the same as in the conventional UKF.

IV. SIMULATION RESULTS

The IMM-UGHF-NJ is tested with simulated data in this section. The simulated scenarios are shown in Fig. 3. Twelve targets move in CV-CT-CV motion with different speeds and ranges. They are categorised into four groups based on the ranges (or distances) to the transmitter and receiver, which are between 0-5km, 5-10km, 10-15km and 15-20km. Each category has three targets with speeds 10m/s, 20m/s and 30m/s, respectively. All targets have two CV legs linked by a CT arc. The durations of the first CV, CT and the second CV are 90s, 45s and 90s, respectively. The CT arc is a 90 o right turn with turn rate 2 o /s. The transmitter and receiver are located at (-3500m,0m) and December 20, 2016 DRAFT

Algorithm 3 UGHF-NJ update x( tk ) = x( tk | tk-1 ) + K k ν(t R k ) P( tk ) = P( tk | tk-1 ) -K k S(t R k )K k where ν(t R k ) = z(t R k ) -ẑ(t R k ) K k = P xz S(t R k ) -1 S(t R k ) = R + P zz P xz = 2nx+1 i=1 w i xi ( ti k | ti k-1 )z i (t R k ) P zz = 2nx+1 i=1 w i [z i (t R k )z i (t R k ) ] zi (t R k ) = ẑi (t R k ) -ẑ(t R k ) xi ( ti k | ti k-1 ) = xi ( ti k | ti k-1 ) -x( tk | tk-1 )
(3500m,0m), respectively. The transmitter emits a CW signal with frequency 1000Hz. The sampling interval of the receiver is T = 1s. The measurement errors of bearing and Doppler shifted frequency at receiver are assumed Gaussians with standard deviations σ b = 1 o and σ f = 0.25Hz, respectively.

The sound propagation speed in water is c p = 1484m/s. The following two algorithms are used in testing:

• IMM-UKF: The mode-matched filters are UKF. They estimate target position and velocity only.

The propagation delay is not taken into consideration at all. The Doppler shifted frequency in the measurement model is based on (37) which is commonly used in multistatic radar tracking December 20, 2016 DRAFT system. The target times are taken as the signal reception times by the receiver.

• IMM-UGHF-NJ: This is the new algorithm proposed in this paper. The propagation delay is taken into consideration in the state estimation, and the target times attached to the target trajectory are estimated from multiple UGHF-NJs.

One CV model and two CT models (CT-L and CT-H) are used in both IMM estimators. The CT-L and CT-H have low and high turn rate process noises, respectively. This setup can provide a fast turn rate adaptation during model switching [START_REF] Huang | Tracking a maneuvering target using two heterogeneous passive sensors on a single stationary platform with IMM estimation[END_REF]. The initial mode probabilities for the three models are 1/3. The probability transition matrix Π 3 is

Π 3 =      
0.950 0.025 0.025 0.025 0.950 0.025 0.025 0.025 0.950

      (41)
The measurement error covariance R is

R = diag[(1 o ) 2 (0.25Hz) 2 ] ( 42
)
In the IMM-UGHF-NJ, the process noise covariances q CV , q CT-L and q CT-H are, respectively,

q CV = diag[(0.1m/s 2 ) 2 (0.1m/s 2 ) 2 (0.1s) 2 ] ( 43
)
q CT-L = diag[(0.1m/s 2 ) 2 (0.1m/s 2 ) 2 (0.1 o ) 2 (0.1s) 2 ] ( 44
)
q CT-H = diag[(0.1m/s 2 ) 2 (0.1m/s 2 ) 2 (1 o ) 2 (0.1s) 2 ] ( 45
)
and κ is set to 1 in all SigPt(•) (see the Appendix), and σ t R k is set to 0.1s in Algorithm 2. The initial state estimates are

xCV (t 0 ) = [r 0 sin b 0 r0 cos b 0 x0 ŷ0 t0 ] (46) xCT-L (t 0 ) = [r 0 sin b 0 r0 cos b 0 x0 ŷ0 0.1 o t0 ] (47) xCT-H (t 0 ) = xCT-L (t 0 ) (48)
where

r0 ∼ N (r R 0 , σ 2 r ) (49) b 0 = b(t R 0 ) (50) x0 ∼ N ( ẋ0 , σ 2 ẋ) (51) ŷ0 ∼ N ( ẏ0 , σ 2 ẏ ) (52) t0 = t R 0 -r0 /c p (53)
with r R 0 the true value of the range from the target at time t 0 to the receiver at time t R 0 , σ r = 400m, and b(t R 0 ) is the measured bearing at time t R 0 , ẋ0 and ẏ0 are the true target velocities, and σ ẋ = σ ẏ = 4m/s. The initial state error covariances for the three models are

P CV (t 0 ) =            
P xx P xy 0 0 0

P yx P yy 0 0 0 0 0 σ 2 ẋ 0 0 0 0 0 σ 2 ẏ 0 0 0 0 0 (σ r /c p ) 2            
(54)

P CT-L (t 0 ) = P CT-H (t 0 ) =                
P xx P xy 0 0 0 0

P yx P yy 0 0 0 0 0 0 σ 2 ẋ 0 0 0 0 0 0 σ 2 ẏ 0 0 0 0 0 0 (0.02 o ) 2 0 0 0 0 0 0 (σ r /c p ) 2                 (55)
where

P xx = (r 0 σ b cos b 0 ) 2 + (σ r sin b 0 ) 2 (56)
P yy = (r 0 σ b sin b 0 ) 2 + (σ r cos b 0 ) 2 (57)
P xy = P yx = (σ 2 r -r2 0 σ 2 b ) sin b 0 cos b 0 (58)
The parameters in the IMM-UKF, including the process noise covariances, initial states and initial state error covariances are the same as the IMM-UGHF-NJ, but the elements corresponding to the target time are removed.

The simulation results present the root mean square errors (RMSE) of the estimated target positions and speeds obtained from 100 Monte Carlo runs. The estimated position and speed errors at time tk are computed by

pos err ( tk ) = [x( tk ) -x( tk )] 2 + [ŷ( tk ) -y( tk )] 2 (59) sp err ( tk ) = [ x( tk ) -ẋ( tk )] 2 + [ ŷ( tk ) -ẏ( tk )] 2 (60)
where x( tk ), ŷ( tk ), x( tk ) and ŷ( tk ) are the estimated target positions and velocities in the x and y coordinates respectively, x( tk ), y( tk ), ẋ( tk ) and ẏ( tk ) are the true target positions and velocities in December 20, 2016 DRAFT the x and y coordinates respectively, and tk is the estimated target time in estimation cycle k.

Tables I and II show the averages of position and speed RMSE for the two algorithms for the twelve simulated targets from the four categories displayed in Fig. 3. Figs. 4-7 show the position RMSE versus time of the two algorithms for four simulated targets, one from each category, respectively.

They are the targets in the range between 0-5km with speed 30m/s, range between 5-10km with speed 10m/s, range between 10-15km with speed 20m/s and range between 15-20km with speed 30m/s. It can be seen that the IMM-UGHF-NJ outperforms the IMM-UKF for all targets. The accuracy improvement is target range and speed dependent. A faster and longer range target has more improvement than a slower one at a shorter range. This is because that estimation error of the IMM-UKF depends on the target speed and propagation delay τ R k (details can be found in Section V-C of [START_REF] Yang | UGHF for acoustic tracking with state-dependent propagation delay[END_REF]). The range from the target to the receiver is proportional to the propagation delay. From the results we can say that the estimation error without considering propagation delay is significant, especially for a long range target or a fast target (such as a speed boat or torpedo). The maneuvering mode probabilities of the two IMM filters are also investigated. Figs. 891011show the sum of the mode probabilities of the two CT models (which represents the target maneuvering probability) versus time for the four targets, respectively. It can be seen that the maneuvering probability for both filters increases when the target is maneuvering. The IMM-UGHF-NJ reacts faster than the IMM-UKF. A delay in the model switching for a long range target (> 5km) is observed.

However, the mode probability does not match the ground truth very well when the target is in CV motion. This is because the turn rate ω in CT models can adapt to a small value when the target is in CV motion. To evaluate the consistency of the IMM-UGHF-NJ and IMM-UKF, the average normalized estimation error squared (NEES) is evaluated. The average (2D) position NEES at time tk for N Monte Carlo runs is [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation: Theory, Algorithms and Software[END_REF]

¯ ( tk ) = 1 2N N i=1 xi 1:2 ( tk ) [P i 1:2,1:2 ( tk )] -1 xi 1:2 ( tk ) (61)
where i the run index, P which are around 90s and 135s. Another one is at the ending part of the near range target (0-5km, 30m/s). When the target is switching between the CV and CT motions, the IMM-UGHF-NJ cannot adapt to the correct model immediately, and this causes short delay in the maneuver start and maneuver end, but these delays are shorter than for the IMM-UKF. For the near range target (0-5km, 30m/s), the NEES is below the lower bound 0.81 at the ending part (t > 160s). We can observe from Fig. 8 that the maneuvering probability is greater than 0.24 when t >160s. It is apparently worse than for the other three targets shown in Figs. 9-11. This is caused by the marginal observability of the CV motion model from the measurements, and results in the maneuvering probability (sum of the probabilities of CT models) not small enough. The error covariance of the combined estimate is too large (pessimistic) when the contribution of the incorrect models (the maneuvering models with probability around 0.25) cannot be overlooked. The small NEES is therefore caused by this large error covariance. a) SigPt(•) generates the sigma points for a random variable x with covariance P x [START_REF] Julier | A new extension of the Kalman filter to nonlinear systems[END_REF].

[x i , w i ] = SigPt(x, P x , κ) i = 1, . . . , 2n x + 1 (64)

where

x 1 = x (65)
x i = x + (n x + κ)P x i-1 i = 2, • • • , n x + 1 (66)
x i = x -(n x + κ)P x i-nx-1 i = n x + 2, • • • , 2n x + 1 (67)

w 0 = κ n x + κ i = 1 (68)
w i = 1 2(n x + κ) i = 2, • • • , 2n x + 1 (69)
where n x is the dimension of x, (n x + κ)P

where x 0 is known. The iteration procedure is

xj+1 1 = xj 1 + (A j ) -1 g(x j 1 , x 0 ) ( 71
)
where j is the iteration index, A j is the Jacobian matrix defined by A j = ∂g[(x j 1 , x 0 )] ∂x j 1 (72) c) NJ(•) calculates the Jacobian (or derivative) H of a function

z = h(x) (73)
at a point of interest x 0 numerically [START_REF] Xiong | The linear fitting Kalman filter for nonlinear tracking[END_REF][13] [START_REF] Xiong | Linear fitting Kalman filter[END_REF]. There is no analytical form for h(•), but z can be obtained through numerical method from a given x. The Jacobian is

H = NJ[{x i }, {z i }, {w i }] (74)
where {x i } is the sigma point set around x 0 generated from a very small covariance, {z i } is its corresponding set after transformation and {w i } is the set of weights. The NJ is implementing through the following steps:

1) Form the sigma point set

X =    x 1 x 2 • • • x 2nx+1 1 1 • • • 1    -    x 1 0    (75)
Z =          z1 z2 . . . zl          = [ z 1 z 2 • • • z 2nx+1 ] (76)
where x 1 = x 0 , and l is the dimension of z.

2) Estimate H using the weighted least squares (WLS) algorithm a j = ( XW X ) -1 XW(z j ) (77)

Ĥ = a 1 a 2 • • • a l (78) Ĥ = Ĥ(1 : l, 1 : n x ) (79)
where W = diag({w i }), j ∈ {1, . . . , l}, and Ĥ is Ĥ without the last column.

Fig. 1 .

 1 Fig. 1. Signal transmission of CW bistatic sonar.

Fig. 2 .

 2 Fig. 2. Time sequences of continuous wave bistatic sonar.

Fig. 3 .

 3 Fig. 3. Test scenarios.

Fig. 4 .

 4 Fig. 4. Position estimate RMSE versus time for the target with speed=30m/s and range less than 5km.

i 1 : 2 , 1 : 2 Fig. 5 .Fig. 6 .

 121256 Fig. 5. Position estimate RMSE versus time for the target with speed=10m/s and range 5-10km.

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Position estimate RMSE versus time for the target with speed=30m/s and range 15-20km.

Fig. 10 .Fig. 11 .Fig. 12 .

 101112 Fig. 10. Maneuvering probability versus time for target with speed=20m/s and range 10-15km.

Fig. 13 .

 13 Fig. 13. Four target position NEES versus time for the IMM-UKF.

Fig. 12

 12 Fig.12shows the average position NEES versus time of the IMM-UGHF-NJ for the four targets

Fig. 13 Fig. 14 .Fig. 15 . 2 (

 1314152 Fig. 13 shows the NEES of the IMM-UKF for the same four targets. All of them are above the upper bound 1.21. Obviously, the IMM-UKF provides biased estimation without considering propagation delay.

  x i * indicates the i * th column of the matrix [•], and κ is a scalar that determines the spread of sigma points. b) GaussN(•) is a Gauss-Newton algorithm to obtain the solution of an implicit equation g(•) = 0 iteratively [11][12] and yields x1 = GaussN[g(x 1 , x 0 )]

December 20, 2016 DRAFT

Although an explicit state transition model for the CV motion can be obtained through solving a

2nd order polynomial equation[START_REF] Jauffret | Doppler-only Target Motion Analysis in a High Duty Cycle Sonar System[END_REF], the GH model is a systematical way which is suitable for both CV and CT motions.December 20, 2016 DRAFT


C Dombry

N Guillotin-Plantard

A functional approach for random walks in random sceneries

A functional approach for the study of the random walks in random sceneries (RWRS) is proposed. Under fairly general assumptions on the random walk and on the random scenery, functional limit theorems are proved. The method allows to study separately the convergence of the walk and of the scenery: on the one hand, a general criterion for the convergence of the local time of the walk is provided, on the other hand, the convergence of the random measures associated with the scenery is studied. This functional approach is robust enough to recover many of the known results on RWRS as well as new ones, including the case of many walkers evolving in the same scenery.

Introduction

Let S = (S n ) n∈N be a random walk in Z d starting from zero. The term "random walk" has to be understood in the broad sense of a discrete process, not necessarily with independent and identically distributed (i.i.d.) increments. Let ξ = {ξ x , x ∈ Z d } be a collection of random variables called the random scenery. The random walk in random scenery (RWRS) is the cumulative process

Z n = n k=0 ξ S k .
Motivated by the construction of a new class of self-similar stationary increments processes, Kesten and Spitzer [START_REF] Kesten | A limit theorem related to a new class of self-similar processes[END_REF] and Borodin [START_REF] Borodin | Limit theorems for sums of independent random variables defined on a nonrecurrent random walk. Investigations in the theory of probability distributions[END_REF][START_REF] Borodin | Limit theorems for sums of independent random variables defined on a recurrent random walk[END_REF] introduced RWRS in dimension one and proved functional limit theorems. More precisely, they considered the case when the random walk has i.i.d. increments belonging to the domain of attraction of a stable distribution with index α ∈ (1, 2], the random scenery is i.i.d. and belongs to the domain of attraction of a stable distribution with index β ∈ (0, 2]; the scenery and the walk are independent. Under these assumptions, the renormalised process (n -δ Z [nt] ) t≥0 is shown to converge to a δ-self-similar stationary increments process where δ = 1α -1 + (αβ) -1 .

Since then, many authors considered RWRS models and related functional limit theorems. Bolthausen [START_REF] Bolthausen | A central limit theorem for two-dimensional random walks in random sceneries[END_REF] studied the case of a recurrent Z 2random walk. Maejima [START_REF] Maejima | Limit theorems related to a class of operator-selfsimilar processes[END_REF] considered the case of a random walk in dimension one and of a multivariate random scenery in link with operator self-similar random processes. All these results require independence properties on the increments of the walk and on the scenery. This constraint can be relaxed by assuming some reasonable dependence structure: in Wang [START_REF] Wang | Weak convergence to fractional Brownian motion in Brownian scenery[END_REF], the increments of the random walk are assumed to be strongly correlated, Gaussian and in the domain of attraction of a fractional Brownian motion; it evolves in a i.i.d. square integrable random scenery. In Lang & Xahn [START_REF] Lang | Strongly correlated random fields as observed by a random walker[END_REF] a functional limit theorem is proved for an i.i.d. increments Z-random walk in a strongly correlated scenery satisfying a non central limit theorem of Dobrushin & Major [START_REF] Dobrushin | Non-central limit theorems for non-linear functionals of Gaussian fields[END_REF]. Weakly dependent sceneries have also recently been studied by Guillotin-Plantard & Prieur [START_REF] Guillotin-Plantard | Central limit theorem for sampled sums of dependent random variables[END_REF][START_REF] Guillotin-Plantard | Limit theorem for random walk in weakly dependent random scenery[END_REF] under a θ-mixing condition. Let us also mention the work of Le Borgne [START_REF] Le Borgne | Exemples de systèmes dynamiques quasihyperboliques à décorrélations lentes[END_REF] in which the random scenery verifies some strong decorrelation condition.

One dimensional RWRS arise in the study of random walks evolving on oriented versions of Z 2 (Guillotin-Plantard & Le Ny [START_REF] Guillotin-Plantard | Transient random walks on 2d-oriented lattices[END_REF][START_REF] Guillotin-Plantard | A functional limit theorem for a 2d-random walk with dependent marginals[END_REF]; Pène [START_REF] Pène | Transient random walk in Z 2 with stationary orientations[END_REF]) and in the context of charged polymers (Chen & Khoshnevisan [START_REF] Chen | From charged polymers to random walk in random scenery[END_REF]). Motivated by models for traffic in a network, Cohen & Samorodnitsky [START_REF] Cohen | Random rewards, fractional Brownian local times and stable self-similar processes[END_REF] introduced a random reward schema consisting of sums of independent copies of RWRS and studied its convergence. Extensions of this work have been considered in Guillotin-Plantard & Dombry [START_REF] Dombry | Discrete approximation of a stable self-similar stationary increments process[END_REF] and Cohen & Dombry [START_REF] Cohen | Convergence of dependent walks in a random scenery to fBm-local time fractional stable motions[END_REF].

In this paper, a functional approach for RWRS is developed and the convergence of the finite-dimensional distributions of RWRS under very general assumptions that cover many of the above results are proved. A functional approach using stochastic integration was suggested by Borodin [START_REF] Borodin | Limit theorems for sums of independent random variables defined on a recurrent random walk[END_REF] and Cadre [START_REF] Cadre | Etude de convergence en loi de fonctionnelles de processus: Formes quadratiques ou multilinéaires aléatoires, Temps locaux d'intersection de marches aléatoires[END_REF] in the case when the random scenery is a martingale. Our approach is rather based on the convergence of random measures to random noises (see [START_REF] Dombry | Convergence to stable noise and applications[END_REF]) and is more adapted when the scenery is not a martingale (including weakly and strongly correlated RS). Under the assumption of independence of the RW and of the RS (assumed in all the above cited works), it allows to study separately the convergence of the local time associated with the random walk on the one hand, and the convergence of the random measures associated with the random scenery on the other hand. The technique is also robust enough to deal with a variant of RWRS: the case when many walkers evolve in a single scenery, which might appear as a more realistic model of the random reward scheme considered in [START_REF] Cohen | Random rewards, fractional Brownian local times and stable self-similar processes[END_REF] since all the users share the same network. To our knowledge, these are the first results in that direction.

Our paper is organized as follows. In Section 2, after introducing a general criterion for the L p -convergence of the local time of the random walk, we verify that it applies in the classical cases above-mentioned. In Section 3 a sequence of random measures is associated with the random scenery and conditions to ensure its convergence to a stable, Gaussian or fractional Gaussian random noise are given. In Section 4 we explain how these results allow to recover limit theorems for RWRS and also consider the case of many walkers evolving in one single random scenery.

2 Convergence in L p of the local times of the random walk

A general criterion

Let S = (S n ) n∈N be a Z-valued random sequence. We define the discrete local time at time n ∈ N and point x ∈ Z by

N(n, x) = card{k ∈ [|0, n|]; S k = x}.
Proposition 2.1. Suppose that the following assumption (RW) is satisfied:

• (RW1) There is some sequence a n verifiying a n → ∞ and n -1 a n → 0 such that the renormalized process a -1 n S [nt] t≥0 converges in D([0, ∞)) to some process (Y t ) t≥0 admitting a local time (L(t, x)) t≥0,x∈R .
• (RW2) There is some p ≥ 1 such that for all M > 0:

lim δ→0 lim sup n→∞ [-M,M ] E|L n (t, x) -L n (t, [x] δ )| p dx = 0
where L n (t, x) is the rescaled discrete local time

L n (t, x) = n -1 a n N([nt], [a n x]) (1)
and

[x] δ = δ[δ -1 x].
Then, for every m ≥ 1, for every θ 1 ∈ R, . . . , θ m ∈ R, for every (t 1 , . . . , t m ) such that 0 < t 1 < . . . < t m , the following convergence holds in L p , p ≥ 1:

m i=1 θ i L n (t i , .) L =⇒ m i=1 θ i L(t i , .) as n → ∞.
Remark: A direct application of Lebesgue's dominated convergence Theorem shows that Assumption (RW2) is a consequence of the following two conditions: (RW2.a) There is some A > 0 such that for all n ≥ 1,

sup x∈R E|L n (t, x)| p ≤ A,
(RW2.b) For every t > 0, for almost every x ∈ R,

lim y→0 lim sup n→∞ E|L n (t, x) -L n (t, x + y)| p = 0.
We notice that if these conditions are verified for some p > 1 then they hold for any p ′ ∈ [1, p].

Proof:

We prove the theorem for m = 1 and t 1 = t > 0. Generalizations to any m ≥ 1 can be easily obtained. We denote by q the conjugate of p.

According to Theorem 2.4 of De Acosta (1970) we need to check two things:

-first that there is a w ⋆ -dense set

D ∈ L q (R) such that for each f ∈ D, R f (x)L n (t, x)dx → R f (x)L(t, x)dx,
-second that the sequence (L n (t, .)) n≥1 is flatly concentrated, that is, for every ε > 0, there is a finite dimensional subspace

F of L p (R) such that lim inf n→∞ P [L n (t, .) ∈ F ε ] ≥ 1 -ε where F ε is the ε-neighborhood of F .
The first assertion follows from the invariance principle (RW1) for a -1 n S [n.] . We take D as the set of continuous functions with bounded support and observe that the definition of the discrete occupation time L n and the uniform continuity of f imply

R f (x)L n (t, x)dx = n -1 R f (a -1 n x)N([nt], [x])dx = n -1 R f (a -1 n [x])N([nt], [x])dx + o(1) = n -1 [nt] k=1 f (a -1 n S k ) + o(1) = t 0 f (a -1 n S [nu] )du + o(1)
From the continuity of the application

D([0, t]) → R, φ → t 0 f • φ,
from the invariance principle (RW1) and the existence of the local time L for Y , we deduce

R f (x)L n (t, x)dx → t 0 f (Y u )du = R f (x)L(t, x)dx.
We turn to the second assertion and consider the partition x j = jδ, j ∈ Z for a given δ > 0. Define the finite dimensional space F of L p (R) as the space of functions that are constant on any interval (x j , x j+1 ), j ∈ Z and that vanish outside [-M, M] (for some M > 0 to be specified later). Define Ln by Ln (x) = L n (t, x j ) if x j ≤ x < x j+1 , i.e. Ln (x) = L n (t, [x] δ ). Notice that Ln ∈ F on the event

{L n (t, .) ≡ 0 on R \ [-M, M]} = max 0≤k≤[nt] |a -1 n S k | ≤ M .
Therefore,

P [L n (t, .) / ∈ F ε ] ≤ P max 0≤k≤[nt] |a -1 n S k | > M + P max 0≤k≤[nt] |a -1 n S k | ≤ M and ||L n (t, .) -Ln || L p ≥ ε
By the invariance principle (RW1),

lim n→∞ P max 0≤k≤[nt] |a -1 n S k | > M = P max 0≤s≤t |Y s | > M
is small for large M. Hence it remains only to show that lim

δ→0 lim sup n→∞ P max 0≤k≤[nt] |a -1 n S k | ≤ M and ||L n (t, .) -Ln || L p ≥ ε = 0.
Using Markov's inequality, we get

P max 0≤k≤[nt] |a -1 n S k | ≤ M and ||L n (t, .) -Ln || L p ≥ ε ≤ ε -p E [-M,M ] |L n (t, x) -Ln (t, [x] δ )| p dx
and the result follows from Assumption (RW2).

Local time of a random walk in the domain of attraction of a stable Lévy motion

Let S = (S n ) n∈N be a Z-random walk such that

     S 0 = 0, S n = n k=1 X k , n ≥ 1,
where the X i are i.i.d. integer-valued random variables belonging to the normal domain of attraction of a strictly stable distribution

S α (σ, ν, 0) with parameters α ∈ (1, 2], σ > 0 and ν ∈ [-1, 1].
This means that the following weak convergence holds:

n -1 α S n L =⇒ n→∞ Z α , (2)
where the characteristic function of the random variable Z α is given by

∀u ∈ R, E exp(iuZ α ) = exp -σ α |u| α (1 -iν tan( πα 2 ) sgn(u)) . (3)
Assumption (RW1) holds with a n = n 1/α and (Y (t)) t≥0 an α-stable Lévy process such that Y (0) = 0 and Y (1) is distributed according to Z α . The local time (N(n, x)) n∈N;x∈Z of the random walk S is defined as in the previous section. Its rescaled local time (L n (t, x)) t≥0,x∈R is defined as

L n (t, x) = n -(1-1/α) N([nt], [n 1/α x]).
We denote by (L(t, x)) t≥0,x∈R the jointly continuous version of the local time of the process Y . We are interested in the convergence in L p , p ≥ 1 of the rescaled local time L n to L. By applying the general criterion of Section 2.1, we get the following result.

Proposition 2.2. For every m ≥ 1, for every

θ 1 ∈ R, . . . , θ m ∈ R, for every (t 1 , . . . , t m ) such that 0 < t 1 < . . . < t m , the following convergence holds in L p , p ∈ [1, 2]: m i=1 θ i L n (t i , .) L =⇒ m i=1 θ i L(t i , .) as n → ∞.

Proof:

From the remark following Proposition 2.1, it is enough to prove that conditions (RW2.ab) are satisfied in our setting with p = 2. Condition (RW2.a): Kesten and Spitzer proved in [START_REF] Kesten | A limit theorem related to a new class of self-similar processes[END_REF] ((2.10) in Lemma 1) that if τ x denotes the hitting time of the point x by the random walk S, then the following inequality holds for any x ∈ Z, r ≥ 0 and n ≥ 1,

P(N(n, x) ≥ r) ≤ P(N(n + 1, 0) ≥ r)P(τ x ≤ n + 1).
It implies that the moment of order 2 of the random variable N(n, x) is uniformly bounded by the moment of order 2 of N(n+1, 0) which is equivalent to Cn 2(1-1/α) for some C > 0 (see (2.12) in Lemma 1 of [START_REF] Kesten | A limit theorem related to a new class of self-similar processes[END_REF]).

Condition (RW2.b):

We know (see Lemma 2 and 3 in [START_REF] Kesten | A limit theorem related to a new class of self-similar processes[END_REF]) that there exists some constant C > 0 independent of x, y ∈ Z and n ≥ 1 s.t.

E(|N(n, [n 1/α x])-N(n, [n 1/α (x+y)])| 2 ) ≤ C [n 1/α x] -[n 1/α (x + y)] α-1 n 1-1/α .
Then, we have

lim y→0 lim sup n→∞ E(|L n (t, x) -L n (t, x + y)| 2 ) ≤ C lim y→0 lim sup n→∞ [n 1/α x] n 1/α - [n 1/α (x + y)] n 1/α α-1 = C lim y→0 |y| α-1 = 0 (α > 1 by hypothesis).

Local time of a strongly correlated walk in the domain of attraction of a fractional Brownian motion

Let S = (S n ) n∈N be a Z-random walk such that

     S 0 = 0, S n = n k=1 X k , n ≥ 1,
where (X k ) k≥1 is a stationary Gaussian sequence with mean 0 and correlations

r(i -j) = E(X i X j ) satisfying as n → ∞, n i=1 n j=1 r(i -j) ∼ n 2H l(n)
with 0 < H < 1 and l is a slowly varying function at ∞. Assumption (RW1) is satisfied with a n = n H l(n)

L n,H (t, x) = n -(1-H) l(n) 1/2 N([nt], [n H l(n) 1/2 x]).
We

m i=1 θ i L n,H (t i , .) L =⇒ m i=1 θ i L H (t i , .) as n → ∞.
Proof: Conditions (RW2.ab) directly follow from Lemma 4.4 and Lemma 4.6 of Wang [START_REF] Wang | Weak convergence to fractional Brownian motion in Brownian scenery[END_REF].

The local time of many independent random walks

We consider independent and identically distributed Z-random walks

S (i) = (S (i) n ) n∈N , i ≥ 1,
verifying assumption (RW). The local time of the i-th random walk is denoted by (N (i) (n, x)) n∈N;x∈Z and the rescaled one by (L (i) n (t, x)) t≥0;x∈R . The independent limit processes of the correctly renormalized random walks S (i) have local time denoted by (L (i) (t, x)) t≥0;x∈R . We prove that Proposition 2.4. Let c n → ∞ as n → ∞. Suppose that the random walks S (i) , i ≥ 1 are i.i.d. and satisfy assumption (RW) for some p ≥ 1. Suppose furthermore that for any fixed t ≥ 0,

E R L (1) (t, x) p dx 1/p < ∞. (4)
Then, for every m ≥ 1, for every θ 1 ∈ R, . . . , θ m ∈ R, for every (t 1 , . . . , t m ) such that 0 < t 1 < . . . < t m , the following convergence holds in L p :

1 c n m j=1 θ j cn i=1 L (i) n (t j , .) L =⇒ m j=1 θ j E(L (1) (t j , .)) as n → ∞. (5)
Proof:

From Proposition 2.1 and Skorohod's representation theorem, there exists a probability space on which are defined independent copies of the processes

L (i)
n and L (i) denoted by L(i) n and L(i) such that for every i ≥ 1, the sequence m j=1 θ j L(i) n (t j , .) converges almost surely in L p to m j=1 θ j L(i) (t j , .) as n tends to infinity. So it is enough to prove the result for these later sequences.

Step 1: First, from triangular inequality, we have

1 c n m j=1 θ j cn i=1 L(i) n (t j , .) -L(i) (t j , .) L p (R) ≤ 1 c n cn i=1 || m j=1 θ j L(i) n (t j , .) -L(i) (t j , .) || L p (R) .
and we prove that this quantity converges to zero as n → ∞. To see this, let

X (i) n = || m j=1 θ j L(i) n (t j , .) -L(i) (t j , .) || L p (R) .
The triangular array of random variables {X

(i) n , n ≥ 1, i ∈ {1, . . . , c n }} is such that: -on each row (fixed n), the random variables X (i) n , i ∈ {1, . . . , c n } are i.i.d.. -on each column (fixed i), (X (i)
n ) n≥1 weakly converges to zero as n → ∞. Under these general assumptions, the following weak law of large numbers holds 1 c n cn i=1

X (i) n L =⇒ 0, (6)
implying the required result for the local time.

We finally give some indications for the proof of the weak law of large numbers. Denote by A n the random variable in the left hand side of [START_REF] Chen | From charged polymers to random walk in random scenery[END_REF]. Using the independence of the X (i) n 's, its characteristic function is given by

E [exp(iuA n )] = E exp(iuc -1 n X (1) n ) cn and hence |E [exp(iuA n )] -1| ≤ E min c n | exp(iuc -1 n X (1) n ) -1|, 2 ≤ E min |u|X (1) n , 2 .
The function x → min (|u|x, 2) is continuous and bounded on R + . Then, since the sequence (X

n ) n≥1 weakly converges to 0, we deduce that (E [exp(iuA n )]) n≥1 converges to 1 as n → ∞ and this proves the weak law of large numbers.

Step 2: From Assumption (4) and the strong law of large numbers in the Banach space L p (R) (see [START_REF] Ledoux | Probability in Banach spaces. Isoperimetry and processes[END_REF] for instance), the sequence

1 c n m j=1 θ j cn i=1 L(i) (t j , .)
converges almost surely in L p to the function m j=1 θ j E( L(1) (t j , .)).

Conclusion:

Combining steps 1 and 2 gives [START_REF] Cadre | Etude de convergence en loi de fonctionnelles de processus: Formes quadratiques ou multilinéaires aléatoires, Temps locaux d'intersection de marches aléatoires[END_REF].

Examples:

The above proposition directly applies in the case of random walks with i.i.d. increments in the domain of attraction of a stable Lévy motion (see Section 2.2) or of a strongly correlated Gaussian random walk in the domain of attraction of a fractional Brownian motion (see Section 2.3). In both cases, equation ( 4) is satisfied for every p ∈ [1, 2]: Lemma 2.1 of [START_REF] Dombry | Discrete approximation of a stable self-similar stationary increments process[END_REF] states that the local time of a stable Lévy motion is L p -integrable and Theorem 3.1 of [START_REF] Cohen | Random rewards, fractional Brownian local times and stable self-similar processes[END_REF] gives a similar result for the fractional Brownian motion.

Convergence of the random measures associated with the random scenery

Let ξ = {ξ x , x ∈ Z} be a family of real random variables. In the sequel, we consider µ h = µ h (ξ) the random signed measure on R absolutely continuous with respect to Lebesgue measure with random density

dµ h dx (x) = γ h h -1 k∈Z ξ k 1 [hk,h(k+1)) (x), (7)
where γ h > 0 is a normalisation constant. For a locally integrable function f ∈ L 1 loc , we want to define

µ h [f ] = γ h k∈Z ξ k h -1 h(k+1) hk f (x)dx. ( 8
)
Denote by F µ h the set of functions f ∈ L 1 loc such that this series is convergent in a sense that will be precised later according to the considered cases: either almost-sure semi-convergence or convergence in L 2 (Ω). Clearly, any integrable function with bounded support belongs to the set F µ h . The following scaling relation is worth noting:

µ h [f (c.)] = γ h cγ ch µ ch [f (.)], (9)
whenever these quantities are well-defined.

We introduce the cumulative scenery (w x ) x∈Z

w x =    x-1 i=0 ξ i if x > 0 -1 i=x ξ k if x < 0 0 if x = 0 , x ∈ Z,
and we extend it to a continuous process by the linear interpolation

w x = w [x] + (x -[x])(w [x]+1 -w [x] ).
Note that

µ h [1 [x 1 ,x 2 ] ] = γ h (w h -1 x 2 -w h -1 x 1 ) = W h (x 2 ) -W h (x 1 ), x 1 , x 2 ∈ R, x 1 ≤ x 2 .
(10) with W h the rescaled cumulative scenery W h (x) = γ h w h -1 x .

Independent and identically distributed scenery

Let 1 < β ≤ 2. We consider the case when the ξ x , x ∈ Z satisfy the following assumption: (RS1) the ξ x 's are i.i.d. random variables in the normal domain of attraction of the stable law S β (σ, ν, 0) with characteristic function given in [START_REF] Borodin | Limit theorems for sums of independent random variables defined on a nonrecurrent random walk. Investigations in the theory of probability distributions[END_REF]. In this case, with the scaling γ h = σ -1 h 1 β , the rescaled cumulative scenery W h converges to a bilateral β-stable Lévy process W such that W (1) is distributed according to the stable distribution S β (1, ν, 0). In [START_REF] Dombry | Convergence to stable noise and applications[END_REF], the convergence of the random measures µ h to a stable white noise was proved. Here a function f ∈ L 1 loc belongs to F µ h if the random series (8) defining µ h [f ] converges almost surely. Proposition 3.1. Suppose that ξ satisfies assumption (RS1) for some 1 < β ≤ 2.

Then for any

h > 0, L β (R) ⊂ F µ h and if f n → f in L β and h n → 0 as n → ∞, µ hn [f n ] L =⇒ R f (x)W (dx)
where W denotes the independently scattered β-stable random noise on R with Lebesgue intensity and constant skewness ν.

A general criterion when the scenery is square integrable

We consider the case of a centered scenery with finite variance in the domain of attraction of some square integrable random process W . Our assumption (RS2) is:

• (RS2.a)
The scenery ξ is centered and square integrable.

• (RS2.b) The finite dimensional distributions of the rescaled cumulative scenery W h converge as h → 0:

(W h (x)) x∈R L =⇒ (W (x)) x∈R
where W is a non-degenerate centered square integrable process.

• (RS2.c)

There is some C 1 > 0 such that for any h > 0, K ≥ 1,

x 1 ≤ • • • ≤ x K+1 , θ 1 , • • • , θ K ∈ R, || K i=1 θ i (W h (x i+1 ) -W h (x i ))|| L 2 (Ω) ≤ C 1 || K i=1 θ i (W (x i+1 ) -W (x i ))|| L 2 (Ω) .
Note that W (0) = 0 and that W non degenerate means that for any pairwise distinct nonzero

x 1 , • • • , x p ∈ R, the random variables W (x 1 ), • • • , W (x p ) are linearly independent in L 2 (Ω).
We now define integration of deterministic functions with respect to W . This construction is classical in the case when W is a Gaussian process (e.g. Brownian or fractional Brownian motion), and the reader can refer to [START_REF] Nualart | Stochastic integration with respect to fractional Brownian motion and applications[END_REF], but we simply assume here that W is centered square integrable and non degenerate. Let E be the class of step functions on R, i.e. f ∈ E if and only if f

= K i=1 θ i 1 ]x i ,x i+1 ] for some θ 1 , • • • , θ K and x 1 < • • • < x K+1 . For f ∈ E, we define the random variable W [f ] ∈ L 2 (Ω) by W [f ] = K i=1 θ i (W (x i+1 ) -W (x i ))
and use the notation W [f ] = R f (x)W (dx) to emphasize the analogy with integration (even if the path W need not to have bounded variations on compact sets). Define the scalar product < ., . > W on E by

< f 1 , f 2 > W = Cov(W [f 1 ], W [f 2 ]) , f 1 , f 2 ∈ E.
This is indeed a scalar product since the process W is non degenerate. With these notations, W :

E → L 2 (Ω), f → W [f ] is a linear isometry. Define the Hilbert space L 2
W as the closure of E with respect to the scalar product < ., . > W . Since L 2 (Ω) is complete and E is dense in L 2 W (by construction), we can extend the isometry into W :

L 2 W → L 2 (Ω).

Examples:

If W is the Brownian motion, the scalar product is given by

< f 1 , f 2 > W = R f 1 (u)f 2 (u)du, f 1 , f 2 ∈ E and then L 2 W = L 2 (R). Then W [f ]
is the Wiener integral of f with respect to Gaussian white noise. If W is the fractional Brownian motion of index H ∈ (1/2, 1), then

< f 1 , f 2 > W = H(2H -1) R×R |u -v| 2H-2 f 1 (u)f 2 (v)dudv, f 1 , f 2 ∈ E
and elements of L 2 W may not be functions but rather distributions of negative order (see [START_REF] Nualart | The Malliavin calculus and related topics[END_REF] or [START_REF] Pipiras | Integration questions related to fractional Brownian motion[END_REF]). To our purpose, we will not have to consider this distributions space, but it will be enough to remark that

L 1 ∩ L 2 ⊂ L 2 W and that ||f || L 2 W ≤ C 2 ||f || L 1 ∩L 2 for some constant C 2 , where L 1 ∩ L 2 is equipped with the norm ||f || L 1 ∩L 2 = max(||f || L 1 , ||f || L 2 ).
In order to compare integration with respect to µ h and to W , we have to suppose that L 2 W contains some suitable subspace of locally integrable functions since F µ h ⊆ L 1 loc is the set of functions f ∈ L 1 loc such that the random series (8) defining µ h [f ] converges in L 2 (Ω). In view of the above examples, we make the following further assumption:

• (RS2.d) L 1 ∩ L 2 ⊆ L 2
W with continuous injection, i.e. there is some

C 2 > 0 such that for any f ∈ L 1 ∩ L 2 , ||f || L 2 W ≤ C 2 ||f || L 1 ∩L 2 .
We then prove the convergence of the random measures µ h defined by ( 7) to the random noise W on the space

L 1 ∩ L 2 ⊆ L 2 W .
Proposition 3.2. Suppose that ξ satisfies assumptions (RS2).

Then for any h > 0,

L 1 ∩ L 2 ⊂ F µ h and if f n → f in L 1 ∩ L 2 and h n → 0 as n → ∞, µ hn [f n ] L =⇒ W [f ] = R f (x)W (dx).
Remark: It is worth noting that in the important case when W is a Brownian motion, assumption (RS2.d) is trivial. Proposition 3.2 applies under assumptions (RS2.bc) and the conclusion can be strengthened: the convergence holds on the whole space L 2 = L 2 W = F µ h rather than only on the subspace L 1 ∩ L 2 .

Proof:

For the sake of clarity, we divide the proof into four steps.

Step 1: First note that in the simple case when

f n ≡ f ∈ E, Proposition 3.2 directly follows from assumption (RS2.b). Let f = K i=1 θ i 1 ]x i ,x i+1 ]
. Using equation [START_REF] Dobrushin | Non-central limit theorems for non-linear functionals of Gaussian fields[END_REF] ,

µ hn [f ] = K i=1 θ i (W hn (x i+1 ) -W hn (x i )).
Then, assumption (RS2.b) states that as n → +∞,

µ hn [f ] L =⇒ K i=1 θ i (W (x i+1 ) -W (x i )) = W [f ].
Step 2: We now prove that for every fixed integer n, L 1 ∩ L 2 ⊂ F µ hn and that µ hn induces a linear functional

L 1 ∩ L 2 → L 2 (Ω) with norm operator bounded by C 1 C 2 .
The linearity of µ hn on E is straightforward. Using equation [START_REF] Dobrushin | Non-central limit theorems for non-linear functionals of Gaussian fields[END_REF], assumption (RS2.c) can be rewritten as

∀f ∈ E, ||µ hn [f ]|| L 2 (Ω) ≤ C 1 ||W [f ]|| L 2 (Ω) = C 1 ||f || L 2
W , and then, in view of assumption (RS2.d), for every f ∈ E,

||µ hn [f ]|| L 2 (Ω) ≤ C 1 C 2 ||f || L 1 ∩L 2 .
Since E is dense in L 1 ∩ L 2 and L 2 (Ω) is complete, the linear functional µ hn can be extended to L 1 ∩ L 2 and its norm operator is bounded by C 1 C 2 .

Step 3: We consider now the case when

f n ≡ f belongs to L 1 ∩ L 2 .
We approximate f by a sequence of simple functions

f m such that f m → f in L 1 ∩ L 2 as m → ∞. Let Y n = µ hn [f ], X m,n = µ hn [f m ], X m = W [f m ], X = W [f ].
We apply Theorem 4.2 in [START_REF] Billingsley | Convergence of probability measures[END_REF] by checking that : -for fixed m, X m,n weakly converges to X m (use step 1 with the simple function f m ), -X m weakly converges to X because using assumption (RS2.d),

||X m -X|| L 2 (Ω) = ||f m -f || L 2 W ≤ ||f m -f || L 1 ∩L 2 → 0, -lim m∞ lim sup n∞ P (|X m,n -Y n | ≥ ε) = 0 for any ε > 0.
The last condition is satisfied since

P (|X m,n -Y n | ≥ ε) = P (|µ hn [f m -f ]| ≥ ε) ≤ ε -2 ||µ hn [f m -f ]|| 2 L 2 (Ω) ≤ (C 1 C 2 ) 2 ε -2 ||f m -f || 2 L 1 ∩L 2 .
Using Theorem 4.2 in [START_REF] Billingsley | Convergence of probability measures[END_REF], we deduce that Y n weakly converges to X as n tends to infinity.

Step 4: We consider the general case f n → f in L 1 ∩ L 2 , and we write

µ hn [f n ] = µ hn [f ] + µ hn [f n -f ]. From step 3, µ hn [f ] weakly converges to W [f ]. From step 2, ||µ hn [f n -f ]|| L 2 (Ω) ≤ C 1 C 2 ||f n -f || L 1 ∩L 2 and hence µ hn [f n -f ] converges to 0 in L 2 (Ω). As a consequence, µ hn [f n ] weakly converges to W [f ].
Note that the above proposition can be applied in the following cases.

Example: Weakly dependent sceneries.

Under the assumption of a L 2 stationary scenery satisfying some θ-mixing condition, a central limit theorem for triangular arrays implying assumption (RS2.b) was proved in [START_REF] Guillotin-Plantard | Central limit theorem for sampled sums of dependent random variables[END_REF] (Theorem 3.1) as well as a suitable control of the covariance implying assumption (RS2.c) (see Lemma 7.1 in [START_REF] Guillotin-Plantard | Central limit theorem for sampled sums of dependent random variables[END_REF]). By applying the techniques developed in [START_REF] Le Borgne | Exemples de systèmes dynamiques quasihyperboliques à décorrélations lentes[END_REF], a Kesten-Spitzer result for a simple Z-random walk in a random scenery coming from a partially hyperbolic dynamical system can be established. It holds when the scenery satisfies some strong decorrelation properties given in Theorem 1 of [START_REF] Pène | Transient random walk in Z 2 with stationary orientations[END_REF]. We just mention that conditions (RS2.ab) should hold under these last assumptions.

Example: Moving averages sceneries.

Stationary sceneries obtained from moving averages of i.i.d. square integrable random variables are considered in [START_REF] Dombry | Convergence to stable noise and applications[END_REF], i.e.

ξ k = i∈Z c k-i η i
where (η i ) i∈Z is an i.i.d. square integrable sequence. In the case when c is summable, then the convergence of the rescaled cumulative random scenery to a Brownian motion is shown and assumptions (RS2.bc) hold.

In the case when c slowly decays, strong correlations persist in the limit. More precisely, if

lim k→+∞ |k| -γ c k = p 1 and lim k→-∞ |k| -γ c k = p 2
with γ = (1/2, 1) and p 1 p 2 = 0, then convergence of the rescaled scenery to a fractional Brownian motion of index H = 3/2γ ∈ (1/2, 1) is shown and assumptions (RS2.bc) are satisfied. Note that in [START_REF] Dombry | Convergence to stable noise and applications[END_REF], the dimension d of the scenery needs not to be equal to one (spatial sceneries) and the innovations η i are more generally assumed to belong to the domain of attraction of a β-stable distribution with 1 < β ≤ 2. In this framework, the convergence of the random measures µ h to a fractional β-stable random noise on the space L β is shown. The method of convergence of random measures to a random noise is thus very robust. We will show how this can be used to obtain new results on RWRS.

Limit theorems for RWRS

The above results can be applied to prove the convergence of RWRS and to obtain limit theorems under fairly general assumptions. The main idea is that the RWRS can be rewritten as a function of the rescaled local time L n (., .) of the random walk S and of the random measures µ h associated with the random scenery ξ. Using the scaling property (9) and introducing the suitable scalings for the local time (equation ( 1)) and for the random measures (equation ( 7)), we get

Z [nt] = [nt] k=0 ξ S k = x∈Z N([nt], x)ξ x = γ -1 1 R N([nt], [x])µ 1 (dx) = n γ a -1 n µ a -1 n [L n (t, .)] Or equivalently, n -1 γ a -1 n Z [nt] = µ a -1 n [L n (t, .)] (11)
with a -1 n → 0. We then use the convergence of the local time and of the random measure to prove the convergence of the RWRS. Two cases are considered and stated in the same following theorem. Theorem 4.1. Suppose that the random walk S and the random scenery ξ are independent and verify either:

1. (First case) ξ satisfies (RS1) for some β ∈ (1, 2] and S satisfies (RW) with p = β.

2. (Second case) ξ satisfies (RS2) and S satisfies (RW) with p = 2.

Then, the RWRS satisfy the following convergence in the sense of the finite dimensional distributions:

(n -1 γ a -1 n Z [nt] ) t≥0 L =⇒ R L(t, x)W (dx) t≥0
where (W (x)) x∈R and (L(t, x)) t≥0;x∈R are defined in (RS) and (RW) respectively and independent of each other.

Proof:

We only prove the theorem in the first case. Thanks to Proposition 2.2 and Skorohod's representation theorem, there exists one probability space on which are defined copies of the processes L n and L denoted by Ln and L such that the sequence ( Ln ) n converges almost surely in L p to L as n tends to infinity. Then, by combining equality [START_REF] Dombry | Convergence to stable noise and applications[END_REF] and Proposition 3.1, we get that conditionally on these local times, for every m ≥ 1, for every

θ 1 ∈ R, . . . , θ m ∈ R, for every 0 ≤ t 1 < . . . < t m , n -1 γ a -1 n m i=1 θ i Z [nt i ] = m i=1 θ i µ a -1 n [ Ln (t i , .)] (12)
converges in distribution to

m i=1 θ i R L(t i , x) W (dx).
The result then easily follows using dominated convergence theorem.

We finally state a result for many walkers in random scenery. We consider c n independent and identically distributed random walks S (i) , i = 1, . . . , c n defined either as in Section 2.2 or as in Section 2.3. They are assumed to evolve in the same random scenery ξ. The i-th random walk induces a RWRS denoted by Z (i) , but the processes Z (1) , • • • , Z (cn) are not independent since all random walkers evolve in the same scenery (see [START_REF] Cohen | Random rewards, fractional Brownian local times and stable self-similar processes[END_REF][START_REF] Cohen | Convergence of dependent walks in a random scenery to fBm-local time fractional stable motions[END_REF][START_REF] Dombry | Discrete approximation of a stable self-similar stationary increments process[END_REF]). The total random reward at time n is then given by Z n,cn = cn i=1 Z (i) n .

Using equation [START_REF] Dombry | Convergence to stable noise and applications[END_REF], we get

c -1 n n -1 γ a -1 n Z [nt],cn = µ a -1 n [c -1 n cn i=1 L (i) n (t, .)]. (13)
For a finite number c n ≡ c of independent walkers in a random scenery the functional approach easily gives us the following convergence

(n -1 γ a -1 n c i=1 Z (i) [nt] ) t≥0 L =⇒ R ( c i=1 L (i) (t, x))W (dx) t≥0 .
When c n → ∞ as n → ∞, thanks to Propositions 2.4 and 3.1, we get the following result, whose proof is very similar to the proof of Theorem 4.1 and is omitted.

Theorem 4.2. Suppose that the random walk S (i) and the random scenery ξ are independent and verify either:

1. (First case) ξ satisfies (RS1) for some β ∈ (1, 2] and the S (i) verify (RW) as well as the integrability condition ( 4) with p = β.

2. (Second case) ξ satisfies (RS2) and the S (i) verify (RW) as well as the integrability condition ( 4) with p = 2.

Then, the total random reward satisfies the following convergence in the sense of the finite dimensional distributions:

(c -1 n n -1 γ a -1 n Z [nt],cn ) t≥0 L =⇒ R E(L(t, x))W (dx) t≥0 .
where (W (x)) x∈R and (L(t, x)) t≥0;x∈R are defined in (RS) and (RW) respectively.

Conclusion and further comments

In this paper a functional approach for the analysis of RWRS was proposed. The method is quite powerful and allows to recover many of the previous results about RWRS as well as new ones. However, our results are limited to the convergence of the finite-dimensional distributions and the question of the functional convergence (for example in the Skohorod topology) naturally arises. The classical techniques used for proving tightness such as Kolmogorov's criterion can be adapted. The reader should refer to [START_REF] Cohen | Convergence of dependent walks in a random scenery to fBm-local time fractional stable motions[END_REF], [START_REF] Dombry | Discrete approximation of a stable self-similar stationary increments process[END_REF], [START_REF] Kesten | A limit theorem related to a new class of self-similar processes[END_REF] or [START_REF] Wang | Weak convergence to fractional Brownian motion in Brownian scenery[END_REF] where tightness results were proved. In a few words, in the case of a random scenery satisfying (RS2), the second order moment of the increments of the RWRS is given for any 0 ≤ t 1 < t 2 by

E |n -1 γ a -1 n Z [nt 2 ] -n -1 γ a -1 n Z [nt 1 ] | 2 = E |µ a -1 n [L n (t 2 , .) -L n (t 1 , .)]| 2 ≤ (C 1 C 2 ) 2 E R |L n (t 2 , x) -L n (t 1 , x)| 2 dx .
The proof of the tightness then relies on the analysis of the self-intersection local time of the random walk (see [START_REF] Kesten | A limit theorem related to a new class of self-similar processes[END_REF] for RW with i.i.d. increments or [START_REF] Wang | Weak convergence to fractional Brownian motion in Brownian scenery[END_REF] for strongly correlated Gaussian RW).

  1/2 and Y the standard fractional Brownian motion (B H (t)) t≥0 of index 0 < H < 1, i.e. the centered Gaussian process such that B H (0) = 0 and E(|BH (t 2 ) -B H (t 1 )| 2 ) = |t 2t 1 | 2H . The local time N H (n, x)with n ∈ N and x ∈ Z is defined as the number of visits of (S n ) n∈N at point x up to time n. Its rescaled local time (L n,H (t, x)) t≥0,x∈R is defined as

  denote by (L H (t, x)) t≥0,x∈R the jointly continuous version of the local time of the process B H . We are interested in the convergence in L p , p ≥ 1 of the rescaled local time L n,H to L H . By applying the general criterion of Section 2.1, we get the following result.

Proposition 2.3. For every m ≥ 1, for every θ 1 ∈ R, . . . , θ m ∈ R, for every (t 1 , . . . , t m ) such that 0 < t 1 < . . . < t m , the following convergence holds in L p , p ∈ [1, 2]:


Claude Jauffret
email: jauffret@univ-tln.fr

Annie-Claude Pérez
email: annie-claude.perez@univ-tln.fr

Pierre Blanc-Benon
email: pierre.blanc-benon@fr.thalesgroup.com

Hervé Tanguy
email: herve.tanguy@fr.thalesgroup.com

Doppler-Only Target Motion Analysis in a High Duty Cycle Sonar System

Keywords: high duty cycle, bistatic sonar, multi-static sonar, Doppler-only TMA, estimation, Fisher Information Matrix, Cramér-Rao lower bound

This paper presents a new method, in an underwater context, of estimating the trajectory of a target assumed to be in constant velocity motion. The observer system is composed of a transmitter (emitting continuously, said also in high duty cycle mode) and a set of passive underwater acoustic sensors (e.g. receiving sonobuoys) not necessarily collocated with the transmitter. Each receiver detects the single-tone signal emitted by the transmitter after its reflection on the target. The two sets of doubly Dopplershifted frequency measurements for each receiver make the target motion analysis (TMA) feasible, even when only one receiver is available. First, we develop the exact model of the bistatic received frequencies and then we compute the Cramér-Rao lower bound (CRLB) when the number of receivers is 1, 2, or 3, the measurements being corrupted by

an additive and time-independent Gaussian noise. Finally, using some realistic scenarios, we show that the maximum likelihood estimated is actually efficient.

I. INTRODUCTION

Continuous transmission was first in use in the 1930s and 1940s for radar applications and electromagnetic (EM) barriers 1 . In 1984, Gough et al. [START_REF] Gough | Continuous Transmission FM Sonar with One Octave Bandwidth and no Blind Time[END_REF] first proposed continuous transmission frequency modulated (CTFM) sonar in order to increase the signal-to-noise ratio and to provide continuous target location updates. The principle of bistatic sonar was pioneered by Henry Cox in his famous 1988 NATO ASI paper [START_REF] Cox | Fundamentals of Bistatic Active Sonar[END_REF]. Cox poses the "fundamentals of bistatic active sonar", introduced by the first sentence of his paper, which summarizes perfectly the context: "the bistatic situation is characterized by the triangle of source, target and receiver positions, and by their velocities". The so-called high duty cycle (HDC) sonar, named also continuous active sonar (CAS), is a logical combination of both aspects. Since the publication of this article by H. Cox 27 years ago, numerous papers have become available in the open literature, dealing with several questions about HDC sonar: How many transmitters and receivers must be used to obtain a certain performance of target detection and target position estimation? How can the positions of drifting receivers and transmitters be estimated? How can ad hoc signal processing be realized? HDC sonar was studied in the late 1990s following radar applications, e.g. CTFM radar.

Among the numerous variations of HDC concepts, DeFerrari studied continuous m-codes in order to extend the temporal coherent integration time compared to single pulse sonar [START_REF] Deferrari | The Application of m-Sequences to Bi-static Active Sonar[END_REF][START_REF] Deferrari | Continuous active pulse compression sonar[END_REF]. Yang [START_REF] Yang | Acoustic Dopplergram for Intruder Defense[END_REF] proposed a continuous tracking for underwater surveillance from what he called Dopplergram, using continuous wave (CW) signals. Blanc-Benon et al. [START_REF] Blanc-Benon | CAS TMA Applied for Sonobuoys in MPA Context[END_REF] studied CAS-TMA Cramér-Rao lower bound (CRLB) in a maritime patrol aircraft context with Doppler and bearing measurements from acoustic sonobuoys. In [START_REF] Liang | On Designing the Transmission and Reception of Multistatic Continuous Active Sonar[END_REF], Liang et al. focused on the signal processing aspect of multi-static continuous active sonar.

In this paper, we study the effect of the propagation delay on Doppler-only TMA in a HDC sonar system. The target has a constant velocity motion. We compare the behavior of the maximum likelihood estimator (MLE) when this delay is taken into account and when it is not. To do this, we consider one transmitter and one or several receiver(s) for various typical scenarios. In each of them, we derive the CRLB, the empirical bias and the empirical covariance matrix of the MLE via numerous Monte-Carlo simulations.

The paper is organized as follows: In Section II, the exact model of a Doppler-shifted frequency is rigorously established and the common approximation in use with the monostatic sonar domain is derived. In Section III, the exact model and an approximated model are compared in the bistatic situation. Section IV is devoted to the Doppleronly TMA estimation in a bistatic situation. In particular, we illustrate by simulations that using the approximated model yields a biased MLE, whereas the exact MLE is efficient in practice. Section V presents an extension of the proposed method to multi-static situations.

The conclusion addresses perspectives on HDC/CAS sonar and future enhancements regarding the sonobuoy drift problem and bearing contribution.

II. THE MODELS OF THE DOPPLER EFFECT

Suppose that a source S is radiating continuously a signal

  t s E . At time t , a receiver R detects the signal   t s R defined by       t t s t s E R   
. The delay   t  is the propagation time taken by the wave carrying the signal to reach the receiver. If

the signal   t s E is a single tone, that is         t f a t s E 0 2 sin , the received signal is               t t f a t s R 0 2 sin ' . It follows that the instantaneous frequency at time t is       t f t f    1 0 .
(1) This is the exact expression of the Doppler-shifted received frequency. Consequently, the time derivative of   t  is needed. The delay is recursively defined by

        c t P t t P t R S      , ( 2
)
where c is the sound wave speed in the medium,  

            c t r f t f  1 0 . (3)
When the source has a constant velocity S V , the time delay

  t  is proven to be (from (2))             2 2 2 2 2 2 S RS S RS T S RS T S V c t P V c t P V t P V t        , (4)
with

      t P t P t P R S RS  
. Consequently, the time derivative of   t  depends also on the motion of the receiver. Because we are concerned in the multi-static configuration by stationary receivers only, we will compute   t  for a motionless receiver. Under this assumption, from (4) we get

              . 2 2 2 2 2 2 2 2 2 2 t P V c t P V V c t P V c V c V t RS S RS T S S RS T S S S         or equivalently       2 2 2 2 2 1 sin 1 1 1 1                              c v c v t c v c v c t r t S S S S     , (5)
where

      t P V t RS S ,    , S S V v 
, and

      t v t V t r S S   cos cos   
(see [START_REF] Jauffret | Target Motion Analysis by Inverse Triangulation[END_REF]). As a consequence, we can propose two approximations of   t  : First order approximation (from (5)): [START_REF] Blanc-Benon | CAS TMA Applied for Sonobuoys in MPA Context[END_REF], we get the approximation

      2 2 2 2 2 1 sin 2 1 1 1                                              c v c v t c v c v c t r t S S S S     . (6) Second approximation ( c v S  ): If we neglect 2       c v S in
    c t r t     . ( 7
)
With this second approximation, we end up with the common expression of the Doppler-shifted frequency:

            c t r f t f  1 0 .
In conclusion, when the source is in CV motion, and the receiver is motionless, the common expression of the Doppler-shifted frequency is the result of two successive approximations.

Note that, in this case

,       c t t r t     .

III. THE DOPPLER EFFECT IN A BISTATIC CONFIGURATION

In a bistatic configuration, a target T, an emitter (or transmitter) E, and a receiver R are assumed to lie on the same horizontal plane. The transmitter sends a signal continuously and the target plays the role of mirror: it reflects the signal emitted by the transmitter toward the receiver (see Fig. 1). Unlike a classic active system, the receiver stands apart from the transmitter. So, the receiver detects the signal

  t s R defined by       t t s t s R 2 Reflected 2     ,
      u u s u s E 1 1 Reflected    
, for any u , 1  accounts for the signal attenuation. In our configuration, E and R are motionless.

      c t t r t 2 2 2    
. The chain rule allows the signal detected by the receiver to be defined as

          t t t t s t s R 2 1 2 E 2 1         
for any t . Again, when E s is a single tone with the frequency 0 f , the instantaneous received frequency is

              t t t t f t f 2 1 2 2 0 1 1              . ( 8
)
Once more, several approximations of   t f can be proposed: The first approximation is based upon the first-order approximation

      t t t 1 2 1        . We get                       t f t t f t t t f t f 1 2 1 0 1 2 2 0 1 1 1 1                   (9).
Another approximation consists of using the approximation

(7) for   t 2  , that is     c t r t 2 2     . We then get         t f c t r c t r f t f 2 2 1 0 1 1                   . ( 10
)
A third approximation is obtained by neglecting the term

    2 2 1 c t r t r  
; we end up with the common approximation used in the bistatic radar systems:

        t f c t r t r f t f 3 2 1 0 1             . ( 11
  t f 3 in magenta for 1300 0   t s.
We also simulated the temporal received signal   t s R , sampled at 3000 Hz, and performed an 8192-point FFT after zero-padding. The outputs of the FFT are plotted in red in Figs. 2 (a) and 2 (b), which is a magnification of Fig. 2 (a). This simulation shows that the exact model fits the frequency track produced by the spectral analysis. It also exhibits the approximation error in the three approximated models. In this example, one can see the double closest-point-of-approach (CPA) inflexion points: the first one regarding the transmitter, the second one regarding the receiver. given by the output of the FFT. Note that in a radar bistatic system, the approximation ( 11) is acceptable thanks to the EM waves' speed (which is huge relative to the aircraft speed) unlike in the sonar domain.

IV. TMA IN A BISTATIC SYSTEM

A. Problem setting

The common configuration of a simple bistatic system consists of a transmitter E and a receiver R, both mounted on sonobuoys. In this study, they are assumed to be stationary. A target T is in CV motion and is detected only by the receiver from the signal coming from T. This signal is actually emitted by E and reflected by T toward R. The target velocity is denoted as T V , and its position at the very beginning of the scenario is   0 T P . The transmitter and the receiver are located at E P and R P , respectively. After an ad hoc spectral analysis, only the set of Doppler-shifted frequency measurements is available:

    k k k m t f t f    , for K k , , 1   and
k  the corrupting noise, which means that the sonobuoy carrying the receiver is not a directive sonobuoy [START_REF] Maranda | The Statistical Accuracy of an Arctangent Bearing Estimator[END_REF] but rather an omnidirectional buoy. Implicitly,  

k t f is completely defined by     T T T T T V P X 0  . Our goal is to estimate X , given   K k t f k m , , 1 ,   .
To do this, we have two options: (i) use the correct model given by [START_REF] Jauffret | Target Motion Analysis by Inverse Triangulation[END_REF]; or (ii) use an approximate model. We choose to put in competition the correct model and the approximated model [START_REF] Xiao | Observability and Performance Analysis of Bi/Multi-static Doppler-Only Radar[END_REF]. Indeed, the correct model is expected to give better results in the estimation process, but the computational burden will definitively be heavy compared to an estimation based upon the model [START_REF] Xiao | Observability and Performance Analysis of Bi/Multi-static Doppler-Only Radar[END_REF].

Approximated mismatch : 3 Hz

B. Observability

For the approximated model [START_REF] Xiao | Observability and Performance Analysis of Bi/Multi-static Doppler-Only Radar[END_REF], X was proven to be observable in [START_REF] Xiao | Observability and Performance Analysis of Bi/Multi-static Doppler-Only Radar[END_REF] by using the rank of the Fisher information matrix (FIM). Obviously, the state vector X is not observable since X is undistinguishable from its image by the axial symmetry around the line (ER) and from its image by the axial symmetry around the perpendicular bisector of the segment [ER]. The line [ER] and the bisector define 4 sectors of the plane. We claim that, in each sector, the condition of observability is that the target must not travel on the line defined by the couple (E, R). This will be the topic of a future paper. For the correct model ( 8), observability has not been proven so far. To bypass this lack in the rest of the paper, we will consider scenarios for which the FIM is non-singular, so X is locally observable [START_REF] Jauffret | Observability and Fisher Information Matrix in Nonlinear Regression[END_REF].

C. Fisher Information Matrix and Cramér-Rao lower bound

The classic expression of the FIM under Gaussian hypothesis

is            K k k T X k X t f t f X F 1 2 1  , with                   . , 1 , , , , 1 , 2 2 1 2 1 2 0 X t X X t t X X t t X t f t f k X k k k k X k k X                   
We have to compute the following four quantities:  

X t k , 1  ,   X t k , 2  ,   X t k , 1  , and   X t k , 2  .
We end up with

      c X t r t P c X t k k k , 1 , 1 1 1    ,             2 2 2 2 2 2 2 2 2 2 , T k T k T T k T T k V c t P V c t P V t P V X t        ,         c X t r t P c t P V X t k T T k , , 1 1 1 1      ,               . , 2 2 2 2 2 2 2 2 2 2 2 2 2 k T k T T T k T T T T k t P V c t P V V c t P V c V c V X t       



The detailed computations are presented in the Appendix.

D. Estimation

We have chosen the MLE, which is also the least squares

estimator X ˆ, that minimizes             K k k k m t f t f X C 1 2 2 1  for the correct model (8), or the approximated MLE app X ˆ, that minimizes             K k k k m app t f t f X C 1 2 3 2 1 
for the approximated model [START_REF] Xiao | Observability and Performance Analysis of Bi/Multi-static Doppler-Only Radar[END_REF]. Both are obtained by a numerical routine (Gauss-Newton). 3) Scenario 3 : These simulations reveal two facts: (i) the bias induced by the approximated model is not negligible; and (ii) the MLE performs properly for the exact model (its sample covariance matrixnot reported hereis very close to the CRLB and its bias is relatively small). A legitimate question must be asked: Is this bias persistent if more than one receiver is employed?

E. Monte-Carlo simulations

3 X X 
We address this question in the next section.

V. TMA IN MULTI-STATIC CONFIGURATION

When ) 1 (  N receivers are available, we can exploit the frequency track collected by each receiver to improve the TMA. The criterion is then

                 N n K k k n k m n n t f t f X C 1 1 2 2 1  (12)
where

    k m n t f
is the frequency measured by the receiver # n at time k t , n K is the number of measurements acquired by the receiver # n , and

    k n t f is the exact model. Consequently, the FIM is                 N n K k k n T X k n X n t f t f X F 1 1 2 1  . Note that,
for the sake of simplicity in the writing of the equations, we have considered that the measurements are synchronous. But we can readily consider TMA when this is not the case: the only thing that matters is that the measurements need to be time stamped. In the following subsections, we compare the respective performances of X ˆ, that minimizes [START_REF]OCEANS'15, Continuous Active Sonar session, MTS/IEEE Proceedings[END_REF], and of app X ˆ, that minimizes [START_REF]OCEANS'15, Continuous Active Sonar session, MTS/IEEE Proceedings[END_REF], in which we replace the correct model by the approximated model ( 11). This comparison is based upon scenario 1.

A. With two receivers Since we benefit from a second set of measurements, the duration of scenario 1 is reduced to 10 min. Fig. 6 shows that the bias is maintained when the approximated model is used. So, again we increase the number of receivers in the next subsection.

B. With three receivers With three receivers, the duration of scenario 1 is now set to 5 min. In Fig. 7, a non-negligible bias, whose components are given in Table V, is visible when the approximated model is used. We note a reduction of the bias with an increasing number of receivers. However, the lesson to be drawn from these simulations is clear: the approximated model provides a biased estimator. For the HDC sonar, we hence recommend to use the correct model to construct the MLE, even if it requires more computational time. The good news is that the MLE estimate for an HDC sonar system has a covariance matrix compatible with an operational exploitation for short favorable scenarios (about one minute), which is an advantage in real situations where the sonobuoys carrying the transmitter and the receivers are set to drift.

VI. CONCLUSION AND PERSPECTIVES Frequency-only TMA in a multi-static sonar configuration has been studied in this paper. Two major results were obtained: first, the approximated model, widely used in multi-static radar systems, induces a non-negligible biased estimator in sonar applications. On the contrary, the performance of MLE based upon the exact model is compatible with the asymptotic performance given by the CRLB: its empirical covariance matrix is close to the bound and it is practically unbiased. The price is an augmentation of the computation time. Finally, in short favorable scenarios, the CRLB allows us to expect TMA results that are compatible with the operational exploitation of a multi-static system (provided the number of receivers is sufficient). As HDC sonar is becoming more and more popular [START_REF]OCEANS'15, Continuous Active Sonar session, MTS/IEEE Proceedings[END_REF], further studies should analyze the impact of drifting of receivers due to ocean currents and in GPS-denied areas. The contribution of bearing measurement should also be analyzed when using directive receivers in terms of duration to get a given target location accuracy compared to the case of frequency-only TMA.

APPENDIX

For convenience, hereafter we will use the notation

  T T T T X V V 0 0 2 2   , 2 2 2 2 2 2 2 T X T T T X V V c c V c V              .                  T k k i T k i T T X V t t P V t P V ,           k i T T X k i T T k i T T X t P V t P V t P V    2 2 .                 k i k k i k i X t P t t P t P 2 2 .             2 2 2 2 2 2 2 2 T X k i k i X T T k i X V t P t P V c V c t P          . Let us define         2 2 2 2 ) ( k i T k i T T k i t P V c t P V t a    .           2 2 2 2 ) ( 2 1 ) ( k i X T k i T T X k i k i X t P V c t P V t a t a        .       2 2 2 2 2 ) ( ) ( ) ( T X k i k i X T k i T X V t a t a V c t a V c                                ) ( 2 2 2 2 2 2 2 2 2 2 2 2 2 k i T X k i T k T T k i T k T T X k i T k T T X t a V c t a V c t P V t a V c t P V t a V c t P V                 
Now, we are able to give the values of the first two gradients:

          . ) ( ) ( , 2 2 2 2 2 2 2 2 2 2 2 T X T k T T k T k T T X k X k X V V c t P V t a V c t P V t a X t                                               k T k T T X T T X k X t a V c t P V c V c V t 2 2 2 2 2 2 2 2 2  . Computation of     X X t t k k X , , 2 1      : First, we write     X X t t k k , , 2 1

  



as the composition of the following two functions:

         . , , , 2 1 1 2 5 4 1 X X t t X g h Y X u X X t t Y X k k h k k g h g                                     R R R
We readily deduce that its gradient is equal to the product of two gradients, that is

        X g h X X t t X k k X       , , 2 1                  X X t t Y Y X k k h g , 2 .  with       4 2 2 , , I X t X X t t g k X k k X X                  and                          h u X u h h X Y ,
. So we get:

                            h u X u h X t X X t t X k X k k X , , , , 4 2 2 1 I          X u h X t u X u h X k X , , , 2         for   X t t u k k , 2              X u X t u X u X X t t X k X k k X , , , , , 1 2 1 2 1                  , for   X t t u k k , 2    .

  where Reflected s is the signal reflected by T, 2  is the propagation time of the targetreceiver range, and 2  accounts for the attenuation due to propagation. The signal Reflected s is the signal emitted by E,

Fig. 1 :

 1 Fig. 1: Standard scenario of the bistatic HDC sonar In the couple (E,T), the source is E; hence     c t r t 1 1  

  , R), T plays the role of the source; therefore,

  Fig. 2(a) and 2(b) depict   t f in blue,   t f 1 in green,   t f 2 in cyan, and

Fig. 2 (

 2 Fig. 2 (a): Frequency tracks with double inflexion CPA points (x axis: time in seconds, y axis: frequency in Hz).

Fig. 2 (

 2 Fig. 2 (b): Magnification of the frequency tracks; the dottedline is given by the output of the FFT. Note that in a radar bistatic system, the approximation (11) is acceptable thanks to the EM waves' speed (which is huge relative to the aircraft speed) unlike in the sonar domain.

Fig. 4 (

 4 Fig. 4 (a): Target's trajectory, initial positions estimated with exact model (blue dots) and with approximated model (red crosses), and the 90% confidence ellipse (magenta).

Fig. 4 (

 4 Fig. 4 (b): A magnification of the clouds of estimates of Fig. 4(a).

Fig. 5 (

 5 Fig. 5 (b): A magnification of the clouds of estimates of Fig. 5(a).

Fig. 6 (

 6 Fig. 6 (a): Target's trajectory and initial positions estimated with exact model (blue dots) and with approximated model (red crosses) in a tri-static configuration.

Fig. 6 (

 6 Fig. 6 (b): A magnification of the clouds of estimates of Fig. 6(a).

Fig. 7 (

 7 Fig. 7 (b): A magnification of the clouds of estimates of Fig.7(a). In Fig.7, a non-negligible bias, whose components are given in TableV, is visible when the approximated model is used. We note a reduction of the bias with an increasing number of receivers. However, the lesson to be drawn from these simulations is clear: the approximated model provides a biased estimator. For the HDC sonar, we hence recommend to use the correct model to construct the MLE, even if it requires more computational time. Table VI and Fig. 8 present results when simulation duration is reduced to 1 minute. TABLE VI: PERFORMANCES OF app X ˆ AND X ˆ AND MINIMAL

Fig. 8 (

 8 Fig. 8 (b): A magnification of the clouds of estimates of Fig. 8(a).

  Our analysis is based on 500 Monte-Carlo runs. In the

	2) Scenario 2 : TABLE II: PERFORMANCES OF app X X  2 X ˆAND X ˆ AND MINIMAL STANDARD DEVIATIONS	scenarios used:  E P 0  T 3500   measurements are time stamped at 1 Hz rate: (m) and   T R P 0 3500   t k	(m). The   t k    1 ,
	X	Bias Approximated	Exact	ˆ Approximated	Exact		CRLB	with choose a reasonable standard deviation of the frequency s t 1   . The frequency 0 f is equal to 1000 Hz. We
	-2500 m	model 61.04	model 0.7	model 17.17	model 17.19	17 .28	measurement: simulated using the correct model. The total duration is 1300s. 25 . 0  Hz. All the measurements are 
	7000 m	40.37	0.06	5.45	5.42	5.49	In this configuration, we consider three targets, whose
	0 m/s	0.09	0.0015	0.02	0.02	0.02	respective trajectories are defined by
	-10 m/s	0.004	0.0003	0.006	0.006	0.006	X	1		 	000 , 10		000 , 25	25	 T 0	(the target sails eastwards),
								X	2		 	500 , 2	000 , 7	0		 T 10	(the target sails southwards),
								and	X	3		 	500 , 2	000 , 14	15		 T 13	(the target sails towards
								south-east). For each target, the respective performances are
								presented in tables, whose first column is the empirical bias of
								app X ˆ, while the second column is the empirical bias of X ˆ. The
								empirical standard deviations of the components of app X ˆ	and
								X ˆ are written in the third and fourth columns. The minimal
								standard deviations (computed from the CRLB with the
								correct model) of the components of X are the elements of the
								fifth column. Figures depict the estimate sets in the plane
								(x,y): the correct estimated positions are represented by blue
								dots, whereas the approximated estimated positions are given
								by red crosses. The 90% confidence ellipse is plotted in
								magenta.
										1) Scenario 1 :	X 	X	1
										TABLE I: PERFORMANCES OF app X ˆ AND X ˆ AND MINIMAL
																	STANDARD DEVIATIONS
										X					Bias Approximated	Exact	ˆ Approximated	Exact		CRLB
																	model	model	model	model
								-10000 m			11037	111.32	2671.4	2734.7 2510.5
								-25000 m			9112.7	103.33	2734.9	878.85	781
									25 m/s				3.77	0.14	1.67	0.23	0.15
										0 m/s					12.07	0.07	3.20	2.65	2.43

Fig.

3

: Target's trajectory, initial positions estimated with exact model (blue dots) and with approximated model (red crosses), and the 90% confidence ellipse (magenta). Note that the estimate clouds are croissant-shaped.

TABLE III

 III

							14000 m	1145.9	4.52	34.67	113.9	113.04
							15 m/s	3.9	0.003	0.19	0.18	0.18
							-13 m/s	4	0.009	0.15	0.26	0.26
							Fig. 5(a): Target's trajectory and initial positions estimated
							with exact model (blue dots) and with approximated model
								(red crosses).
		: PERFORMANCES OF app X ˆAND X ˆ	AND MINIMAL
		STANDARD DEVIATIONS
	X	Bias Approximated	Exact	ˆ Approximated	Exact		CRLB
		model	model	model	model
	-2500 m	4626.4	5.97	205.77	255.37 254.16

TABLE IV :

 IV PERFORMANCES OF

				X ˆ	AND	X ˆ	AND MINIMAL
				app
		STANDARD DEVIATIONS
		Bias			ˆ
	X	Approximated	Exact	Approximated	Exact		CRLB
		model	model	model		model
	-10000 m	1619.6	11.28	379 .12		452.41 449.98
	-25000 m	585.71	3.82	343.63		340.82 338.88
	25 m/s	0.55	0.001	0.21		0.21	0.20
	0 m/s	1.85	0.007	0.31		0.37	0.37

TABLE V :

 V PERFORMANCES OF

				X ˆ	AND	X ˆ	AND MINIMAL
				app
		STANDARD DEVIATIONS
		Bias			ˆ
	X	Approximated	Exact	Approximated	Exact		CRLB
		model	model	model		model
	-10000 m	376.69	3.31	59.98		58.67	59.48
	-25000 m	186.52	11.06	198.56		198.11 197.18
	25 m/s	0.06	0.006		0.11		0.11	0.11
	0 m/s	0.81	0.002		0.07		0.06	0.06
	Fig. 7 (a): Target's trajectory and initial positions estimated
	with exact model (blue dots) and with approximated model
	(red crosses) in a quadri-static configuration.

  Table VI and Fig. 8 present results when simulation duration is reduced to 1 minute.

	TABLE VI: PERFORMANCES OF app X ˆ AND	X ˆ AND MINIMAL
		STANDARD DEVIATIONS
	X	Bias Approximated	Exact	ˆ Approximated	Exact		CRLB
		model	model	model		model
	-10000 m	389.44	5.72	199.3		230.8	230.1
	-25000 m	130.05	14.72	1340.4		1337.9 1338.3
	25 m/s	0.02	0.004	1.02		1.02	1.02
	0 m/s	0.78	0.004	0.16		0.2		0.2
	Fig. 8(a): Target's trajectory and initial positions estimated
	with exact model (blue dots) and with approximated model
		(red crosses) in quadri-static configuration.

See https://en.wikipedia.org/wiki/History_of_radar