Antoine Lebon
email: antoine-lebon@etud.univ-tln.fr

Annie-Claude Perez
email: annie-claude.perez@univ-tln.fr

Claude Jauffret
email: claude.jauffret@univ-tln.fr

Dann Laneuville
email: dann.laneuville@naval-group.com

TMA from cosines of conical angles acquired by a towed array

Keywords: Target motion analysis, observability, Fisher information matrix, Cramér-Rao lower bound, conical angles, nonlinear estimation

This paper deals with the estimation of the trajectory of a target in constant velocity motion at an unknown constant depth, from measurements of conical angles supplied by a linear array. Sound emitted by the target does not navigate necessarily along a direct path toward the antenna, but can bounce on the sea bottom and/or on the surface. Observability is thoroughly analyzed to identify the ghost targets before proposing an efficient way to estimate the trajectory of the target of interest and of the ghost targets when they exist.

Introduction

Bearings-only target motion analysis (BOTMA) is a problem that has been widely studied and various solutions have been proposed in the literature: batch [START_REF] Nardone | Fundamental Properties and Performance of Conventional Bearings-Only Target Motion Analysis[END_REF][START_REF] Song | Practical guidance for homing missiles with bearings-only measurements[END_REF][START_REF] Clavard | Target Analysis of a Source in a Constant Turn from a Nonmaneuvering Observer[END_REF][START_REF] Farina | Target Tracking with Bearings-Only Measurements[END_REF][START_REF] Jauffret | Bearings-Only Maneuvering Target Motion Analysis from a Nonmaneuvering Platform[END_REF] or recursive filter (such as extended Kalman filter [START_REF] Zhang | Bearings-Only Filtering Using Uncorrelated Conversion Based Filters[END_REF][START_REF] Aidala | Utilization of Modified Polar Coordinates for Bearings-Only Tracking[END_REF][START_REF] Arulampalam | Performance of the Shifted Rayleigh Filter in Single-Sensor Bearings-Only Tracking[END_REF], unscented Kalman filter [START_REF] Laneuville | Recursive Bearings-Only TMA via Unscented Kalman Filter: Cartesian vs. Modified Polar Coordinates[END_REF], particle filter [START_REF] Arulampalam | Bearings-Only Tracking of Maneuvring Targets Using Particle Filters[END_REF], modified instrumental variable [START_REF] Zhang | Bearings-Only Target Motion Analysis via Instrumental Variable Estimation[END_REF][START_REF] Doğançay | On the efficiency of a bearings-only instrumental variable estimator for target motion analysis[END_REF][START_REF] Chan | Bearings-only tracking using data fusion and instrumental variables[END_REF]), or a mix of recursive and batch methods [START_REF] Kirubarajan | Bearings-only tracking of maneuvering targets using a batch-recursive estimator[END_REF]. Citing all the papers dealing with this topic is now a hard task. Among the abundant literature, most papers share the same assumption: the target is moving in a straight line with a constant speed, while the passive observer is maneuvering adequately in order to ensure observability of the target [START_REF] Jauffret | Observability in Passive Target Motion Analysis[END_REF][START_REF] Le Cadre | Discrete-time observability and estimability analysis for bearings-only target motion analysis[END_REF][START_REF] Nardone | Observability Criteria for Bearings-Only Target Motion Analysis[END_REF]. The bearings are the measurements.

In this paper we are concerned with the same problem except that the available measurements are the cosine of the relative bearings, called also conical angles because the target belongs to the cone of ambiguity whose revolution axis is the line along which the towed array is moving (see [START_REF] Urick | Principles of underwater sound[END_REF] p. 39). Implicitly, we consider a target moving in 3D at a constant and unknown depth in near field; in this case, the two more energetic rays are the direct and the reflected paths (bottom or surface). Most of the cases, the sound is bounced by the sea bottom. Therefore, we extend our analysis to surface and sea bottom bounced rays. Indeed, the array detects the cosine of the relative angle of the direction of arrival by a suitable spatial filtering such as beamforming, or more sophisticated techniques (see [START_REF] Van Trees | Optimum Array Processing. Part IV of Detection, Estimation, and Modulation Theory[END_REF]). In the near field, sound can propagate to the sensor array along the direct path and/or the bottom-reflected path, and/or the surface-reflected path. Most of the time, at most two rays coming from the same target are detected [START_REF] Urick | Principles of underwater sound[END_REF][START_REF] Abraham | Underwater Acoustic Signal Processing: Modeling, Detection, and Estimation; Modern Acoustic and Signal Processing[END_REF].

Unlike Gong [START_REF] Gong | Multipath target motion analysis: Properties and implication of the multipath process[END_REF] and Blanc-Benon [START_REF] Blanc-Benon | Target Motion Analysis from Bearing and Time-delay Measurements: The Use of Multipath[END_REF], who address the three-dimensional target motion analysis (TMA) from a sequence of time differences of arrival (TDOA) of a signal traveling by two different paths coupled with a sequence of azimuths, we assume in this paper that the available measurements are the cosines of the conical angles only. In [START_REF] Oh | Batch Processing through Particle Swarm Optimization for Target Motion Analysis with Bottom Bounce Underwater Acoustic Signals[END_REF], a similar problem was addressed, but observability was not studied. We will consider two situations: the first case is devoted to TMA when sound propagates along a non-direct path at each sampling time. This will be the topic of section 3: we will conduct observability analysis and identify all the ghost targets, given a set of noise-free measurements. We will prove that an assumption on the target's depth makes the target's trajectory observable, but not estimable (in the sense that the asymptotic performance given by the Cramér-Rao lower bound -CRLB -of the estimator of the depth is out of the physical constraints, that is, the source is navigating between the surface and the sea bottom).

In the fourth section, we will consider scenarios in which the antenna changes its own route. We will prove that the trajectory of the target is almost certainly observable.

In the fifth section, we will assume that sound will propagate along the direct path and the bottom-reflected path. The two rays will be assumed detected. Observability analysis will reveal that only three ghost targets at most exist without maneuvering of the antenna. We will check that in this case the depth is not "estimable". We will give a palliative allowing us to propose an estimator which is operationally acceptable, the price being a small bias. Convincing simulations will be given at the end of this section, proving that, even when the duration of the scenario is short, the estimated trajectory is very close to the true one. A conclusion ends the paper.

Notation and problem formulation

We consider two underwater vehicles moving at their own constant depth. The first mobile is a surface vessel or a submarine towing a horizontal sensor array, and the second one is the target of interest. Given a Cartesian coordinate system, the acoustic center of the array is located at time t at (𝑥 𝑂 (𝑡) 𝑦 𝑂 (𝑡) 𝑧 𝑂 ) 𝑇 . At the same time, the target is at (𝑥 𝑇 (𝑡) 𝑦 𝑇 (𝑡) 𝑧 𝑇 ) 𝑇 . The respective horizontal positions of the target of interest and of the center of the array at time t are denoted by 𝑃 𝑇 (𝑡) = (𝑥 𝑇 (𝑡) 𝑦 𝑇 (𝑡)) 𝑇 and 𝑃 𝑂 (𝑡) = (𝑥 𝑂 (𝑡) 𝑦 𝑂 (𝑡)) 𝑇 . The sea bottom depth (assumed to be a constant) is denoted D. The source is said to be endfire to the line array if its trajectory is in the same line as the array (which implies that the array and the source are at the same depth, and share the same route). It is broadside to the antenna if it navigates in the vertical plane orthogonal to the line array and passing by the acoustic center of the array. The sensors array detects the line of sight of the target; more precisely, an ad-hoc array processing (or spatial filtering) delivers at time t the cosine of the conical angle 𝑐 𝑎 (𝑡) given by cos(𝑐 𝑎 (𝑡)) = cos(𝜃(𝑡) -ℎ(𝑡)) cos(𝜙(𝑡)) ≜ 𝑚(𝑡) , where 𝜃(𝑡) and 𝜙(𝑡) are respectively the azimuth (or bearing) and the elevation of the path along which the sound emitted by the source propagates. The angle ℎ(𝑡) is the heading of the sensor array. Denoting the relative position coordinates of the source with reference to the acoustic center of the array by 𝑥 𝑂𝑇 (𝑡) = 𝑥 𝑇 (𝑡) -𝑥 𝑂 (𝑡) and 𝑦 𝑂𝑇 (𝑡) = 𝑦 𝑇 (𝑡) -𝑦 𝑂 (𝑡) , we have 𝜃(𝑡) = 𝑎𝑟𝑐𝑡𝑎𝑛(𝑥 𝑂𝑇 (𝑡), 𝑦 𝑂𝑇 (𝑡)). Figure 1 displays the different angles and the two actors (the observer reduced to the linear array, and the target). The ray of the sound (or signal) emitted by the source can be reflected by the bottom and/or the surface or travels in the surface or deep channel. The sound-speed profile makes the paths curve. In this paper, we will consider that the target is in the near field (the distance between the source and the array is less than 20 km), and the bottom depth is in the range [2000 m, 5000 m]. Due to the large curvature of the ray (about 80 km), we will approximate the path of the sound as a set of zigzags defined by the reflections on the bottom or on the surface. So, we implicitly use the Snell law widely employed in geometrical optics. An image-source is created whose depth 𝜁 𝑇 will be called "image-depth". A path is then defined by the triplet (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ), where We have to consider the depth difference between the array and the image-source defined by 𝜁 𝑂𝑇 ≜ 𝜁 𝑇 -𝑧 𝑂 if the ray has been reflected (by the sea bottom or by the surface), or 𝜁 𝑂𝑇 ≜ 𝑧 𝑇 -𝑧 𝑂 if the sound wave uses the direct path.

A general expression of 𝜁 𝑂𝑇 based on the triplet (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ) is given by 𝜁 𝑂𝑇 (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ) = -2𝛿𝑛 𝐵 (-1) 𝑛 𝑆 +𝑛 𝐵 𝐷 -𝑧 𝑂 + (-1) 𝑛 𝑆 +𝑛 𝐵 𝑧 𝑇 . Note that, given the path, the link between 𝜁 𝑂𝑇 (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ) and 𝑧 𝑇 is linear: 𝜁 𝑂𝑇 (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ) = 𝑎𝑧 𝑇 + 𝑏, the constants being a function of the triplet (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ), 𝐷, and 𝑧 𝑂 . Moreover, 𝜁 𝑂𝑇 (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ) is null if and only if the antenna and the target are navigating at the same depth (𝑧 𝑇 = 𝑧 𝑂 ), and sound is traveling in the direct path. In this case, cos(𝜙(𝑡)) = 1. For the sake of simplicity of the notations, we will simply subsequently denote 𝜁 𝑂𝑇 instead of 𝜁 𝑂𝑇 (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ).

For the above examples, we have 𝜁 𝑂𝑇 (1,0,0) = 𝑧 𝑇 -𝑧 𝑂 (direct path), 𝜁 𝑂𝑇 (+1,1,0) = 2𝐷 -(𝑧 𝑇 + 𝑧 𝑂 ) (bottom-reflected path), and 𝜁 𝑂𝑇 (+1,2,1) = 4𝐷 -(𝑧 𝑇 + 𝑧 𝑂 ) (bottom-surface-bottom reflected path). Note that 𝜁 𝑂𝑇 (𝛿, 𝑛 𝐵 , 𝑛 𝑆 ) can be negative (the image-source is above the surface). Consequently, the cosine of the elevation is

cos(𝜙(𝑡)) = √𝑥 𝑂𝑇 2 (𝑡)+𝑦 𝑂𝑇 2 (𝑡)
√𝑥 𝑂𝑇 2 (𝑡)+𝑦 𝑂𝑇 2 (𝑡)+𝜁 𝑂𝑇 2 (𝛿,𝑛 𝐵 ,𝑛 𝑆 )

.

Figure 3 (a) displays the cone of ambiguity, defined by the set of sources sharing the same cos(𝜙(𝑡)). In Figure 3 We assume that the source is moving in constant velocity (CV) motion during the scenario. Our challenge is to estimate its trajectory, i.e. the state vector defining it, 𝑋 ≜ (𝑥 𝑇 (𝑡 * ) 𝑦 𝑇 (𝑡 * ) 𝑧 𝑇 𝑥̇𝑇 𝑦̇𝑇) 𝑇 , for a chosen 𝑡 * , from noisy measurements.

We consider two situations:

1. Only one ray is detected by the array during the scenario; in this case, we have at each time t a measurement 𝑚(𝑡), given the path along which the wave propagates.

2. Two rays (traveling on two different paths) arrive at the sensors antenna. In this case, the available measurement at time t is a couple of measurements, say (𝑚 1 (𝑡), 𝑚 2 (𝑡)), given the two paths along which the wave propagates.

After the spatial filtering, the antenna supplies a noisy measurement of 𝑚(𝑡) or a noisy measurement of (𝑚 1 (𝑡), 𝑚 2 (𝑡)). The noisy measurements are regularly acquired at 𝑡 𝑘 = (𝑘 -1)Δ𝑡, 𝑘 ∈ {1, … , 𝑁}, for a fixed sampling time Δ𝑡.

Before attempting to estimate 𝑋, we must answer several questions: 1. Is the vector 𝑋 observable from the set of measurements {𝑚(𝑡), 𝑡 ∈ [0, 𝑇]}? Note that, in TMA problems, observability is often analyzed in continuous time (see [START_REF] Jauffret | Observability in Passive Target Motion Analysis[END_REF] and [START_REF] Nardone | Observability Criteria for Bearings-Only Target Motion Analysis[END_REF], for example), even though the noisy measurements are given in discrete time.

2.

If not, what are the ghost targets (those which could be detected at the same set of measurements {𝑚(𝑡), 𝑡 ∈ [0, 𝑇]})?

3. How do we make 𝑋 observable or with which new information? 4. Is the vector 𝑋 observable from the set of couples {(𝑚 1 (𝑡), 𝑚 2 (𝑡)), 𝑡 ∈ [0, 𝑇]}?

For the cases where 𝑋 is observable, we have then to compute the asymptotical performance of an unbiased estimator (given by the CRLB [START_REF] Steven | Fundamentals if statistical signal processing -Estimation theory[END_REF]), and the performance of our estimators in terms of bias and the covariance matrix. It is worth noting that using the FIM to prove observability can lead to a wrong conclusion [START_REF] Jauffret | Observability and Fisher Information Matrix in Nonlinear Regression[END_REF]. This why we use an analytic approach.

TMA from one ray

In this section, we consider the case where the array collects the cosine of a conical angle, the path of the ray being known by the operator. We start by analyzing the observability of the trajectory of the source of interest.

Observability analysis

Theorem 1. Let a linear antenna measure the cosine of a conical angle in the direction of a source, both in CV motion. The path of the sound emitted by the source is known, as is also the sea bottom depth.

1. If the target is broadside to the antenna, then the set of ghost targets is composed of virtual sources broadside to the antenna.

2. If the target is endfire to the antenna, the set of ghost targets is composed of virtual sources endfire to the antenna.

3. If the target has the same heading as the array (but is not endfire to it), then the set of ghost targets is composed of virtual targets having the same heading as the antenna. More precisely, the ghost image of each ghost target is moving on a cylinder whose axis is the antenna axis, and whose radius is a positive scalar 𝛽. The relative ghost target velocity is equal to 𝛽 times the target's velocity. The initial distance between the ghost image and the center of the antenna is equal to 𝛽 times the initial distance between the target-image and the center of the antenna. [START_REF] Farina | Target Tracking with Bearings-Only Measurements[END_REF] Preamble: In the following proof, we choose 𝑡 * = 0. Instead of working with the state vector 𝑋 = (𝑥 𝑇 (0) 𝑦 𝑇 (0) 𝑧 𝑇 𝑥̇𝑇 𝑦̇𝑇) 𝑇 , we will use the relative state vector of the image source, which is 𝑌 ≜ (𝑥 0𝑇 (0) 𝑦 𝑂𝑇 (0) 𝜁 𝑂𝑇 𝑥̇𝑂 𝑇 𝑦̇𝑂 𝑇 ) 𝑇 . The reason is that we are able to recover 𝑋 from 𝑌 without ambiguity.

We will prove this theorem in the special case where the heading of the antenna is equal to 0°, and the value 𝑦 𝑂𝑇 (𝑡) is positive. This can be easily obtained with an ad-hoc rotation of the whole scenario. This will simplify the expression of the measurement, without loss of generality.

Proof of Theorem 1. We are seeking the ghost target whose horizontal position at time t is (𝑥 𝐺 (𝑡) 𝑦 𝐺 (𝑡)) 𝑇 , detected in the same cosine of the conical angle, that is

𝑦 𝑂𝑇 (𝑡) √𝑥 𝑂𝑇 2 (𝑡)+𝑦 𝑂𝑇 2 (𝑡)+𝜁 𝑂𝑇 2 = 𝑦 𝑂𝐺 (𝑡) √𝑥 𝑂𝐺 2 (𝑡)+𝑦 𝑂𝐺 2 (𝑡)+𝜁 𝑂𝐺 2 , with 𝑥 𝑂𝐺 (𝑡) = 𝑥 𝐺 (𝑡) -𝑥 𝑂 (𝑡) ,
𝑦 𝑂𝐺 (𝑡) = 𝑦 𝐺 (𝑡) -𝑦 𝑂 (𝑡), and 𝜁 𝑂𝐺 is the image-depth of the ghost target. This equality is equivalent to

𝑦 𝑂𝑇 2 (𝑡) 𝑥 𝑂𝑇 2 (𝑡)+𝑦 𝑂𝑇 2 (𝑡)+𝜁 𝑂𝑇 2 = 𝑦 𝑂𝐺 2 (𝑡) 𝑥 𝑂𝐺 2 (𝑡)+𝑦 𝑂𝐺 2 (𝑡)+𝜁 𝑂𝐺 2 . ( 1
)
Note that because the target is moving (as is the ghost target also), the denominators of the left term and of the right term of (1) are two polynomial functions of degree 2.

Case 1: 𝑦 𝑂𝑇 (𝑡) is a zero function, i.e. ∀𝑡 𝑦 𝑂𝑇 (𝑡) = 0.

This means the source is broadside to the antenna: 𝑌 𝑇 = (𝑥 0𝑇 (0) 0 𝜁 𝑂𝑇 𝑥̇𝑂 𝑇 0) 𝑇 .

In this case, 𝑦 𝑂𝑇 (𝑡) = 0, ∀𝑡 ∈ [0, 𝑇]. Hence, the set of ghost targets is composed of the virtual targets broadside to the antenna: 𝑌 𝐺 = (𝑥 0𝐺 (0) 0 𝜁 𝑂𝐺 𝑥̇𝑂 𝐺 0) 𝑇 .

Case 2: 𝑦 𝑂𝑇 (𝑡) is not a zero function.

If 𝑦̇𝑂 𝑇 = 0, then 𝑦 𝑂𝑇 (𝑡) is a constant. To respect the degrees of the terms of (1), 𝑦 𝑂𝐺 (𝑡) is a constant too.

If 𝑦̇𝑂 𝑇 ≠ 0, then there is a root, say 𝑡 , such as 𝑦 𝑂𝑇 (𝑡) = 0, since 𝑦 𝑂𝑇 (𝑡) is a polynomial function of degree 1. Consequently, 𝑦 𝑂𝐺 (𝑡) = 0, and ∀𝑡 ≠ 𝑡 , 𝑦 𝑂𝐺 (𝑡) ≠ 0.

We deduce that, in both cases (𝑦̇𝑂 𝑇 = 0, and 𝑦̇𝑂 𝑇 ≠ 0), there exists a positive value 𝛽 such that 𝑦 𝑂𝐺 (𝑡) = 𝛽𝑦 𝑂𝑇 (𝑡).

(

) ⇔ {[𝑥 𝑂𝐺 2 (𝑡) + 𝑦 𝑂𝐺 2 (𝑡) + 𝜁 𝑂𝐺 2 ] -𝛽 2 [𝑥 𝑂𝑇 2 (𝑡) + 𝑦 𝑂𝑇 2 (𝑡) + 𝜁 𝑂𝑇 2 ]}𝑦 𝑂𝑇 2 (𝑡) = 0 ⇔ 𝑥 𝑂𝐺 2 (𝑡) + 𝜁 𝑂𝐺 2 = 𝛽 2 [𝑥 𝑂𝑇 2 (𝑡) + 𝜁 𝑂𝑇 2 ] 1
The quantity 𝑥 𝑂𝑇 2 (𝑡) + 𝜁 𝑂𝑇 2 can be equal to zero at any time, or at one time or never.

Subcase 1: ∀𝑡, 𝑥 𝑂𝑇 2 (𝑡) + 𝜁 𝑂𝑇 2 = 0.

Then 𝑥 𝑂𝑇 (𝑡) = 0, ∀𝑡 and 𝜁 𝑂𝑇 = 0. Note that this case is the one when the target is traveling in the endfire to the array and at the same depth as the antenna and the path is the direct one. For the same reason, 𝑥 𝑂𝐺 (𝑡) = 0, ∀𝑡 and 𝜁 𝑂𝐺 = 0. The set of ghost targets is hence composed of virtual targets traveling in the endfire to the array and at the same depth as the antenna.

Subcase 2: ∃𝑡 ̆ such that 𝑥 𝑂𝑇 2 (𝑡) + 𝜁 𝑂𝑇 2 ≠ 0.

We deduce from (2) that

{ 𝑥 𝑂𝐺 2 (0) = 𝛽 2 𝑥 𝑂𝑇 2 (0) + 𝛽 2 𝜁 𝑂𝑇 2 -𝜁 𝑂𝐺 2 𝑥 𝑂𝐺 (0)𝑥̇𝑂 𝐺 = 𝛽 2 𝑥 𝑂𝑇 (0)𝑥̇𝑂 𝑇 𝑥̇𝑂 𝐺 2 = 𝛽 2 𝑥̇𝑂 𝑇 2 (3-a) (3-b) (3-c) If 𝑥̇𝑂 𝑇 = 0, then 𝑌 𝐺 = (±√𝛽 2 𝑥 𝑂𝑇 2 (0) + 𝛽 2 𝜁 𝑂𝑇 2 -𝜁 𝑂𝐺 2
𝛽𝑦 𝑂𝑇 (0) 𝜁 𝑂𝐺 0 𝛽𝑦̇𝑂 𝑇 ) 𝑇 , for any positive constant 𝛽 and any positive constant 𝜁 𝑂𝐺 less than √𝛽 2 𝑥 𝑂𝑇 2 (0) + 𝛽 2 𝜁 𝑂𝑇 2 . Note that, when 𝑦̇𝑂 𝑇 = 0, the target is motionless relative to the center of the array (both have the same velocity); and when 𝑦̇𝑂 𝑇 ≠ 0, the target has the same heading as the array.

If 𝑥̇𝑂 𝑇 ≠ 0, then squaring the elements of (3-b), and using (3-c), we draw from (3-a) that 𝛽 2 𝜁 𝑂𝑇 2 = 𝜁 𝑂𝐺 2 . If 𝜁 𝑂𝑇 = 0, then 𝜁 𝑂𝐺 = 0, and the scalar 𝛽 can take any positive value; else

𝛽 = |𝜁 𝑂𝐺 | |𝜁 𝑂𝑇 | .
In both cases, the trajectory of a ghost target is defined by the state vector 𝑌 𝐺 = (±𝛽𝑥 𝑂𝑇 (0) 𝛽𝑦 𝑂𝑇 (0) 𝛽𝜁 𝑂𝑇 ±𝛽𝑥̇𝑂 𝑇 𝛽𝑦̇𝑂 𝑇 ) 𝑇 . □ Remarks:

1. When the source and the observer are at the same depth, and the path is direct, Theorem 1 recovers the conclusions given in [START_REF] Jauffret | Leg-by-leg Bearings-Only TMA without Observer Maneuver[END_REF].

2. The cases (1), ( 2) and (3) of Theorem 1 are "rare events" since the events of dealing with a source in endfire, broadside or having the same heading as the antenna during the scenario occur with a probability equal to 0. However, when the target has a trajectory close to one of these special cases, the estimates will have a poor behavior.

3. For cases (4), when the detected ray is not a direct path, for example, when the ray is bottom-reflected, a hypothesis about the source is sufficient to get one solution, corresponding to a ghost target. Indeed, if we suppose that the depth of the target is 𝑧 𝐴𝑠 (whereas the true value is 𝑧 𝑇 ), then we have 𝛽 = is too wide to be useful.

If the target and the antenna are at the same depth, 𝛽 can take any positive value.

Estimation of the trajectory

We run 500 Monte Carlo simulations for a typical scenario described as follows:

The observer starts from (0 0) 𝑇 at the depth 𝑧 𝑂 = 200m. Its speed and heading are respectively 5 m/s and 0°. The initial position of the target is (5000 7000) 𝑇 and its depth is 𝑧 𝑇 = 100m. Its route is 45° and its speed is 4 m/s.



The measurements are collected every 4 seconds (∆𝑡 = 4s). The scenario lasts 20 min.



The sea bottom depth is 4000 m. The detected ray is a bottom-reflected ray.



The assumed target depth is 𝑧 𝐴𝑠 = 200m (whereas the true one is 100 m).



First, the measurements have been corrupted with an additive Gaussian noise whose standard deviation is 𝜎 = 1.7 10 -2 .

Then we choose the least squares estimator, which is identical to the maximum likelihood estimator with these assumptions. Note that, in open literature about TMA, the confidence regions are given by the confidence ellipsoid obtained with the covariance matrix of the estimate. Since the maximum likelihood estimate is asymptotically efficient under nonrestrictive conditions, we use here the Cramér-Rao lower bound to compute such confidence regions.

The result of the simulation is presented in Table 1 and illustrated in Figure 4. Obviously, even if the assumption made on the target's depth makes the state vector observable, it remains inestimable: the hugeness of the diagonal elements of the CRLB does not allow this kind of TMA to be employed. We note in Figure 4 that the cloud of horizontal estimates is hyperbola-shaped. This is because the state vector is "weakly" estimable. The parametric equation of this hyperbola is

{ 𝑥(𝜔) = 𝜁 𝐴𝑠 sinh(𝜔) 𝑦(𝜔) = 𝜁 𝐴𝑠 𝑚 √1-𝑚 2 cosh(𝜔)
, with 𝜁 𝐴𝑠 = 2𝐷 -(𝑧 𝐴𝑠 + 𝑧 𝑂 ), and 𝑚 = 𝑦 𝑂𝑇 (0)

√𝑥 𝑂𝑇 2 (0)+𝑦 𝑂𝑇 2 (0)+𝜁 𝑂𝑇 2 .
We further reduced the standard deviation to 𝜎 = 1.7 10 -4 in order to appreciate the behavior of the MLE. With this (unrealistic) value, the MLE is efficient, as shown in Table 2 and in Figure 5 (which validates our observability analysis).

Our conclusion is that the state vector is not estimable, even though it is observable with an assumption on the target's depth. This is why we propose to maneuver the antenna in order to render the state vector observable with no assumption on the target's depth, and to augment the information about it.

TMA with one ray when the array maneuvers

In this section, the antenna maneuvers, i.e. it changes its own heading. We start by proving that the state vector is observable (without any assumption on the target's depth).

Then we have recourse to Monte Carlo simulations to evaluate the performance of the MLE.

Observability Analysis

Theorem 2. Suppose the antenna's trajectory is composed of two successive legs at constant velocity (however with the same speed). Let the target be in CV motion. The linear array acquires the conical angles of the wave emitted from the target, the path of the ray being known as well as the sea bottom depth. If the target is broadside or endfire to the antenna during a leg, then there is at most a ghost target. Otherwise, there is no ghost target.

Due to its length, the proof of this theorem is given in the Appendix.

Estimation

In this subsection, we present the result of 500 Monte Carlo simulations that are run to illustrate the behavior of the proposed estimators. First, we give the scenario used here.

The center of the array and the initial position of the source are respectively at (0 0 200) 𝑇 and (5000 7000 100) 𝑇 at the very beginning of the scenario. The speed of the array is a constant along the scenario and is equal to 5 m/s. The trajectory of the array is composed of two legs linked by an arc of a circle. The first leg lasts 1 min 40 s, during which the array's heading is 135°. Then the array turns to the right with a turn rate equal to 20°/min to adopt a new heading equal to 270°. The duration of the maneuver is hence equal to 6 min 44 s. The second leg lasts 5 min, so the total duration of the scenario is 13 min and 20 s. Meanwhile, the target is navigating with a heading equal to 45° and a speed of 4 m/s. The bottom depth is 𝐷 = 4000m. Measurements are acquired every ∆𝑡 = 4s, with 𝑡 𝑘 = (𝑘 -1)∆𝑡.

The noise vector 𝜀 𝑘 is assumed to be Gaussian, 0-mean and its standard deviation equal to 𝜎 = 1.7 10 -2 . The vectors 𝜀 𝑘 are also assumed to be temporally independent.

Again, we choose the least squares estimator.

Estimation of 𝑋

The 500 obtained estimates of the initial horizontal position are plotted in Figure 6, together with the trajectory of the target, the 90%-confidence ellipse and the trajectory of the array. Again, the view is from the sky.

The performance of the estimator (bias and standard deviation of each component) is presented in Table 3. A convenient way to evaluate the behavior of an estimator is to compute the so-called normalized estimation error squared (NEES) [START_REF] Bar-Shalom | Estimation and Tracking: Principles, Techniques and Software[END_REF], defined as 𝑁 𝑙 = (𝑋 ̂𝑙 -𝑋)

𝑇 𝐹(𝑋 ̂𝑙 -𝑋),

where 𝐹 is the FIM, and 𝑋 ̂𝑙 is the estimate computed at the 𝑙-th simulation. If 𝑋 ̂𝑙 is Gaussian-distributed with 𝑋 as the mathematical expectation and the CRLB as the covariance matrix, then 𝑁 𝑙 is Chi-square distributed with d degrees of freedom (𝜒 𝑑 2 ), where d is the dimension of 𝑋 (here 5). From the central limit theorem, the averaged NEES 𝑁 𝑆 ≜ .

From our simulations, we obtain 𝑁 𝑆 = 5.34.

In conclusion, the estimator can be declared efficient. However, the minimum standard deviation of the target's depth is not compatible with the physical constraints: with the standard deviation given in Table 3, the target could be up above the sea surface! Therefore, a palliative of this is to impose a depth on the target. Indeed, we saw in subsection 3.1 that a supposed depth creates a small bias in estimation of the horizontal position of the target. 4.2.2. Estimation of 𝑋 reduced when the depth of the target is fixed. Now, the third component of 𝑋 has not to be estimated. The new state vector is the denoted 𝑋 𝑟 = (𝑥 𝑇 (0) 𝑦 𝑇 (0) 𝑥̇𝑇 𝑦̇𝑇) 𝑇 . We impose that 𝑧 𝐴𝑠 = 200m (whereas the true depth is still 100 m). Hence, we introduce a bias. Figure 7 displays the position's estimates in the same manner as Figure 6. The bias is not visible to the naked eye. However, Table 4 reveals this bias, which may be acceptable in a real situation. Even though the averaged NEES (=7.31) is out of its 90% confidence interval, its value remains acceptable. The main interest of assuming the depth to be known is to economize on the CPU time, and reduce the standard deviation of the remaining components to estimate. We are in the presence of the well-known bias-variance tradeoff.

Estimation of the reduced state vector by the conventional BOTMA

In such a scenario, the conventional BOTMA can be run by neglecting the site effect, so by imposing that cos(𝜙(𝑡)) = 1, ∀𝑡. The (incorrect) noise-free measurement model is then

cos(𝛼(𝑡)) = cos(𝜃(𝑡) -ℎ(𝑡)).

The results are plotted in Figure 8. Obviously a huge bias appears, leading to an averaged NEES equal to 1960. More precisely, the bias on the components of the reduced state vector is (-3062.8 -2319.9 14.4 15.8) 𝑇 , rendering the BOTMA inoperative. Clearly, the conventional BOTMA cannot be recommended for the near field. This justifies a posteriori the interest in taking the site effect and the nature of the wave ray into account, as previously pointed out in the introduction of [START_REF] Oh | Batch Processing through Particle Swarm Optimization for Target Motion Analysis with Bottom Bounce Underwater Acoustic Signals[END_REF].

TMA from the direct path and the bottom-reflected path

We assume in this section that the sound wave emitted by the target travels on the direct path and the bottom-reflected path.

Observability

Theorem 3. Let a linear antenna and a source both be in CV motion.

The antenna acquires the cosines of the conical angles of the direct path and of the bottom-reflected path.

If the target is broadside to the array, then the set of ghost targets is uncountable: it is com-

posed of all the (virtual) targets at broadside to the array.

2. If the target is endfire to the antenna, the set of ghost targets is composed of virtual sources at endfire to the antenna. . If the depth of the antenna is equal to the depth of the source, then there is one single ghost target given by 𝑃 𝑂𝐺 (𝑡) = 𝑺𝑃 𝑂𝑇 (𝑡).

Proof of Theorem 3. With no loss of generality, we will again assume that the axis of the sensors array is toward North and that the target is in the half-space where the second component y of any vector is positive. A convenient rotation helps us to be in this case.

under the constraint that 𝑧 𝐺 is in [0, 𝐷].

Case 1: 𝑦 𝑂𝑇 (𝑡) is a zero function, i.e. ∀𝑡 𝑦 𝑂𝑇 (𝑡) = 0.

The target is broadside to the antenna, so any ghost targets will be too (see Case 1 in the proof of theorem 1).

Case 2: 𝑦 𝑂𝑇 (𝑡) is not a zero function.

From Case 2 of the proof of theorem 1, there is a positive scalar 𝛽 such that 𝑦 𝑂𝐺 (𝑡) = 𝛽𝑦 𝑂𝑇 (𝑡).

(

) 4
(5) ⟺ [𝑥 𝑂𝐺 2 (𝑡) + 𝑦 𝑂𝐺 2 (𝑡) + [2𝐷 -(𝑧 𝐺 + 𝑧 𝑂 )] 2 ] = 𝛽 2 [𝑥 𝑂𝑇 2 (𝑡) + 𝑦 𝑂𝑇 2 (𝑡) + [2𝐷 -(𝑧 𝑇 + 𝑧 𝑂 )] 2 ] (7)
Subtracting ( 7) from ( 6), we get

𝑧 𝑂𝐺 2 -[2𝐷 -(𝑧 𝐺 + 𝑧 𝑂 )] 2 = 𝛽 2 [𝑧 𝑂𝑇 2 -[2𝐷 -(𝑧 𝑇 + 𝑧 𝑂 )] 2 ].
Now, we simplify the expressions of these two terms:

𝑧 𝑂𝐺 2 -[2𝐷 -(𝑧 𝐺 + 𝑧 𝑂 )] 2 = -4(𝐷 -𝑧 𝐺 )(𝐷 -𝑧 𝑂 ) 𝑧 𝑂𝑇 2 -[2𝐷 -(𝑧 𝑇 + 𝑧 𝑂 )] 2 = -4(𝐷 -𝑧 𝑇 )(𝐷 -𝑧 𝑂 ).
We deduce from this that 𝛽 = √ 𝐷 -𝑧 𝐺 𝐷 -𝑧 𝑇 [START_REF] Arulampalam | Performance of the Shifted Rayleigh Filter in Single-Sensor Bearings-Only Tracking[END_REF] Note that 𝛽 = 1 iif 𝑧 𝐺 = 𝑧 𝑇 .

(

) ⟺ 𝑥 𝑂𝐺 2 (𝑡) + 𝑦 𝑂𝐺 2 (𝑡) + 𝑧 𝑂𝐺 2 = 𝛽 2 [𝑥 𝑂𝑇 2 (𝑡) + 𝑦 𝑂𝑇 2 (𝑡) + 𝑧 𝑂𝑇 2 ] ⟺ 𝑥 𝑂𝐺 2 (𝑡) -𝛽 2 𝑥 𝑂𝑇 2 (𝑡) = 𝛽 2 𝑧 𝑂𝑇 2 -𝑧 𝑂𝐺 2 6
Since 𝑥 𝑂𝐺 2 (𝑡) -𝛽 2 𝑥 𝑂𝑇 2 (𝑡) is a polynomial function of degree 2, ( 9) is equivalent to Reporting this in (10-a), we obtain finally 𝛽 2 𝑧 𝑂𝑇 2 = 𝑧 𝑂𝐺 2 , i.e.

{ 𝑥 𝑂𝐺 2 (0) = 𝛽 2 𝑥 𝑂𝑇 2 (0) + 𝛽 2 𝑧 𝑂𝑇 2 -𝑧 𝑂𝐺 2 𝑥 𝑂𝐺 (0)𝑥̇𝑂 𝐺 = 𝛽 2 𝑥 𝑂𝑇 (0)𝑥̇𝑂 𝑇 𝑥̇𝑂 𝐺 2 = 𝛽
𝛽 2 = ( 𝑧 𝑂𝐺 𝑧 𝑂𝑇 ) 2 (11)
If 𝑧 𝑇 = 𝑧 𝑂 , then 𝑧 𝐺 = 𝑧 𝑂 . In this case, 𝛽 = 1 , and consequently 𝑦 𝑂𝐺 (𝑡) = 𝑦 𝑂𝑇 (𝑡) and 𝑥 𝑂𝐺 2 (𝑡) = 𝑥 𝑂𝑇 2 (𝑡) from ( 9). The source's trajectory is observable up to the axial symmetry around the (Oy)-axis.

Eq. ( 8) and [START_REF] Zhang | Bearings-Only Target Motion Analysis via Instrumental Variable Estimation[END_REF] give us

𝐷-𝑧 𝐺 𝐷-𝑧 𝑇 = ( 𝑧 𝑂𝐺 𝑧 𝑂𝑇 ) 2 .
The unknown 𝑧 𝐺 is hence a root of the following equation of degree 2:

(𝑧 𝐺 -𝑧 𝑂 ) 2 - (𝑧 𝑇 -𝑧 𝑂 ) 2
𝐷-𝑧 𝑇 (𝐷 -𝑧 𝐺 ) = 0 which can be expanded as follows:

𝑧 𝐺 2 + 𝑧 𝐺 [-2𝑧 𝑂 + (𝑧 𝑇 -𝑧 𝑂 ) 2 𝐷-𝑧 𝑇 ] - 𝐷(𝑧 𝑇 -𝑧 𝑂 ) 2 𝐷-𝑧 𝑇 + 𝑧 𝑂 2 = 0.
Of course, 𝑧 𝑇 is a root of this equation. For this value, 𝑧 𝐺 = 𝑧 𝑇 , hence 𝛽 = 1.

The second root (𝑧 𝐺 itself) is hence 2𝑧 𝑂 -𝑧 𝑇 - (which is negative).

(
We deduce from this that:

1. when the target's depth is larger than the array's depth, there is a ghost whose depth is smaller than the array's depth, and vice-versa.

2. 𝛽, which is a positive coefficient, is equal to

𝐷-𝑧 𝑂 𝐷-𝑧 𝑇
, or 1.

Therefore, we have identified three ghost targets:

the first one is defined by 𝑋 𝑂𝐺 = (-𝑥 𝑂𝑇 (0) 𝑦 𝑂𝑇 (0) 𝑧 𝑂𝑇 -𝑥̇𝑂 𝑇 𝑦̇𝑂 𝑇 ) 𝑇 , the second is defined by 𝑋 𝑂𝐺 = (𝛽𝑥 𝑂𝑇 (0) 𝛽𝑦 𝑂𝑇 (0) -𝛽𝑧 𝑂𝑇 𝛽𝑥̇𝑂 𝑇 𝛽𝑦 ̇𝑂𝑇 ) 𝑇 , and the third by 𝑋 𝑂𝐺 = (-𝛽𝑥 𝑂𝑇 (0) 𝛽𝑦 𝑂𝑇 (0) -𝛽𝑧 𝑂𝑇 -𝛽𝑥̇𝑂 𝑇 𝛽𝑦 ̇𝑂𝑇 ) 𝑇 . □ Remark:

Most of the time, the depth of a submarine vehicle is under the operational constraint: values of 𝑧 𝑇 are in [0, 𝑧 𝑀𝑎𝑥 ] and 𝑧 𝑀𝑎𝑥 ≪ 𝐷 . For example, 𝑧 𝑀𝑎𝑥 = 400 m, while 𝐷 = 4000m.

The proof of the previous theorem must be adapted to this new constraint. If 𝑧 𝑂 > 𝐷 -𝐷√1 -𝑧 𝑀𝑎𝑥 𝐷 , then 𝑧 𝐺 > 𝑧 𝑀𝑎𝑥 . In this case, there is a unique ghost target given by 𝑋 𝐺 = (-𝑥 𝑇 (0) 𝑦 𝑇 (0) 𝑧 𝑇 -𝑥̇𝑇 𝑦̇𝑇) 𝑇 .

If 𝑧 𝑂 ≤ 𝐷 -𝐷√1 -𝑧 𝑀𝑎𝑥 𝐷 , then 𝑧 𝐺 ≤ 𝑧 𝑀𝑎𝑥 . In this case, there are three ghost targets: one is defined by 𝑋 𝐺 = (-𝑥 𝑇 (0) 𝑦 𝑇 (0) 𝑧 𝑇 -𝑥̇𝑇 𝑦̇𝑇) 𝑇 , the second is defined by 𝑋 𝐺 = (𝛽𝑥 𝑇 (0) 𝛽𝑦 𝑇 (0) 𝑓(𝑧 𝑇 ) 𝛽𝑥̇𝑇 𝛽𝑦̇𝑇) 𝑇 , and the third by 𝑋 𝐺 = (-𝛽𝑥 𝑇 (0) 𝛽𝑦 𝑇 (0) 𝑓(𝑧 𝑇 ) -𝛽𝑥̇𝑇 𝛽𝑦̇𝑇) 𝑇 .

Note that the operational constraint allows us to benefit from the following range:

𝐷-𝑧 𝑀𝑎𝑥 𝐷 ≤ 𝛽 ≤ 𝐷 𝐷-𝑧 𝑀𝑎𝑥
. For example, when 𝑧 𝑀𝑎𝑥 = 𝐷 10 , 0.9 ≤ 𝛽 ≤ 1.11. Consequently, the ghost target is very close to the target of interest.

Estimation of the trajectory

This section is devoted to the estimation of the target's trajectory, or in other words, the estimation of 𝑋 with 𝑡 * = 0 (the first time). Before going into detail, we compute the so-called Cramér-Rao lower bound to evaluate the asymptotical performance of any unbiased estimator.

We have considered two typical scenarios. In both, the array is assumed motionless (or, more realistically, all the mobiles are referenced to it) at the depth 𝑧 𝑂 = 200m, and the state vector defining the target's trajectory is given by the state vector 𝑋 = (5000 7000 100 2.83 2.83) 𝑇 . The standard deviation of the measurement is 𝜎 = 1.7 10 -2 . The total duration of the scenario is 5 min, and the sampling time is Δ𝑡 = 4s; consequently, the number of measurement couples is 𝑁 = 75.

In the first scenario, the bottom depth is 𝐷 = 2000 m, while in the second 𝐷 = 4000m.

Note that in the first scenario 𝛽 = 0.89 and in the second one 𝛽 = 0.97. The ghost target is hence very close to the target of interest.

Estimability

As pointed in section 1, the state vector 𝑋 is "estimable" if its asymptotical performance given by the CRLB is compatible with the physical constraints. Typically, if the minimum standard deviation defined by the square root of the third diagonal element of the CRLB (hence of the depth) is much larger than the depth, then 𝑋 is declared nonestimable.

First scenario

For this scenario, the square root of the diagonal of the CRLB 𝜎 𝐶𝑅𝐿𝐵 = ( 1. [START_REF] Le Cadre | Discrete-time observability and estimability analysis for bearings-only target motion analysis[END_REF] In both cases, the minimum standard deviations are huge. We can conclude that the state vector is not estimable. Computations of minimum standard deviations were made for various scenarios; in all, the state vector is not estimable.

A palliative of this is to fix the depth of the source at an arbitrary and realistic value, say 𝑧 𝐴𝑠 and compute the CRLB of the reduced state vector 𝑋 𝑟 ≜ (𝑥 𝑇 (0) 𝑦 𝑇 (0) 𝑥̇𝑇 𝑦̇𝑇) 𝑇 when we assume that 𝑧 𝑇 = 𝑧 𝐴𝑠 . For example, for 𝑧 𝐴𝑠 = 300m, the minimum standard deviations are 𝜎 𝐶𝑅𝐿𝐵 = (281.17 319.37 1.78 2.02) 𝑇 for the first scenario, and 𝜎 𝐶𝑅𝐿𝐵 = (130.1 115.3 0.80 0.71) 𝑇 for the second one. Therefore, we propose to estimate the state vector with this hypothesis (𝑧 𝐴𝑠 = 300m). In so doing, we introduce a bias. The next subsection gives us the result of the 500 Monte Carlo simulations.

Monte Carlo simulations

The computation of the maximum likelihood estimator (MLE) is made with the Gauss-Newton routine. No numerical issue was encountered.

First scenario

The performance of the MLE is summarized in Table 5. We have numerically computed the bias and the empirical standard deviation (given respectively in the second and third column of the table). We can see that the empirical standard deviation is very close to the one given by the CRLB. However, as expected, the MLE is biased (of course, there is no bias if we choose 𝑧 𝐴𝑠 = 𝑧 𝑇 ). In Figure 9, the 90%-confidence ellipse has been drawn, together with the cloud of the 500 estimates (in pink). Again, the performance is presented in Table 6. The bias of the estimator is similar to the one obtained for the first scenario. Only the empirical standard deviations of (𝑥 𝑇 (0) 𝑦 𝑇 (0)) 𝑇 are larger than the one computed from the CRLB. However, figure 10 shows us that the cloud of estimates is close to the true value and not spread. What is remarkable is the short duration and still the very good performance (in terms of accuracy) of the result. Numerous simulations (not reported here) were performed; all confirm the correct performance of the MLE. The shortness of the scenario is crucial, because everything that we propose here works properly under the condition that the sea bottom is a plane. During a short scenario, this assumption is likely.

Conclusion

In this paper, conical-angle TMA has been addressed, and various multipaths of sound have been taken into account. The sensor is a line array. Observability was analyzed deeply, allowing all the existing ghost targets to be identified. The main results are that, if the array detects one ray (corresponding to one path), the trajectory is not observable: the set of ghost targets is composed of trajectories that are homothetic to the trajectory of the target of interest, and their symmetrical images by the axial symmetry around the line array. If the array detects two rays (corresponding to two different paths), the number of ghost targets is reduced to three (except when the target is endfire or broadside to the antenna). When the antenna maneuvers, the target's trajectory is observable (apart from the special scenario where there is one single ghost target). Even for "observable" scenarios, the depth of the target is not estimable (its asymptotical standard deviation is huge). In these cases, we give a non-restrictive palliative that allows us to provide estimates close to the truth.

In the future, in this context, many problems remain to be faced: identification of the paths, maneuvering targets, and fusion of data collected by other sensors as in [START_REF] Payan | Passive Target Motion Analysis by Fusion of Linear Arrays and Sonobuoys in a Cluttered Environment[END_REF]. The problem of seeking a "good" maneuver of the observer, as it was solved in 2D environment [START_REF] Le Cadre | Optimization of the observer motion for bearings-only target motion analysis[END_REF][START_REF] Fawcett | Effect of course maneuvers on bearings-only range estimation[END_REF][START_REF] Passerieux | Optimal observer maneuver for bearings-only tracking[END_REF] will be addressed in the future. Some of these are already under investigation.

Appendix A

Proof of Theorem 2. The proof is made when the first leg is towards North (as previously).

Hence From Theorem 1, we have to consider the four following cases for each leg:

 the target is broadside to the antenna,  the target is endfire to the antenna,  the target has the same heading as the array (but is not endfire to it),  the other cases.

Note that if the target is in case (1) during the first leg, then in case (2) during the second one (provided that this situation is possible), the conclusion about observability will be the same as if the target is in case (2) during the first leg, then in case (1) during the second leg. To be convinced of this, we have just to reverse the time in the equation. This remark allows us to shorten the proof.

Case 1: the target is broadside to the antenna during the first leg.

Hence We skip the case where the target is in case (3) during the second leg. This will be treated later. Therefore, we have now to consider the other cases during the second leg. There are two possibilities for the ghost targets: those whose trajectories are defined by (i) 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑃 𝑂𝑇 (𝑡), and those whose trajectories are given by (ii) 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑺 2 𝑃 𝑂𝑇 (𝑡), both for 𝑡 ≥ 𝜏.

The derivative of (i) is

𝑉 𝐺 -𝑉 2 = 𝛽𝑉 𝑇 -𝛽𝑉 2 , hence 𝑉 𝐺 = 𝛽𝑉 𝑇 + (1 -𝛽)𝑉 2 . ⟺ [ 𝑐 𝐺 𝑣 ] = 𝛽 [ 𝑐 𝑇 𝑣 ] + (1 -𝛽)𝑣 [ sin(𝛼) cos(𝛼) ],
which implies that (1 -𝛽) cos(𝛼) = 1 -𝛽. Since cos(𝛼) ≠ 1, 𝛽 = 1. There is no ghost given by (i).

The derivative of (ii) is 𝑉 𝐺 -𝑉 2 = 𝛽𝑺 We deduce that 𝛽 = 𝑣 Again, we skip the case where the target is in case (3) during the second leg. This will be treated later. So, we have now to consider the other cases during the second leg. There are two possibilities for the ghost targets: those whose trajectories are defined by (i) 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑃 𝑂𝑇 (𝑡) and those whose trajectories are given by (ii) 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑺 2 𝑃 𝑂𝑇 (𝑡), both for 𝑡 ≥ 𝜏.

The derivative of (i) is 𝑉 𝐺 -𝑉 2 = 𝛽𝑉 𝑇 -𝛽𝑉 This implies that 𝛼 = 0, which must be rejected by assumption. There is no ghost.

The other cases:

In the other cases, the motion of ghost targets is defined during the first leg by

Eq. ( 16) implies that 𝑃 𝑂𝑇 (𝜏) is an eigenvector of 𝑺 2 , with the eigenvalue Hence, the set of ghost targets is reduced to those whose positions at time 𝜏 are given by (14-a) or (14-b), and (15-a). Now suppose that a ghost target satisfies (14-b) and (15-a). By

Figure 1 .

 1 Figure 1. A typical scenario, viewed from the sky.



  𝛿 indicates the direction of the path of the sound emitted by the source: if the path is toward the surface, 𝛿 = -1, otherwise 𝛿 = +1,  𝑛 𝐵 is the number of bottom reflections, and  𝑛 𝑆 is the number of surface reflections.

Figure 2 Figure 2 .

 22 Figure 2 illustrates three different paths

Figure 3 .

 3 Figure 3 (a) displays the cone of ambiguity, defined by the set of sources sharing the same cos(𝜙(𝑡)). In Figure 3 (b), we have plotted a direct ray and a bottom-bounced ray, which allows us to figure out the various angles with which we will work.

Figure 4 .

 4 Figure 4. The cloud of estimated position (in green), a piece of the hyperbola (intersection of the cone of ambiguity and the plane 𝑧 = 𝑧 𝐴𝑠 (= 200 𝑚), for 𝜎 = 1.7 10 -2 .

Figure 5 .

 5 Figure 5. The cloud of estimated position (in green) for 𝜎 = 1.7 10 -4 . The cloud is no longer hyperbola-shaped. The small black segment is the 90%-confidence ellipsoid.

  5000 7000 100 2.83 2.83) 𝑇 . The array is assumed to measure the cosines of the conical angles of the bottom-reflected path given by 𝑚(𝑡 𝑘 ) = 𝑦 𝑂𝑇 (𝑡 𝑘 ) √𝑥 𝑂𝑇 2 (𝑡 𝑘 ) + 𝑦 𝑂𝑇 2 (𝑡 𝑘 ) + [2𝐷 -(𝑧 𝑇 + 𝑧 𝑂 )] 2 + 𝜀 𝑘 .

Figure 6 .

 6 Figure 6. The cloud of the 500 initial positions estimates and the 90%-confidence ellipse.

Figure 7 .

 7 Figure 7. The cloud of the 500 initial positions estimates with the reduced state vector and the 90%-confidence ellipse.

Figure 8 .

 8 Figure 8. The cloud (in green) of the 500 initial positions estimates given by the classic BOTMA together with the 90%-confidence ellipse.

3 .

 3 If the route of the antenna and the route of the target are parallel, then the set of ghost targets is uncountable: at each depth 𝑧 𝐺 , there are two ghost targets moving on a cylinder whose axis is the antenna axis, and the radius is a positive scalar 𝛽 = √ 𝐷-𝑧 𝐺 𝐷-𝑧 𝑇 . The relative ghost target velocity is equal to 𝛽 times the target's velocity. The initial distance between the ghost image and the center of the antenna is equal to 𝛽 times the initial distance between the ghost image and the center of the antenna. 4. If the route of the antenna and the route of the target are not parallel, then there are three ghost targets whose motion relative to the antenna is 𝑃 𝑂𝐺 (𝑡) = 𝑺𝑃 𝑂𝑇 (𝑡), 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑃 𝑂𝑇 (𝑡), and 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑺𝑃 𝑂𝑇 (𝑡), where 𝑺 is the matrix of the axial symmetry around the line of the antenna, and 𝛽 ≜ 𝐷-𝑧 𝑂 𝐷-𝑧 𝑇

  of a ghost target.Second case ẋO T ≠ 0Using (10-c), and squaring the terms of (10-b), we get 𝑥 𝑂𝐺 2 (0) = 𝛽 2 𝑥 𝑂𝑇 2 (0).

2 𝐷

 2 𝑧 𝑇 -𝑧 𝑂 ) 2 𝐷-𝑧 𝑇 ≜ 𝑧 𝐺 . We can check readily that 𝑧 𝐺 -𝑧 𝑂 = 𝑧 𝑂 -𝑧 𝑇 -(𝑧 𝑇 -𝑧 𝑂 ) 2 𝐷-𝑧 𝑇 = (𝑧 𝑂 -𝑧 𝑇 )(𝐷-𝑧 𝑇 )-(𝑧 𝑇 -𝑧 𝑂 )

First, we use

  the fact that the function 𝑢 ↦ 𝑓(𝑢) ≜ 2𝑧 𝑂 -𝑢 -(𝑢-𝑧 𝑂 ) 2 𝐷-𝑢 is an involution, i.e. 𝑓(𝑓(𝑢)) = 𝑢.Since 𝑓(0) = 2𝑧 𝑂 -Now the question is: what are the values of 𝑧 𝑂 for which the following inequality holds: 𝑀𝑎𝑥 , the greatest value of 𝑧 𝑂 guaranteeing that 2𝑧 𝑂less than 𝑧 𝑀𝑎𝑥 ).

Figure 9 . 2 .

 92 Figure 9. The location of the sensor array (in black), the cloud of the 500 estimates and the 90%confidence ellipse when 𝐷 = 2000m, 𝑧 𝐴𝑠 = 300m, and 𝑧 𝑇 = 100m. The symmetrical cloud is plotted too. 2. Second scenario: Bottom depth 𝐷 = 4000m.

Figure 10 .

 10 Figure 10. The location of the sensor array (in black), the cloud of the 500 estimates of the initial positions, and the 90%-confidence ellipse when 𝐷 = 4000m, 𝑧 𝐴𝑠 = 300m, and 𝑧 𝐴𝑠 = 100m, together with the symmetrical cloud.

  assume 𝛼 ≠ 𝑘𝜋.

  when 𝑡 ≤ 𝜏, { 𝑃 𝑂𝐺 (𝑡) = 𝛽 1 𝑃 𝑂𝑇 (𝑡) or 𝑃 𝑂𝐺 (𝑡) = 𝛾 1 𝑺 1 𝑃 𝑂𝑇 (𝑡) (12-a) (12-b) and during the second leg by when 𝑡 ≥ 𝜏, { 𝑃 𝑂𝐺 (𝑡) = 𝛽 2 𝑃 𝑂𝑇 (𝑡) or 𝑃 𝑂𝐺 (𝑡) = 𝛾 2 𝑺 2 𝑃 𝑂𝑇 (𝑡) (13-a) (13-b) Hence at time 𝜏, the position of a ghost target is 𝑃 𝑂𝐺 (𝜏) = 𝛽 1 𝑃 𝑂𝑇 (𝜏) (14-a) or 𝑃 𝑂𝐺 (𝜏) = 𝛾 1 𝑺 1 𝑃 𝑂𝑇 (𝜏) (14-b) and 𝑃 𝑂𝐺 (𝜏) = 𝛽 2 𝑃 𝑂𝑇 (𝜏) (15-a) or 𝑃 𝑂𝐺 (𝜏) = 𝛾 2 𝑺 2 𝑃 𝑂𝑇 (𝜏) (15-b) Of course, (14-a) and (14-b) are not compatible, and neither are (15-a) and (15-b). Now, let us show that (14-a) is not compatible with (15-b): Indeed, if 𝑃 𝑂𝐺 (𝜏) = 𝛽 1 𝑃 𝑂𝑇 (𝜏) = 𝛾 2 𝑺 2 𝑃 𝑂𝑇 (𝜏), then 𝛽 1 𝛾 2 𝑃 𝑂𝑇 (𝜏) = 𝑺 2 𝑃 𝑂𝑇 (𝜏)

  eigenvalue is equal to 1, i.e. 𝛾 2 = 𝛽 1 . Hence 𝑃 𝑂𝑇 (𝜏) is in the second leg.

  𝑧 𝑇 +𝑧 𝑂 ) where 𝑧 𝑀𝑎𝑥 is the largest depth of a submarine vehicle. Typically, in deep water, 𝐷 ≥ 4000m. A reasonable choice of 𝑧 𝑀𝑎𝑥 could be 400 m. We can then have a range of 𝛽 : [𝛽 𝑀𝑖𝑛 , 𝛽 𝑀𝑎𝑥 ] = [ For instance, when the depths of the antenna and the target are respectively 200 m and 100 m, we have [𝛽 𝑀𝑖𝑛 , 𝛽 𝑀𝑎𝑥 ] = [0.974,1.013]. In this way, we bound the set of ghost targets, and we can expect that the bias induced by a wrong choice of 𝑧 𝐴𝑠 is very low. 𝑧 𝐴𝑠 -𝑧 𝑂 has the same sign as 𝑧 𝑇 -𝑧 𝑂 : if 𝑧 𝑇 > 𝑧 𝑂 , then 𝑧 𝑂 < 𝑧 𝐴𝑠 ≤ 𝑧 𝑀𝑎𝑥 , and [𝛽 𝑀𝑖𝑛 , 𝛽 𝑀𝑎𝑥 ] = ]0, In both cases, the range [𝛽 𝑀𝑖𝑛 , 𝛽 𝑀𝑎𝑥 ]

	𝛽 𝑀𝑎𝑥 = ]. 4. For cases (4) again, with a direct path, if the target is not at the same depth as the 2𝐷-𝑧 𝑂 2𝐷-(𝑧 𝑀𝑎𝑥 +𝑧 𝑂 ) 2𝐷-(𝑧 𝑇 +𝑧 𝑂 ) , and the minimum value is 𝛽 𝑀𝑖𝑛 = 2𝐷-(7600-𝑧 𝑂 8000-(𝑧 𝑇 +𝑧 𝑂 ) , 8000-𝑧 𝑂 8000-(𝑧 𝑇 +𝑧 𝑂 ) antenna, 𝛽 = 𝑧 𝐴𝑠 -𝑧 𝑂 𝑧 𝑇 -𝑧 𝑂 . Because 𝛽 is a positive number, 𝑧 𝑀𝑎𝑥 -𝑧 𝑂 ]; if 𝑧 𝑇 < 𝑧 𝑂 , 𝑧 𝑇 -𝑧 𝑂 then 0 ≤ 𝑧 𝐴𝑠 < 𝑧 𝑂 , and [𝛽 𝑀𝑖𝑛 , 𝛽 𝑀𝑎𝑥 ] = [0, 𝑧 𝑂 𝑧 𝑇 -𝑧 𝑂 ].

2𝐷-(𝑧 𝐴𝑠 +𝑧 𝑂 )

2𝐷-(𝑧 𝑇 +𝑧 𝑂 ) , whose biggest value

Table 1 .

 1 Performance of the estimator of the reduced state vector when 𝜎 = 1.7 10 -2 , in terms of bias, sample standard deviation and the one given by the square root of the di-

	agonal of the CRLB
	𝑿 𝒓	Bias	𝝈 𝒔𝒂𝒎𝒑	𝝈 𝑪𝑹𝑳𝑩
	5000 m	-3525	6962	13356
	7000 m	-2367	4052	5599
	2.83 m/s	-1.37	1.81	4.35
	2.83 m/s	0.53	1.62	2.75

Table 2 .

 2 Performance of the estimator of the reduced state vector with 𝜎 = 1.7 10 -4 .

	𝑿 𝒓	Bias	𝝈 𝒔𝒂𝒎𝒑	𝝈 𝑪𝑹𝑳𝑩
	5000 m	60.40	138.67	133.56
	7000 m	88.20	58.42	55.99
	2.83 m/s	0.043	0.044	0.044
	2.83 m/s	0.037	0.028	0.028

Table 3 .

 3 Performance of the estimator of the plain state vector

	𝑿	Bias	𝝈 𝒔𝒂𝒎𝒑	𝝈 𝑪𝑹𝑳𝑩
	5000 m	-44.77	854.72	868.12
	7000 m	-68.16	1162.1	1173.60
	100 m	7.14	558.55	545.99
	2.83 m/s	0.092	1.67	1.65
	2.83 m/s	0.194	2.72	2.68

Table 4 .

 4 Performance of the estimator of the reduced state vector.

	𝑿 𝒓	Bias	𝝈 𝒔𝒂𝒎𝒑	𝝈 𝑪𝑹𝑳𝑩
	5000 m	65.40	655.61	606.41
	7000 m	93.05	831.75	762.56
	2.83 m/s	0.034	1.71	1.58
	2.83 m/s	0.063	2.76	2.52

  10 6 1.59 10 6 8.22 10 5 637.9 646.1 ) 𝑇 .

	2. Second scenario
	With	the	bottom	depth,	things	are	not	much	better,	since
	𝜎 𝐶𝑅𝐿𝐵 = ( 6.59 10 5 8.96 10 5 9.73 10 5 352.6 362.2 ) 𝑇 .

Table 5 .

 5 Performance of the estimator of the reduced state vector.

	𝑿 𝒓	Bias	𝝈 𝒔𝒂𝒎𝒑	𝝈 𝑪𝑹𝑳𝑩
	5000 m	401.12	281.85	281.17
	7000 m	557.24	330.87	319.37
	2.83 m/s	0.13	1.58	1.78
	2.83 m/s	0.12	1.81	2.02

Table 6 .

 6 Performance of the estimator.

	𝑿 𝒓	Bias	𝝈 𝒔𝒂𝒎𝒑	𝝈 𝑪𝑹𝑳𝑩
	5000 m	306.81	219.28	130.08
	7000 m	432.46	276.61	115.26
	2.83 m/s	0.18	0.74	0.80
	2.83 m/s	0.18	0.66	0.71

  Can the target be endfire to the antenna? If so, the target has the same heading as the antenna during the second leg or, in other words, 𝑉 𝑇 -𝑉 2 = 𝜆𝑉 2 , and 𝑃 𝑂𝑇 (𝑡), which is equal to 𝑃 𝑂𝑇 (𝑡) = 𝑃 𝑂𝑇 (𝜏) + (𝑡 -𝜏)(𝑉 𝑇 -𝑉 2 ), is collinear with 𝑉 2 , whenever 𝑡 ≥ 𝜏.

	𝑃 𝑇 (0) = [ 𝑥 𝑇 (0) 0 𝑃 𝑂𝑇 (𝑡) = [ 𝑥 𝑇 (0) + 𝑡𝑐 𝑇 ] for 𝑡 ≤ 𝜏 . The ghost targets are also in the broadside, hence ] , and 𝑉 𝑂𝑇 = [ 𝑐 𝑇 ] during the first leg, which implies 0 0 𝑃 𝑂𝐺 (𝑡) = [ 𝑥 𝐺 (0) + 𝑡𝑐 𝐺 0 ] and 𝑉 𝑂𝐺 = [ 𝑐 𝐺 0 ]for 𝑡 ≤ 𝜏.
	condition cannot be satisfied since 𝑉 𝑇 = [	𝑐 𝑇 𝑣	], and 𝑉 2 = [ ±𝑣 0	]. There is no ghost.	The first

  2 𝑉 𝑇 -𝛽𝑉 2 , hence 𝑉 𝐺 = 𝛽𝑺 2 𝑉 𝑇 + (1 -𝛽)𝑉 2 .

	⟺ [	𝑐 𝐺 𝑣	] = 𝛽 [ -𝑐 𝑇 cos(2𝛼) + 𝑣 sin(2𝛼) 𝑐 𝑇 sin(2𝛼) + 𝑣 cos(2𝛼)	] + (1 -𝛽)𝑣 [ cos(𝛼) sin(𝛼)	].

  𝑡 ≤ 𝜏 . During this first leg, the ghost targets are also endfire to the antenna, so

			1-cos(𝛼)
			𝑐 𝑇 sin(2𝛼)+𝑣 cos(2𝛼)-𝑣 cos(𝛼)
	Hence 𝑃 𝑇 (0) = [	0 𝑦 𝑇 (0) ], and 𝑉 𝑂𝑇 = [ 𝑐 𝑇 0	], which implies that 𝑃 𝑂𝑇 (𝑡) = [	0 𝑦 𝑇 (0) + 𝑡𝑐 𝑇	] for
	𝑃 𝑂𝐺 (𝑡) = [ 𝑦 𝐺 (0) + 𝑡𝑐 𝐺 0	]for 𝑡 ≤ 𝜏, and 𝑉 𝐺 -𝑉 1 = [ 𝑐 𝐺 0	].

.

One ghost exists if 𝑐 𝑇 sin(2𝛼) + 𝑣 cos(2𝛼) -𝑣 cos(𝛼) is a positive quantity. If so, we then compute 𝑐 𝐺 . There is one ghost at most. Case 2: the target is endfire to the antenna during the first leg.

  2 , hence 𝑉 𝐺 = 𝛽𝑉 𝑇 + (1 -𝛽)𝑉 2 . 𝑉 2 = 𝛽𝑺 2 𝑉 𝑇 -𝛽𝑉 2 , hence 𝑉 𝐺 = 𝛽𝑺 2 𝑉 𝑇 + (1 -𝛽)𝑉 2 . 𝑐 𝑇 sin(2𝛼) + 𝑣 sin(2𝛼) -𝑣 sin(𝛼) is a negative quantity. If so, we then compute 𝑐 𝐺 . There is one ghost at most. Case 3: the target has the same heading as the array (but is not endfire to it) As in case (2), 𝑉 𝑂𝑇 = [ 0 𝑐 𝑇 ], but here, the first component of 𝑃 𝑂𝑇 (𝑡) is not zero. Hence, Can the target be broadside to the antenna? The answer is positive if 𝑉 2 ⊥ 𝑉 𝑇 and 𝑉 𝑂𝑇 is collinear to 𝑃 𝑂𝑇 (𝑡) , when 𝑡 ≥ 𝜏 . The first condition implies that 𝑉 2 = [ ±𝑣 0 ] . Since 𝑃 𝑂𝑇 (𝑡) = 𝑃 𝑂𝑇 (𝜏) + (𝑡 -𝜏)(𝑉 𝑇 -𝑉 2 ), the second condition is satisfied if 𝑃 𝑂𝑇 (𝜏) is collinear to 𝑉 𝑇 -𝑉 2 . This is not the case when the first component of 𝑃 𝑂𝑇 (𝜏) is zero, while the first component of 𝑉 𝑇 -𝑉 2 is ±𝑣. So, we have now to consider the other cases during the second leg. There are two possibilities for the ghost targets: those whose trajectories are defined by (i) 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑃 𝑂𝑇 (𝑡) and those whose trajectories are given by (ii) 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑺 2 𝑃 𝑂𝑇 (𝑡), both for 𝑡 ≥ 𝜏. The derivative of (i) is 𝑉 𝐺 -𝑉 2 = 𝛽𝑉 𝑇 -𝛽𝑉 2 , hence 𝑉 𝐺 = 𝛽𝑉 𝑇 + (1 -𝛽)𝑉 2 or, in other words, If 𝑃 𝑂𝐺 (𝑡) = 𝛽𝑺 2 𝑃 𝑂𝑇 (𝑡), then 𝑉 𝐺 = 𝛽𝑺 2 𝑉 𝑇 + (1 -𝛽)𝑉 2

	We	deduce	that			𝛽 = -𝑣	sin(𝛼) 𝑐 𝑇 sin(2𝛼)+𝑣 sin(2𝛼)-𝑣 sin(𝛼)	.	One	ghost	exists	if
	𝑉 𝑇 = [ this case, 𝑉 𝑂𝐺 = [ 0 ], and the target cannot be endfire to the antenna during the second leg. In 𝑐 𝑇 + 𝑣 0 𝛽𝑐 𝑇 ], hence 𝑉 𝐺 = [ 0 𝛽𝑐 𝑇 + 𝑣 ].
	[ 𝛽𝑐 𝑇 + 𝑣 0	] = [ 𝛽(𝑐 𝑇 + 𝑣) 0	sin(𝛼) cos(𝛼) ] + (1 -𝛽)𝑣 [
	[ 𝛽𝑐 𝑇 + 𝑣 0	] = 𝛽(𝑐 𝑇 + 𝑣) [	sin(2𝛼) cos(2𝛼)	] + (1 -𝛽)𝑣 [	sin(𝛼) cos(𝛼) ].
	⟺ [	0 𝑐 𝐺 + 𝑣	] = 𝛽 [	0 𝑐 𝑇 + 𝑣	] + (1 -𝛽)𝑣 [	sin(𝛼) cos(𝛼) ].
	We deduce that 𝛽 = 1. There is no ghost.
	Now, 𝑉 𝐺 -⟺ [ differentiating (ii) gives us 0 𝑐 𝐺 + 𝑣 ] = 𝛽(𝑐 𝑇 + 𝑣) [ sin(2𝛼) cos(2𝛼) ] + (1 -𝛽)𝑣 [ sin(𝛼) cos(𝛼) ].

⟹ 𝛽(𝑐 𝑇 + 𝑣) sin(2𝛼) + (1 -𝛽)𝑣 sin(𝛼) = 0. ] . We conclude that 𝛽 = 1 , i.e. there is no ghost.

the same computation, we conclude that 𝑃 𝑂𝑇 (𝜏) is in the first leg, which is impossible since 𝑃 𝑂𝑇 (𝜏) is in the second leg.

We have proven that (14-a) and (15-a) only are compatible. It follows that a ghost target verifies these two equalities (given by (12-a) and (13-a)):

Hence, 𝛽 1 = 𝛽 2 . Now taking the derivative of the two members of (12-a) and of (13-a), we obtain 𝑉 𝐺 = 𝛽 1 (𝑉 𝑇 -𝑉 1 ) + 𝑉 1 = 𝛽 1 (𝑉 𝑇 -𝑉 2 ) + 𝑉 2 , which is equivalent to

Putting this value in (12-a) or in (13-a), we finally get 𝑃 𝐺 (𝑡) = 𝑃 𝑇 (𝑡). The "ghost" is the target of interest. In conclusion, there is no ghost target. □


Jérémy Payan
email: jeremy-payan@etud.univ-tln.fr

Antoine Lebon
email: antoine-lebon@etud.univ-tln.fr

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

Claude Jauffret
email: jauffret@univ-tln.fr

Dann Laneuville
email: dann.laneuville@naval-group.com

Passive Target Motion Analysis by Fusion of Linear Arrays and Sonobuoys in a Cluttered Environment

Keywords: Fisher information matrix, multipath, observability, probabilistic data association with maximum likelihood estimation (ML-PDA), sonobuoy, target motion analysis, time difference of arrival

This paper is devoted to the analysis of the motion of a target using information from two kinds of cooperative maritime sensors: A wireless network of sonobuoys detecting a signal emitted by a source in motion with constant velocity; Vertical antennas measuring the cosine of the elevation angles of the received signal. We prove that the trajectory of the source is observable, under non-restrictive assumptions concerning the scenario. After thresholding, the data are surrounded by false alarms; therefore, a probabilistic data association model is employed. The joint exploitation of the measurements of the time differences of arrival together with the measurements of the cosines of the elevation angles allow estimating the trajectory of the source. The empirical performance of the maximum likelihood estimator (MLE), evaluated by extensive simulations, reaches the asymptotic performance given by the Cramér-Rao lower bound. Finally, we extend our study to the case where the energy of the data is available, after thresholding. Again, the MLE is efficient and its performance is significantly improved.

I. INTRODUCTION

In this paper, we are concerned with passive threedimensional submarine target tracking. Passive target tracking consists of estimating the parameters that define the trajectory of the target, from sensors which do not emit any signal. Each sensor receives the signal emitted by the moving target. Passive target tracking is of significant importance in coastal waters or submarine surveillance systems [START_REF] Sittler | An optimal data association problem in surveillance theory[END_REF].

Passive target tracking may be performed by using one sensor, for example in [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], but a network of sensors may also be used. Generally, the sensors are spatially separated from each other. Then, different kinds of sensors may be used in the network. We will focus specifically on two of them: vertical linear arrays and sonobuoys [START_REF] Urick | Principles of Underwater Sound[END_REF]. Sonobuoys measure the time difference of arrival (TDOA) between the signal originated by the target and that originated by a reference sonobuoy [START_REF] Arnold | Target parameter estimation using measurements acquired with a small number of sensors[END_REF]. The measurements acquired by a vertical linear array are the cosines of the elevation angles of the target originated wave.

Each of these two kinds of sensors has its own advantages and disadvantages regarding the asymptotic performance of the estimation. The asymptotic performance is given by the Cramér-Rao lower bound (CRLB) [START_REF] Van Trees | Detection, Estimation, and Modulation Theory, Part I[END_REF]. Due to their rather similar depths, the sonobuoys are very unlikely to estimate the depth of the target with relatively low variance [START_REF] Lo | Broadband passive acoustic technique for target motion parameter estimation[END_REF], unless there is a large number of sonobuoys. But vertical linear arrays can estimate the target depth with low variance. In view of these considerations, it would be interesting to combine the measurements coming from both kinds of sensors, to obtain 'the best of both worlds', especially for the estimation of the three-dimensional position vector. By merging these different kinds of sensors, we obtain what we call a 'mixed sensor network'.

In a realistic situation, a detection step is applied to the raw data. One then has to deal with a cluttered environment due to background noise during the tracking. The tracking may then be carried out using the probabilistic data association (PDA) algorithm [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Bar-Shalom | Tracking in a cluttered environment with probabilistic data association[END_REF], [START_REF] Bar-Shalom | Tracking and Data Association[END_REF], which takes into account the presence of false alarms (which compose the clutter), as well as a non-unity detection probability. It is then possible to estimate the parameters of interest when true detections are drowned in false alarms and do not surely appear at each sampling time. In our study, we will use the PDA with maximum likelihood estimation (ML-PDA) [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Jauffret | Trajectographie passive, observabilité et prise en compte des fausses alarmes[END_REF], [START_REF] Ciuonzo | Tracking the tracker from its passive sonar ML-PDA estimates[END_REF]. ML-PDA is a batch technique using all available measurements at once. It has been shown that this technique allows the construction of statistically efficient estimators [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Jauffret | Trajectographie passive, observabilité et prise en compte des fausses alarmes[END_REF]. Later, an extension to the ML-PDA technique using the amplitude information (AI) of each detection was created [START_REF] Kirubarajan | Passive ranging of a low observable ballistic missile in a gravitational field[END_REF], [START_REF] Kirubarajan | Low observable target motion analysis using amplitude information[END_REF]. It allows tracking in situations with very low signal-to-noise ratio (SNR).

With all these considerations, the computation of a maximum likelihood estimator to achieve the target tracking in a passive sensor network will be discussed. We place ourselves in a realistic context due to the statistical assumptions of the ML-PDA. We use both the PDA techniques mentioned above: ML-PDA and ML-PDA with AI. Our paper has 6 main sections. Section II formulates the problem as well as the notation to be used. In Section III, the observability is detailed. Section IV outlines the interest of the sensor fusion by determining the relevant asymptotic performances. Section V presents briefly the ML-PDA technique. Section VI presents some numerical results. Then follow the Conclusion, appendices, and reference list.

II. PROBLEM FORMULATION AND NOTATIONS

The mixed sensor network is composed of N S sonobuoys located on the surface of the sea and N A vertical linear arrays. At each sampling time k = 1 : K, each sonobuoy gets the TDOA of the target signal with respect to a chosen reference sonobuoy. This may be expressed as a range difference, by assuming a constant speed in the sea. Each linear array measures the cosine of the elevation angle: the angle between the direct path and the reflected path. We assume that the waves propagate in straight lines and are only reflected once [START_REF] Blanc-Benon | TMA from bearings and multipath time delays[END_REF], [START_REF] Gong | Multipath target motion analysis: Properties and implication of the multipath process[END_REF]. The reflected path corresponds to the reflection from the seabed. The positions of the sensors in the network are uniformly distributed over a square. However, each position is exactly known and the sensors are motionless during the whole acquisition time.

The target is supposed to travel at a constant velocity (CV) and at a constant depth. Consequently, the state vector defining the target trajectory is given by

X T = (x; y; z; ẋ; ẏ) (1)
where (x; y; z) is the initial position (given in meters) in a 3dimensional Cartesian space and ( ẋ; ẏ) is the velocity of the target (given in meters per second) in the horizontal plane. So, we may define the target position at each sampling time k = 1 : K with respect to X:

  x k y k z k   =   1 0 0 k∆t 0 0 1 0 0 k∆t 0 0 1 0 0   X (2)
with ∆t being the sampling period. Note that z k = z. Now we can specify the measurements acquired by each sensor. Each sonobuoy is located at the coordinates (x S i ; y S i ; z S i ), where z S i = 0. r i,k is the range between the target and the sonobuoy i at sampling time k:

r i,k = x k -x S i 2 + y k -y S i 2 + z 2 (3)
Then we get

∆r i,k = r i,k -r i0,k + ε i,k (4)
∆r i,k is the measurement coming from the sonobuoy i at sampling time k. i 0 denotes the reference sonobuoy 1 , and ε S i,k is a zero-mean additive white Gaussian noise of standard deviation σ ∆ . The measurements collected by the vertical arrays are the following:

c D n,k = cos φ D n,k + ε D n,k = z -z A n (x k -x A n ) 2 + (y k -y A n ) 2 + (z -z A n ) 2 + ε D n,k (5)
c B n,k = cos φ B n,k + ε B n,k = 2z B -(z A n +z) (x k -x A n ) 2 +(y k -y A n ) 2 +(z+z A n -2z B ) 2 +ε B n,k (6)
Equation ( 5) defines the measurement coming from the direct path, whereas [START_REF] Lo | Broadband passive acoustic technique for target motion parameter estimation[END_REF] refers to the reflected one. Both are relative to antenna n at sampling time k. Each antenna is located at the coordinates (x A n ; y A n ; z A n ), z A n being the acoustic 1 We arbitrarily choose the nearest sonobuoy from the center of the network as the reference. are zero-mean additive white Gaussian noises of standard deviation σ A . Figure 1 illustrates the elevation angles of the target. Note that the propagation delay is neglected in this study (this assumption was made in [START_REF] Arnold | Target parameter estimation using measurements acquired with a small number of sensors[END_REF], [START_REF] Lo | Broadband passive acoustic technique for target motion parameter estimation[END_REF]).

In the rest of the paper, we will consider a network composed of N S = 3 sonobuoys and N A = 2 vertical linear arrays. The locations of the sensors are uniformly distributed in a square with sides of length 20 km. The number of sampling times is K = 100 and the sampling period is ∆t = 4 s. We will consider two cases of the target velocity: a slow target linked to the true state X L and a fast target linked to X H . The true states are the following:

X T L = (-5000; 3000; -300; 4; 3)

X T H = (-5000; 3000; -300; 25; 12)

The standard deviation of the measurements coming from the sonobuoys is σ ∆ = 30 meters, corresponding to a trade-off between pessimistic [START_REF] Arnold | Target parameter estimation using measurements acquired with a small number of sensors[END_REF] and optimistic [START_REF] Moose | Adaptive underwater target tracking using passive multipath time-delay measurements[END_REF] cases. The standard deviation of the measurements collected by the antennas is σ A = 0.017 (implying a standard deviation of the angle greater than 1 degree 2 ) The acoustic center of the linear arrays z A n is supposed to be 50 m deep. The depth of the seabed z B is assumed to be -2000 m.

Fig. 2 depicts an example of the sensor network as well as the target track, for the fast target defined above. 2 For a measurement m = cos(α) + εm, the corresponding standard

deviation of α is approximatively σ(α) = σ(m) | sin(α)| .

III. ANALYSIS OF OBSERVABILITY

The following analysis will be conducted step-by-step. First, we suppose there are only two vertical antennas. For the sake of simplicity, the two antennas are located at 0; y A ; z A T and 0; -y A ; z A T . For convenience, we define A 1 = (0; y A ) and A 2 = (0; -y A ). The noise-free measurements collected by them, i.e. cos Φ D

1,k , cos Φ D 2,k and cos Φ B 1,k , cos Φ B 2,k
, correspond to the state vector X T = (x, y, z, ẋ, ẏ), and the one defining the trajectory of a ghost target, say X T G = (-x, y, z, -ẋ, ẏ). The trajectory is then observable up to the axial symmetry around the axis (0; y). Now, we add two sonobuoys allowing us to measure the time differences of arrival (TDOA). We define S 1 = (x S 1 ; y S 1 ) and S 2 = (x S 2 ; y S 2 ). The question is: do these new measurements eliminate the ghost? To answer this question, we adopt the approach proposed in [START_REF] Pérez | Observability in target motion analysis from the sum or the difference of ranges with two stationary sensors[END_REF], i.e. the measurements are taken in continuous time:

r 1 (t) = x(t) -x S 1 2 + y(t) -y S 1 2 + z 2 (9) r 2 (t) = x(t) -x S 2 2 + y(t) -y S 2 2 + z 2 (10) r G,1 (t) = x(t) + x S 1 2 + y(t) -y S 1 2 + z 2 (11) r G,2 (t) = x(t) + x S 2 2 + y(t) -y S 2 2 + z 2 (12)
Suppose that r

1 (t) -r 2 (t) = r G,1 (t) -r G,2 (t). Then [r 1 (t) -r 2 (t)] 2 = [r G,1 (t) -r G,2 (t)] 2 ⇔ r 2 1 (t) + r 2 2 (t) -2r 1 (t)r 2 (t) = r 2 G,1 (t) + r 2 G,2 (t) -2r G,1 (t)r G,2 (t) (13)
⇔ 2 x S 1 + x S 2 x(t) = r G,1 (t)r G,2 (t) -r 1 (t)r 2 (t) (14)
We deduce from [START_REF] Gong | Multipath target motion analysis: Properties and implication of the multipath process[END_REF] two remarks:

• The difference r G,1 (t)r G,2 (t)-r 1 (t)r 2 (t) is a polynomial function. • Clearly, r 1 (t)r 2 (t) is a polynomial function iff r G,1 (t)r G,2 (
t) is a polynomial function. Our analysis is based on this remark. We examine if the existence of a ghost target is compatible with the different assumptions discussed in the following. Discussion:

A. x S 1 + x S 2 = 0 Multiplying the terms of (14) by r G,1 (t)r G,2 (t)+r 1 (t)r 2 (t) we get [r G,1 (t)r G,2 (t) + r 1 (t)r 2 (t)] x(t) = [r G,1 (t)r G,2 (t)] 2 -[r 1 (t)r 2 (t)] 2 2 x S 1 + x S 2 (15) Hence, [r G,1 (t)r G,2 (t) + r 1 (t)r 2 (t)] x(t) is a polynomial.
1) ẋ = 0, x = 0, i.e. x(t) = x = 0: From ( 14), we have that r G,1 (t)r G,2 (t) -r 1 (t)r 2 (t) is a polynomial, and r G,1 (t)r G,2 (t) + r 1 (t)r 2 (t) as well by [START_REF] Moose | Adaptive underwater target tracking using passive multipath time-delay measurements[END_REF]. This implies that r 1 (t)r 2 (t) and r G,1 (t)r G,2 (t) are two polynomials of degree 2. We define now three other polynomials (of degree 2):

P (t) r 2 1 (t) Q(t) r 2 2 (t) R(t) r 1 (t)r 2 (t) (16)
We have R 2 = P Q, which is equivalent to R 2 P Q = 1. Hence P and Q divide R. Since the degree of these three polynomials is equal to 2, there is a non-zero number, say α, such that

P = αR and Q = 1 α R. We deduce that P = α 2 Q. The coefficient of t 2 is equal to ẏ2 in P and in Q. Hence α 2 = 1. Consequently, P = Q. Now the equality P -Q = 0 is equivalent to 2y -y S 1 -y S 2 + 2t ẏ y S 1 -y S 2 + 2x -x S 1 -x S 2 + 2t ẋ x S 1 -x S 2 = 0 ⇔      y S 1 = y S 2 and x S 1 = x S 2 or x = x S 1 +x S 2 2 (17)
The case where x S 1 = x S 2 and y S 1 = y S 2 must be discarded since the sonobuoys are not in the same place. Hence, we retain y S 1 = y S 2 and

x = x S 1 + x S 2 2 (18)
Now, if we define three new polynomials

P G (t) r 2 G,1 (t) Q G (t) r 2 G,2 (t) R G (t) r G,1 (t)r G,2 (t) (19)
and if we carry out similar computations, we obtain

x = - x S 1 + x S 2 2 (20)
Equalities ( 18) and ( 20) are incompatible. We conclude that when x S 1 + x S 2 = 0, the trajectory of the source is observable.

2) ẋ = 0, i.e. x(t) = x + t ẋ: First, from (14), we deduce that r G,1 (t)r G,2 (t)-r 1 (t)r 2 (t) is a polynomial. Let us rewrite x(t) as x(t) = (t -t * ) ẋ, with t * = -x ẋ (15) ⇔ 2(t -t * ) x S 1 + x S 2 [r G,1 (t)r G,2 (t) + r 1 (t)r 2 (t)] ẋ = [r G,1 (t)r G,2 (t)] 2 -[r 1 (t)r 2 (t)] 2 (21)
The time t * is hence a zero of the polynomial

[r G,1 (t)r G,2 (t)] 2 -[r 1 (t)r 2 (t)] 2 .
It can be written as

[r G,1 (t)r G,2 (t)] 2 -[r 1 (t)r 2 (t)] 2 = (t -t * )S(t), where S(t) is a polynomial. Equality (21) gives us r G,1 (t)r G,2 (t) + r 1 (t)r 2 (t) = 1 2 ẋ x S 1 + x S 2 S(t) (22)
We deduce from ( 22) that r G,1 (t)r G,2 (t) + r 1 (t)r 2 (t) is again a polynomial. Hence r G,1 (t)r G,2 (t) and r 1 (t)r 2 (t) are both polynomials.

The same reasoning as in 1) yields again the equalities P = Q, which is equivalent to

(x-x S 1 ) 2 +(y-y S 1 ) 2 +z 2 = (x-x S 2 ) 2 +(y-y S 2 ) 2 +z 2 2(x-x S 1 ) ẋ+2(y-y S 1 ) ẏ = 2(x-x S 2 ) ẋ+2(y-y S 2 ) ẏ ⇔ (2x-x S 1 -x S 2 )(x S 1 -x S 2 )+(2y-y S 1 -y S 2 )(y S 1 -y S 2 ) = 0 (x S 1 -x S 2 ) ẋ + (y S 1 -y S 2 ) ẏ = 0 (23)
We also get

P G = Q G , which is equivalent to (-2x-x S 1 -x S 2 )(x S 1 -x S 2 ) + (2y-y S 1 -y S 2 )(y S 1 -y S 2 ) = 0 -(x S 1 -x S 2 ) ẋ + (y S 1 -y S 2 ) ẏ = 0 (24)
Equations ( 23) and ( 24) can be rewritten as the following system

    2(x S 1 -x S 2 ) 2(y S 1 -y S 2 ) 0 0 -2(x S 1 -x S 2 ) 2(y S 1 -y S 2 ) 0 0 0 0 (x S 1 -x S 2 ) (y S 1 -y S 2 ) 0 0 -(x S 1 -x S 2 ) (y S 1 -y S 2 )         x y ẋ ẏ    =     (x S 1 ) 2 -(x S 2 ) 2 + (y S 1 ) 2 -(y S 2 ) 2 (x S 1 ) 2 -(x S 2 ) 2 + (y S 1 ) 2 -(y S 2 ) 2 0 0     (25)
Because the target is moving (since ẋ = 0), the determinant of the above matrix, which is 16

(x S 1 -x S 2 )(y S 1 -y S 2 )
2 , must be zero (otherwise ẋ = ẏ = 0). We have hence to consider two possibilities: either x S 1 = x S 2 or y S 1 = y S 2 (we cannot have simultaneously the two equalities, since the sonobuoys are not collocated). If

x S 1 = x S 2 , then y = y S 1 +y S 2 2
and ẏ = 0. In this case, the trajectory of the target is not observable since the two vectors X T = x,

y S 1 +y S 2 2
, z, ẋ, 0 and X T G = -x,

y S 1 +y S 2 2
, z, -ẋ, 0 give the same measurements.

If

y S 1 = y S 2 , then x = x S 1 +x 2 S 2
and ẋ = 0. This must be discarded.

B. x S 1 + x S 2 = 0 1) x S 1 = x S 2 = 0: The measurements are r 1 (t) = x 2 (t) + y(t) -y S 1 2 + z 2 = r G,1 (t) (26) r 2 (t) = x 2 (t) + y(t) -y S 2 2 + z 2 = r G,2 (t) (27)
The two state vectors X and X G give the same measurements. The trajectory of the target is unobservable.

2) x S 1 = -x S 2 = 0: The measurements are

r 1 (t) = x(t) -x S 1 2 + y(t) -y S 1 2 + z 2 (28) r 2 (t) = x(t) + x S 1 2 + y(t) -y S 2 2 + z 2 (29) r G,1 (t) = x(t) + x S 1 2 + y(t) -y S 1 2 + z 2 (30) r G,2 (t) = x(t) -x S 1 2 + y(t) -y S 2 2 + z 2 (31)
And, from ( 14), we have r G,

1 (t)r G,2 (t) -r 1 (t)r 2 (t) = 0, i.e. r1(t)r2(t) r G,1 (t) = r G,2 (t), which implies that r1(t)r2(t) r G,1 (t) 2 = r 2 G,2 (t), i.e. P Q P G = Q G (
with the definitions ( 16) and ( 19)). Note that the two polynomials

P (t) = x(t) -x S 1 2 + y(t) -y S 1 2 +z 2 and P G (t) = x(t) + x S 1 2 + y(t) -y S 1 2 +
z 2 are coprime. Hence, thanks to the Gauss Lemma, P G divides Q. Because they have the same degree, there is a real non-zero number, say β, such that P G = βQ. Consequently, P = βQ G . The coefficients of t 2 in P G and βQ being equal, we deduce that

β = 1. Consequently, r 1 (t) = r G,2 (t) and r 2 (t) = r G,1 (t), which implies that r 1 (t) -r 2 (t) = -[r G,1 (t) -r G,2 (t)].
Hence, the two state vectors X and X G do not give the same measurements. The trajectory of the target is observable. Now, we can summarize this discussion by the following theorem Theorem 1:

• If the line (S 1 , S 2 )
is not in the vertical plane containing A 1 and A 2 , then the trajectory of any target travelling in this axis is observable. Otherwise, the trajectory of any target travelling out of this axis is observable up to the axial symmetry around this line. • If the lines (A 1 , A 2 ) and (S 1 , S 2 ) are parallel, but not contained in the same vertical plane, then the trajectory of any target travelling in the perpendicular bisector of [S 1 , S 2 ] is observable up to the axial symmetry around the line (A 1 , A 2 ). • In any other cases, the trajectory of the target is observable. Consequence: If we have three sonobuoys forming a triangle, the trajectory of any target is observable.

IV. ASYMPTOTIC PERFORMANCE

In this section, we present the main concern about the fusion of the sensors in the sense of the asymptotic performance. The asymptotic performances are known from the Cramér-Rao Lower Bound (CRLB), or, similarly, from the Fisher Information Matrix (FIM). We recall the link between the CRLB for unbiased estimators and the FIM:

B = F -1 ( 32
)
where F is the FIM and B is the CRLB for unbiased estimators. First, we show the expressions for the FIM in a clean environment for each kind of sensor:

F ∆ = N S i=1 i =i0 K k=1 1 σ 2 ∆ [∇ X r i,k -∇ X r i0,k ] [∇ X r i,k -∇ X r i0,k ] T (33)
F A = N A j=1 K k=1 1 σ 2 A ∇ X cos φ D j,k ∇ T X cos φ D j,k +∇ X cos φ B j,k ∇ T X cos φ B j,k (34)
F ∆ and F A are then the FIMs for the sonobuoys and linear arrays, respectively. The FIM F after the fusion of the sensors is then

F = F ∆ + F A (35)
To illustrate the benefit of fusing the sensors, the following procedure is applied. For the slow target and the numerical values defined in Section II, 5000 random mixed sensors networks are created. Indeed, we recall that the positions of the sensors are uniformly distributed in a square with sides of length 20 km. For each network, F ∆ , F A and F are computed. Then, the empirical mean value of the diagonal elements of each FIM is extracted to allow the comparison. The result is shown in Table I.

It is noticeable that the amount of information coming from the sonobuoys on the depth of the target is very poor. But the other components would be estimated accurately. On the other hand, the linear arrays do not bring much information about the position in the plane nor on the velocity of the target, compared to the sonobuoys. But, the amount of information on the depth is rather important. Then, the mixed sensor network benefits from the best of both kinds of sensors, and thus allows an accurate estimation of the whole state vector.

Let us note that in some cases where the linear arrays are too distant, the gain of information on the depth of the target may be negligible.

V. THE PROBABILISTIC DATA ASSOCIATION MODEL

A. Statistical Assumptions of the Probabilistic Data Association Model 1) General Case: The consideration of a cluttered environment is justified by the technological roughness of the sonobuoys. These sensors cannot embed a sophisticated and complex signal processing chain. Thus, these sensors alone cannot extract the target-originated measurements.

The classic scenario is the following: a sensor gets at time k a collection s k = (s 1,k , • • • s m k ,k ) T of measurements. Some statistical assumptions are necessary to formulate an ML-PDA model. Following [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Bar-Shalom | Tracking and Data Association[END_REF], these are:

• The vectors s k are independent conditionally on X for k = 1 : K. • The detections due to the target are corrupted by a zeromean additive white Gaussian noise. The power of the noise is σ 2 (a subscript, A or ∆, will be added when applied to the mixed sensor network). These detections are called 'true detections'. • The clutter is composed of the false alarms and is distributed according to an uniform law over a scan space u.

• The number of false alarms at a specific sample k follows a Poisson law µ f a of parameter λu (the expected number of false alarms). • True detections appear at most once at each sampling time, with a detection probability P d . It is then possible to express the likelihood of X given s k , due to the total probability theorem:

L(X|s) = K k=1 1 -P d u m k µ f a (m k ) + P d u m k -1 µ f a (m k -1) m k × m k i=1 1 √ 2πσ exp - 1 2 s i,k -h k (X) σ 2 (36)
with h k the state model at time k.

2) Extension to the Use of Amplitude Information: The ML-PDA model using AI was introduced in [START_REF] Kirubarajan | Passive ranging of a low observable ballistic missile in a gravitational field[END_REF], [START_REF] Kirubarajan | Low observable target motion analysis using amplitude information[END_REF]. All the previously stated assumptions still hold. The amplitude associated with each detection is then considered. In a similar manner to the previous part and with the same notation, we may express the likelihood of X given s k :

L(X|s, R) = K k=1 1 -P d u m k µ f a (m k ) m k i=1 p τ 0 (R i,k ) + P d u m k -1 µ f a (m k -1) m k m k i=1 p τ 0 (R i,k ) × m k i=1 1 √ 2πσ ρ i,k exp - 1 2
s i,k -h k (X) σ 2 (37)
with R i,k the amplitude associated with the detection, τ the detection threshold, and

p τ 0 (R i,k ) = 1 P f a p 0 (R i,k ) (38)
p τ 1 (R i,k ) = 1 P d p 1 (R i,k ) (39)
ρ i,k = p τ 1 (R i,k ) p τ 0 (R i,k ) (40)
where p 0 and p 1 are the probability density functions (pdf) associated to the amplitudes under both hypotheses (there is only noise or there is signal and noise, respectively) and P f a is the probability of a false alarm.

B. Application to the Mixed Sensor Network

The general likelihoods under both architectures have been presented. According to the literature [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Bar-Shalom | Tracking and Data Association[END_REF], [START_REF] Kirubarajan | Passive ranging of a low observable ballistic missile in a gravitational field[END_REF], [START_REF] Kirubarajan | Low observable target motion analysis using amplitude information[END_REF], the log-likelihood function is chosen as the criterion to be maximized (up to a constant term). The noise is supposed to be independent through the samples and between sensors. The criteria will now be presented by applying the ML-PDA model to our study with the notation described in Section II:

• When the amplitude information is not available, we define the quantities

C 1 (X|∆r i,k ) = log 1 -P d∆ + P d∆ λ ∆ × m i,k j=1 1 √ 2πσ ∆ exp - 1 2 ∆r i,j,k -r i,k +r i0,k σ ∆ 2    (41) C 2 (X|c D l,k ) = log 1 -P dA + P dA λ A × m l,k j=1 1 √ 2πσ A exp   - 1 2   c D l,j,k -cos φ D l,k σ A   2         (42) C 3 (X|c B n,k ) = log 1 -P dA + P dA λ A × m n,k j=1 1 √ 2πσ A exp   - 1 2   c B n,j,k -cos φ B n,k σ A   2         (43)
where σ 2 is the power of the noise, λ is the expected number of false alarms per unit area, and P d is the detection probability. The subscript ∆ indicates the sonobuoys, and the subscript A indicates the linear arrays.

The criterion to optimize is then

C(X) = N S i=1 i =i0 K k=1 C 1 (X|∆r i,k ) + N A l=1 K k=1 C 2 (X|c D l,k ) + N A n=1 K k=1 C 3 (X|c B n,k ) (44)
• Similarily, when the amplitude information is available, we define

C 1 (X|∆r i,k , R i,k ) = log 1 -P d∆ + P d∆ λ ∆ × m i,k j=1 ρ i,j,k √ 2πσ ∆ exp - 1 2 ∆r i,j,k -r i,k +r i0,k σ ∆ 2    (45) C 2 (X|c D l,k , R D l,k ) = log 1 -P dA + P dA λ A × m l,k j=1 ρ l,j,k √ 2πσ A exp   - 1 2   c D l,j,k -cos φ D l,k σ A   2         (46) C 3 (X|c B n,k , R B n,k ) = log 1 -P dA + P dA λ A × m n,k j=1 ρ n,j,k √ 2πσ A exp   - 1 2   c B n,j,k -cos φ B n,k σ A   2         (47)
The criterion to optimize is then

C(X) = N S i=1 i =i0 K k=1 C 1 (X|∆r i,k , R i,k ) + N A l=1 K k=1 C 2 (X|c D l,k , R D l,k ) + N A n=1 K k=1 C 3 (X|c B n,k , R B n,k ) (48)
Fig. 3 presents an example of the measurements obtained from each sensor.

Remark: Following [START_REF] Lo | Broadband passive acoustic technique for target motion parameter estimation[END_REF], we neglect the correlation between the TDOA measurements.

C. Numerical Computation and Optimization

All the methods explained hereafter are applied whatever the architecture (ML-PDA or ML-PDA with AI) is.

The maximization of the criterion is carried out by using a quasi-Newton (variable metric) method, as given in [START_REF] Kirubarajan | Low observable target motion analysis using amplitude information[END_REF]. The aim is to compute an approximation of the inverse of the Hessian matrix. The update is done with the Davidon-Fletcher-Powell [START_REF] Lenard | Practical convergence conditions for the Davidon-Fletcher-Powell method[END_REF]- [START_REF] Davidon | Variable metric method for minimization[END_REF] method. The matrix is updated after each iteration of the algorithm. This guarantees that the matrix is symmetric and positive definite. The state is then updated employing the classical Newton-Raphson method:

Xi = Xi-1 -p i H -1 i ∇C( Xi ) (49)
where

H -1 i
is the approximation of the inverse of the Hessian matrix at step i and p i is the optimal step-size scaling. Note that H -1 i i converges to the inverse of the Hessian matrix. The log-likelihood function being a highly non-convex and multi-modal function [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Kirubarajan | Low observable target motion analysis using amplitude information[END_REF], several techniques are used to perform the maximization of the criterion. Indeed, due to the multi-modality aspect, the algorithm could converge to a local maximum. To avoid this, a multi-pass approach was performed in [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF] and will be now described.

This technique achieves the maximization in M steps instead of one. The algorithm is first initialized by a rough grid search method over the state space. At step m, instead of maximizing C(X), the maximization is applied to C (m) (X), which is equal to C(X) defined either by ( 44) or (48) where the standard deviations are multiplied by a m (see [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], Section III). The numerical series a m is decreasing such that a 1 = M , a 2 = M -1, ..., a M = 1. The resulting state after maximization at step m is used as the initialization at step m + 1. This method allows the algorithm to have a better 'sight' of the true detections, and ensures convergence towards the global maximum.

Moreover, to avoid convergence to an unrealistic state, a constraint (penalty function) φ(X) is added to C (m) (X) [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF]. This penalty function is about the velocity of the target and is defined as follows:

φ(X) = - 1 2 v(X) - v σ v 2 (50) where v(X) = X(4) 2 + X(5) 2 (51)
v is a typical value of the target velocity and σ v is the standard deviation around v. At step m = M , the penalty function is relaxed.

D. The Fisher Information Matrix

An expression for the FIM in a clean environment is given in (35). In the presence of clutter, the computation of the FIM will use the statistical assumptions of the ML-PDA model. First, the measurement space is restricted to a gate around the true measurements [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Bar-Shalom | Tracking and Data Association[END_REF]. The size of the gate 3 is denoted by v g : 3 For a one-dimensional measurement.

v g = 2gσ (52)
The restriction to a gate around the true measurements allows the inclusion of only the most informative measurements. Indeed, with g = 5, more than 99% of a Gaussian population is within 5 standard deviations of the mean. The measurements located outside the gate contain no additional information about the state vector. By using these assumptions as well as the likelihood computed only with the measurements inside the gate, it is possible to find an expression for the FIM in the presence of clutter. Then, the latter appears to be the FIM in a clean environment multiplied by a scalar factor: the information reduction factor. This factor, called q 2 in the literature on the ML-PDA architecture and q 3 in the ML-PDA with AI architecture, lies between 0 and 1 and represents the loss of information due to the presence of clutter and the non-unity detection probability [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Bar-Shalom | Tracking and Data Association[END_REF]. The information reduction factor only depends on P d and λv g , the expected number of false alarms in the gate of size v g . It is then clear that the loss of information will increase the CRLB. Here are expressions for q 2 and q 3 :

q 2 (λv g , P d ) = ∞ m=1 2P d √ 2πg m-1 µ f a (m -1) × g 0 • • • g 0 exp(-ξ 2 1 )ξ 2 1 (1-P d ) √ 2πλvg 2gP d + m i=1 exp -1 2 ξ 2 i dξ 1 • • •dξ m (53)
q 3 (λv g , P d ) = ∞ m=1 2P d √ 2πg m-1 µ f a (m -1) × g 0 • • • g 0 ∞ τ • • • ∞ τ m j=1 p τ 0 (R j )ρ 2 1 exp -ξ 2 1 ξ 2 1 (1-P d ) √ 2πλvg 2gP d + m j=1 ρ j exp -1 2 ξ 2 j × dR 1 • • •dR m dξ 1 • • •dξ m (54)
The derivation of q 2 is detailed in [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF] and in [START_REF] Jauffret | Trajectographie passive, observabilité et prise en compte des fausses alarmes[END_REF], where it is generalized to a multi-dimensional measurement. The derivation of the q 3 is given in Appendix A, based on [START_REF] Kirubarajan | Passive ranging of a low observable ballistic missile in a gravitational field[END_REF], [START_REF] Kirubarajan | Low observable target motion analysis using amplitude information[END_REF].

Then, it is possible to express the FIM in the presence of clutter in both architectures:

F 1 = q 2 (λ ∆ 2gσ ∆ , P d∆ )F ∆ + q 2 (λ A 2gσ A , P dA )F A (55) F 2 = q 3 (λ ∆ 2gσ ∆ , P d∆ )F ∆ + q 3 (λ A 2gσ A , P dA )F A (56)
where F 1 and F 2 are the FIM in ML-PDA and ML-PDA with AI architectures respectively.

These factors are computed numerically by using Monte Carlo runs [START_REF] Robert | L'analyse statistique bayésienne[END_REF].

E. Test for Track Acceptance

Once the algorithm returns an estimate, it is necessary to decide if the estimated track is correct or not. Indeed, whether there is a target or not, the algorithm computes an estimate. If there is no correct detection, the estimate must be rejected as well as when the algorithm converges towards a local maximum. So we have to choose between two hypotheses:

• H 0 : There is a track and X is the global maximum of the likelihood function. • H 1 : There is no track, or X is a local maximum. The natural reflex would be to use the so-called generalized likelihood ratio (GLR) to construct a test between these two hypotheses. This needs the statistical distributions of the GLR (or a monotonic function of it) under each hypothesis. Determining these distributions is a tough task (from [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], Section IV, "the statistical distribution of [the GLR] is not available"). Following [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], we propose the palliative described below to circumvent this. We give the details for the case where the amplitude information is not available. We consider the following statistic 4 :

T 0,1 ( X) = T 1 + T 2 + T 3 √ V 1 + V 2 + V 3 (57)
with

T 1 = N S i=1 i =i0 K k=1 C 1 X|∆r * i,k -E C 1 X|∆r * i,k |H 0 (58)
T 2 = N A l=1 K k=1 C 2 X|c D, * l,k -E C 2 X|c D, * l,k |H 0 (59)
T 3 = N A n=1 K k=1 C 3 X|c B, * l,k -E C 3 X|c B, * l,k |H 0 (60) V 1 = N S i=1 i =i0 K k=1 Var C 1 X|∆r * i,k |H 0 (61) V 2 = N A l=1 K k=1 Var C 2 X|c D, * l,k |H 0 (62) V 3 = N A l=1 K k=1 Var C 3 X|c B, * l,k |H 0 ( 63
)
Since X is supposed to be the global maximum, we get the following approximations.

E C 1 X|∆r * i,k |H 0 E C 1 X|∆r * i,k |H 0 = µ 0,1 (64) E C 2 X|c D, * l,k |H 0 E C 2 X|c D, * l,k |H 0 = µ 0,2 (65)
E C 3 X|c B, * l,k |H 0 E C 3 X|c B, * l,k |H 0 = µ 0,3 (66)
4 * is for the measurements inside the gate.

Var C 1 X|∆r * i,k |H 0 Var C 1 X|∆r * i,k |H 0 = σ 2 0,1 (67) Var C 2 X|c D, * l,k |H 0 Var C 2 X|c D, * l,k |H 0 = σ 2 0,2 (68) Var C 3 X|c B, * l,k |H 0 Var C 3 X|c B, * l,k |H 0 = σ 2 0,3 (69)
The statistic (57) can be rewritten as

T 0,1 ( X) = C( X) -(N S -1)Kµ 0,1 -N A K(µ 0,2 + µ 0,3 ) (N S -1)Kσ 2 0,1 + N A K(σ 2 0,2 + σ 2 0,3 ) (70
) It has been shown empirically that T 0,1 follows the standard Gaussian distribution by using the Central Limit Theorem [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF]. With all these considerations, the test is the following: if T 0,1 > c α H 0 , is accepted. Otherwise, H 0 is rejected. c α is then the threshold with a level of significance α.

When the amplitude information is available, only the values of µ 0,i and σ 0,i will change. The expressions for µ 0,i and σ 0,i are detailed in [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Jauffret | Trajectographie passive, observabilité et prise en compte des fausses alarmes[END_REF] for the ML-PDA architecture and in Appendix B for the ML-PDA with AI architecture. low. The detection threshold is placed such that P f a = 0.02 and P d = 0.6. The value of K is still 100. For this configuration, 100 Monte Carlo runs are performed. The result of these runs is shown in Tables IV andV. The 95% confidence interval of the NEES for 100 simulations is [4.37; 5.63]. The mean value of the NEES is in the confidence interval and near the theoretical value of 5. Hence, the estimator is efficient.

VII. CONCLUSION

In this paper, we treated the problem of estimating the trajectory of a source in CV motion when the available measurements are a set of TDOA provided by a sonobuoy network, and the cosines of elevation angles acquired by vertical antennas. To begin with, the observability was studied before proposing an estimate. Then, we chose a convenient model for noisy measurements. Due to the simplicity of the system embarked in a sonobuoy (which performs very basic signal processing functions such as thresholding the signal energy), the data were assumed to be composed of correct detections and false alarms. Consequently, we employed the PDA model, well adapted to this situation. The proposed estimator (the MLE) was found to be efficient (the empirical covariance matrix reaches the CRLB) even with only three sonobuoys and two vertical linear arrays. Coupling their energy to these thresholded data improve the accuracy of the MLE, which remains efficient in this case. Nevertheless, two points were neglected: the propagation delay of the waves carrying the signal, and the correlation between the TDOA measurements. Future research will be focused on these two points.

APPENDIX A DERIVATION OF THE INFORMATION REDUCTION FACTOR q 3

We derive q 3 according to the principles detailed in [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Jauffret | Trajectographie passive, observabilité et prise en compte des fausses alarmes[END_REF], [START_REF] Kirubarajan | Passive ranging of a low observable ballistic missile in a gravitational field[END_REF], [START_REF] Kirubarajan | Low observable target motion analysis using amplitude information[END_REF]. We first express the likelihood at time k by considering only the measurements inside the gate of size v g . We use the same notation as in Section III.

L(X|s) = 1 -P d u m µ f a (m) m i=1 p τ 0 (R i ) + P d u m-1 µ f a (m -1) m m i=1 p τ 0 (R i ) × m i=1 1 √ 2πσ ρ i exp - 1 2 s i -h(X) σ 2 (71)
After factorization we get

L(X|s) = µ f a (m) v m g m i=1 p τ 0 (R i ) 1 -P d + P d λ × m i=1 1 √ 2πσ ρ i exp - 1 2 s i -h(X) σ 2 (72)
We define

Φ(X|s) = 1 -P d + P d λ m i=1 ρ i √ 2πσ exp - 1 2 s i -h(X) σ 2 (
73) Then, the information reduction factor q 3 is

q 3 = ∞ m=1 1 2πσ 2 P d λ 2 µ f a (m) v m g × g -g • • • g -g ∞ τ • • • ∞ τ m i=1 p τ 0 (R i ) Φ(X|s) × m i=1 m j=1 ρ i ρ j exp - 1 2 s i -h(X) σ 2 - 1 2 s j -h(X) σ 2 s i -h(X) σ s j -h(X) σ × dR 1 • • • dR m ds 1 • • • ds m ( 74
)
Jérémy Payan was born in Toulon, France in 1997. He received his masters degree "Vision Signal Trajectographie Automatique" from the Université de Toulon, France in 2019.

He is currently pursuing a Ph.D. degree in signal processing in the Université de Toulon, France, with the french laboratory IM2NP, and with the company Naval Group Research, Ollioules, France. His work is about target motion analysis methods for passive sonar systems.

Antoine Lebon was born in Briey, France in 1994. He received his masters degree "signal et trajectographie" from the Université de Toulon, France in 2018.

He is currently pursuing a Ph.D. degree in signal processing in the Université de Toulon, France, with the french laboratory IM2NP, and with the company Naval Group Research, Ollioules, France. His work is about target motion analysis methods for passive sonar systems.

Fig. 1 .

 1 Fig. 1. Elevations of the direct path and of the reflected path.

Fig. 2 .

 2 Fig. 2. Example of a mixed-sensor network. Dots: sonobuoys. Vertical lines: linear arrays. Stars: Acoustic center of linear array. Arrow: Target.

Fig. 3 .

 3 Fig. 3. Measurements obtained for each sensor in one run. P d = 0.6, λu = 8.

  Annie-Claude Pérez was born in France on November 10, 1965, received the Diplôme d'Etudes Appronfondies in Optics and Image Processing, from the Université de Toulon, France, in 1988, and the title of "Docteur de l'Université" in 1991 from the Université de Toulon, France.Since Sept. 1994, she has been at the University de Toulon, where she teaches electronic systems. Her research is focused on signal processing, applied to biomedical systems before turning them to target motion analysis. Claude Jauffret was born in France on March 29, 1957, and received the Diplôme d'Etudes Appronfondies in Applied Mathematics, from Saint Charles University, Marseille, France, in 1981, the Diplôme d'Ingénieur from Ecole Nationale Supérieure d'Informatique et de Mathématiques Appliqués de Grenoble, Grenoble, France, in 1983, the title of "Docteur de l'Université" in 1993, and the "Habilitation à Diriger des Recherches" from the Université de Toulon, France. From November 1983 to November 1988, he worked on passive sonar systems, more precisely, on target motion analysis at GERDSM, France. After a sabbatical year at the University of Connecticut (from November 1988 to December 1989), during which he worked on tracking problems in a cluttered environment, he developed research in tracking, data fusion, and extraction in CERDSM. Since Sept. 1996, he has been at the Université de Toulon, where he teaches statistical signal processing. His current research is about observability and estimation in nonlinear systems, as it appears in tracking problems. Dann Laneuville graduated from "Ecole Supérieure d'Electricité" (a.k.a. Supélec) in 1987 and got a Ph.D. in automatic and signal processing from the Paris XI Orsay University in 1998. He is currently R& D engineer in the Advanced Algorithms and Architectures department within Naval Group Research, Naval Group's research and technology centre based in Nantes, France. Before joining Naval Group in 2003, he has been for ten years within the Image and Signal Processing Laboratory of Matra Group and with MBDA's department of guidance and control for five years. Before joining Naval Group Research at its creation in 2012, he has spent a year as a visiting scholar within the Electrical and Computing Engineering department at the University of Connecticut where he studied maneuvering target tracking and video extraction algorithms. He has published over forty papers in international renowned conferences and workshops (Conference on Decision and Control, Aerospace Conference, International Conference on Information Fusion, IET seminar, SDF workshop, MAST, UDT. . . ), been a key actor during fifteen years within the "collège de Polytechnique" multisensor fusion seminar. His areas of interest are multitarget tracking, maneuvering target tracking, bearing only tracking and data fusion.

TABLE II RESULTS

 II OF 200 MONTE CARLO RUNS WITH ML-PDA ARCHITECTURE FOR THE SLOW TARGET (X = X L )

	(P d ; λu)	(0.8; 2) (0.6; 8)
	x		-5004	-5005.5
	ȳ		2998.9	2995.6
	z	-298.09	-298.55
	x		4.036	4.013
	ȳ		3.015	3.010
	NEES		5.39	5.37
	Acceptance rate	94.5%	84%
		TABLE III
	RESULTS OF 200 MONTE CARLO RUNS WITH ML-PDA ARCHITECTURE
	FOR THE FAST TARGET (X = X H )
	(P d ; λu)	(0.9; 2) (0.8; 6) (0.6; 8)
	x	-5001.7	-5004.4	-4918.9
	ȳ	2999.1	2999.3	3013.8
	z	-299.77	-300.69	-298.03
	x	24.993	25.011	24.936
	ȳ	11.999	12.002	11.958
	NEES	5.18		5.13	5.36
	Acceptance rate	96.5%	94%	93%

TABLE IV RESULTS

 IV OF 100 MONTE CARLO RUNS WITH ML-PDA WITH AI ARCHITECTURE FOR THE SLOW TARGET (X = X L )

	X	NEES Acceptance rate
	(-5005.9; 2995.5; -303.39; 4.069; 3.043)	5.48	99%
	TABLE V
	RESULTS OF 100 MONTE CARLO RUNS WITH ML-PDA WITH AI
	ARCHITECTURE FOR THE FAST TARGET (X = X H )
	X	NEES Acceptance rate
	(-4999.2; 3002.7; -300.19; 24.997; 11.990)	5.01	99%

VI. NUMERICAL RESULTS

In this section, some numerical results of the estimator are presented.

A. Results in ML-PDA Model

We first present the numerical results obtained for the ML-PDA architecture in different configurations. The following assumptions are made:

• P d = P d∆ = P dA • λ ∆ and λ A are set such that the expected number of false alarms λu is the same for each sensor. Each configuration is defined by specific values of (P d ; λu) and by the specific positions of the sensors (according to the uniform law in a square with sides of 20 km). We recall that the number of sampling times is K = 100. Then, 200 Monte Carlo runs are performed for each configuration. The results of the Monte Carlo runs are presented in Tables II andIII.

X is the mean of the estimates. The mean value of the NEES 5 [START_REF] Bar-Shalom | Tracking and Data Association[END_REF] is computed. For 200 simulations, and a 5-dimensional parameter, the confidence interval of the NEES at 95% is [4.55; 5.45]. In each configuration, the NEES is close to the theoretical value [START_REF] Van Trees | Detection, Estimation, and Modulation Theory, Part I[END_REF] and belongs to the confidence region. This indicates that the estimator is efficient, that is, unbiased and meeting the CRLB.

B. Results in ML-PDA with AI Model

As we use the AI, we have to define the pdf of the amplitude under both hypotheses. It is assumed that the energy of the background noise follows a central χ 2 2 distribution, and the target originated signal follows a non-central χ 2 2 distribution. The non-centrality parameter depends on the SNR, which is 5 Normalized Estimation Error Square.

We define

Reintroducing in (74) these notations, we obtain

Then, we use the following principles to simplify the expression for q 3 : it is an odd function with respect to ξ i for ξ j fixed. The cross terms then vanish from the integrals. Next, we use the parity of the function with respect to ξ i . After some manipulations, we finally obtain the expression given in (54).

APPENDIX B DERIVATION OF THE MOMENTS OF THE LOG-LIKELIHOOD FUNCTION UNDER TARGET HYPOTHESIS IN ML-PDA WITH AI We derive µ

(n) 0 in the ML-PDA architecture. The derivation is based upon [START_REF] Kirubarajan | Passive ranging of a low observable ballistic missile in a gravitational field[END_REF], [START_REF] Kirubarajan | Low observable target motion analysis using amplitude information[END_REF]. The criterion is the following:

We want to compute (80) with the likelihood given in (71):

For m = 0 we have

We use the same change of variables given in (75), ( 76), (77) and we get

We then use the parity with respect to ξ i and notice that the integrals are the same for each considered index. After some manipulations we get