About Times of Arrival Estimation for Sources and Buried Objects Localization

: Sources localization is an important theme of antenna and signal processing researches since decades. It has a lot of applications for source characterization, which leads to object detection, used in submarine ﬁeld as well as seismic and remote sensing studies. The particularity of array processing is that we are looking at spatio-temporal vector signals which are generated by both radiating sources or reﬂected by buried objects and ambient medium. These signals are generally corrupted by an additive noise. When the number of sources is greater than the number of the sensors of the smart antenna system, classical methods become unusable. Recently, a new subspace method called incoherent sensor by sensor processing (ISSP) for estimation of directions of arrival (DOAs) is proposed. It is based on a study sensor by sensor to estimate ﬁrstly times of arrival(TOAs), then the DOAs of the sources are estimated. The main drawback of this method is the computational load needed to estimate accurately the DOAs. In this paper, we propose new faster methods thanks to the numerical complexity reduction of the existing algorithms. This last point is important nowadays to make embedded systems faster.


Introduction
Detection and localization of sources, in the underwater acoustic environment, have received considerable attention.Many studies have been recently developed.In order to perform better analysis, a lot of fields rely on it and link it to decision theory.Applications of localization can be found in telecommunications, submarine detections, astronomy, RADAR, medical imaging and diagnostic, building maintenance, earthquake analysis.Systems have to be increasingly accurate and accelerated with using fewer resources.They have to be more autonomous too, as we can see today for embedded systems, such as drones, whose development is growing steadily.Since the half of the last century, array processing research arouses interest all around the world.Many high resolution (HR) methods are developed, and one of the most interesting is the MUSIC algorithm [1,2].It uses a subspace-based method in order to perform eigendecomposition and identify both signal subspace and noise subspace.It was firstly used for frequency estimation then used for objects and sources localization [3].ESPRIT algorithm [4] is also a method based on subspace separation, but use rotational invariance of the smart antenna system which is linear and uniform.These methods have been studied to go further [5], get more results [6].They are constantly improved [7].In these HR methods, incoherent methods are suitable for wide-band sources localzation [8].The HR methods are, generally, developed under the main assumption which assumes that the number of the radiating sources is small compared to the number of the sensors of the antenna.Nowadays, in many cases, our aims are to detect several objects in our field of view.But on the other hand, our specifications will not allow us to exceed a certain sensor network length.Moreover, the manufacturing cost increases rapidly.Therefore, it would be useful to develop methods without constraint on the sensor array length.In order to provide an answer of this problematic, a new subspace-based method, Incoherent sensor by sensor processing (ISSP), was developed [9,10].This method brings a new point of view : instead of using the number of sensor, which leads to a spatial diversity; the process uses the frequential content of each sensor on the network.By using a subspace decomposition, it is possible to get an estimation of delays between each sensor or times of arrival (TOA).In [11], it was proposed to use this estimation to localize the emitting sources.For this, it was necessary to link TOA of each sensor to each source.In [9] was presented two algorithms.The first one proposed to match each TOA by looking sensor by sensor and source by source, but its computation time is too heavy.The second one reduced the complexity by using an a priori first DOA.The purpose of this paper is to review the first algorithm and to present new hypothetical ways to compute TOAs without any knowledge on DOA, with a complexity gain.We will begin to explain the signal model then we will move on the analysis of the TOA association algorithm.Otherwise we will introduce both thought which could reduce the computation time.
The remainder of the paper is organized as follows : Section 2 presents the data model.Section 3 gives the TOA estimation by HR methods.Section 4 presents the detailed description of our proposed for the association of estimated TOAs.Some algorithms and comparative results are contained in Section 5. Finally, section 6 concludes the paper.Within the scope of this paper, scalar is denoted by x, vector by x, matrix by X, estimation of A by Â, the transposition by (.) T , the transconjugation by H and M N,P (R) for the real matrix space of size (N × P ), N rows and P columns, C can also be used.

Signal Model
We begin to consider P sources emitting a signal which impinge a linear antenna system composed of N sensors.These sources are characterized by their direction to the array main sensor θ p where p = 1, ..., P .Figure 1 illustrates the situation.At a source p is associated the temporal signal s(t).The ambient environment acts on it as an amplitude modulation c p,n , which depend of both sensor and source and assumed to be independent of time.The whole event is experienced with an additive Gaussian noise b(t).So, for sensor n = 1, ..., N , we receive such signal r n (t) : τ p,n represents the time delay of arrival (TOA) from the source to the sensor.It is also equal to : With ρ p,n the euclidean distance source-sensor and v the celerity in the current medium.Whatever, we can write the Fourier transform as the signal will be sampled in order to get its Discrete Fourier Transform.So we have : Then we use M frequencies (m = 1, ..., M ) to sample r n and : We compute all sampling frequencies into one vector : where we have these relations : • S = diag(s(f 1 ), ...s(f M )) the diagonal matrix of original signal at each working frequency, • c n = [c 1,n , ..., c P,n ] T amplitude coefficients vector, • A n ∈ M M,P (C) the transfer matrix of our system, • b n , independent Gaussian background noise, A n is the transfer matrix "sources to sensors per working frequency".It contains all phases terms : Instead of using spatial notions, as DOA estimation with HR algorithm [3], we will use frequency notions.These notions were already used in hybrid method called MUSICAL [12] which aimed to compute DOA and distances array-source.Method developed in [9] suggest to use an estimation of TOA sensor by sensor with HR methods, then move forward to sources localization by using them to calculate DOA.

TOA estimation by HR methods
Therefore, it is necessary to use signals statistics to apply HR algorithms on sensor n.It needs to write the covariance matrix Γ n ∈ M M,M (C) of this sensor : That leads us to : The following text is written assuming that background noise is uncorrelated and white.So, Γ noise = σI M .
To take the analysis one step further, it can be possible to do frequential smoothing, whitening and look at sub-bands influence [11].Under hypothesis of full-rank matrix, Γ n has M eigenvalues γ m , m = 1, ..., M and also M eigenvectors v m , m = 1, ..., M .Assuming that the number of sources P is known, we can sort the P largest eigenvalues and their eigenvectors concatenated in a matrix V = [v 1 , ..., v P ].Let Π b be the projector onto noise orthogonal subspace.So Π b = I M − VV H , with normalised eigenvectors.The MUSIC cost function can be re-arranged to : Where a(τ p ) is a column vector of A n .Estimated TOA for sensor n is found by maximization of this function.By moving on each sensor from 1 to N we get a set of TOA called (τ p,n ) n=1,...,N p=1,...,P .Accuracy of TOA estimation calculation on nth sensor can be found in [11,13,14].By writing τ p = τ p + ∆τ p , it is possible to get a relation between ∆τ p and ∆Γ n .By computing HR methods on each sensor, we get all TOAs τ p,n and now it is necessary to link each p family in regards to the N sensors.

Times of arrival association: Overview of ISSP algorithms
We are working under this following assumption : • Far field: waves are supposed to be plane waves so we are dealing with the fact that sources are far from the antenna.According to [9], this assumption involves a linear relation between each sensor.With d the spacing inter-sensors, v the celerity of the ambient environment and θ 1 the first DOA, we have the phase shift : Let define a matrix P rows (for P sources) and N columns (for sensors), to have a representation of the TOAs : TOA matrix with P = 3 and N = 4.
Without knowledge on TOAs links, it is a priori impossible to say if, for a TOA on sensor "n", the next TOA on sensor "n + 1" will be really linked, even if they are very close from each other.Maybe that is true, maybe both TOA refer two strictly different sources.
It needs to substitute the index "p" by "q n " (for sensor n) to bring out the fact that we got a completely random mix of times of arrival.How can we know which sequences {τ qn,n | n ∈ 1, 2, ..., N ; q n ∈ 1, 2, ..., P } are referencing real wavefronts/real signals ?
Figure 2: The TOA association processing, without information, can be difficult.Here is illustrated an example with P = 4 sources and N sensors.
In [9,11], was presented two main algorithms which performs TOA association.The first one consists on working on a two-by-two TOA comparisons possible for the (q n , n) which leads to q n = p for sensor n.It uses the minimization of the following comparison criterion : where We give a depiction of this algorithm :

Algorithm 1 Two-by-two comparison algorithm
Require: T OAs, N ∈ N, P ∈ N for q 1 from 1 to P with step 1 do for all combinations of (q 2 , ..., q N ) ∈ [1, ..., P ] N −1 do Compute the f criterion end for Take the sequence (q 2 , ..., q N ) which minimizes f Associate it with source q 1 end for Ensure: Give all combinations.
This algorithm performs the TOA association as it was expected in [9].Nevertheless, its time complexity is about P N .Its computation time increases rapidly with P even if a small number of sensors is used.It was decided to build a predictive algorithm.This second algorithm works with an information on a first DOA θ q1,1 and a first combination (τ q1,1 , τq2,2 ).The DOA is estimated in [−π, π].In the end, this algorithm works with a lower time complexity than the first one.It is equal to P 2 N 2 N θ where N θ is the number of values used to calculate θq1,1 .This method is interesting but needs one more estimation for the first direction of arrival.What we wanted to do was to build a new algorithm which is free of any DOA and have a lower complexity than the two-by-two comparison.In the following sections we will present both methods to expand.

Proposed methods: Matching criterion
Looking the TOA matrix as we have seen in (11), written below, gives a priori no relation between each terms composing it.
Indeed a sensor's sort, therefore a column's sort is not a good way to associate TOA n with its neighbor.It is not impossible that two close times have no relation.Under considered hypothesis, let us make a TOA association criterion.By equation ( 10), working with N sensors (1, ..., n, ..., N ) and P sources (1, ..., p, ..., P ) (or (1, ..., q n , ..., P ) for the unknown sequence), we have : We apply it between the n and n + 1 (or n − 1) sensors with reference to the q 1 = p term : But we keep a dependence of DOA θ q1,1 .In order to avoid it, we define the ratio : Notice the fact that we suppose to find the right τ for sensor n.Therefore, if relation ( 17) is satisfied for a q n+1 TOA : we must find the right time of arrival for sensor n + 1.This indicator, or its equivalent |R p qn+1 − 1 n |, which we need to minimize, seems to highlight a relation between TOA.Except for the second column of the TOA matrix, where we have (n = 1) : for all values of q 2 .This issue can be resolved in a low complexity P 2 by working on the three first column of the TOA matrix with a scanning for τ p,1 on all values of columns 2 and 3. We define an algorithm and apply it to a reference time τ p,1 .In the following, States is an array storing the matched TOAs.
index of TOA q 1 = p q 2 = 3 q 3 = 2 index of Sensor 1 2 3  Require: for t from 1 to P with step 1 do for q from 1 to P with step 1 do Calculate (τ (t,3)−τ (q,2) Return States; Ensure: Give the three first states for source p. q 1 = p q 2 = q r q 3 = t r 1 2 3 Table 2: States array at the output of the first matching (P=3 and N=4) With this setting, we can compute the criterion for the others TOAs in a new algorithm.In the worst case, its time complexity is equal to P (P 2 + N P ).The concatenation of the first one, this one and the last presented below does the complete expected work for one source p.We will have to repeat it P times : Algorithm 3 TOA association columns 3...N for n from 1 to N − 1 with step 1 do 1 → index; tcurrent → qcurrent; for q from 1 to P with step 1 do Calculate Return States; Ensure: Give missing states for source p.
From the States array we get the right TOAs with : Algorithm 4 States to TOA Require: N ∈ N, τ ∈ Mn,p(R), States for i from 1 to N with step 1 do τ (States(i, 1), i) → τ List; end for Return τ List; Ensure: Give TOAs for starting source p. q 1 = p q 2 = q r q 3 = t r q current = 2 1 2 3 4 τ List = {τ 1,p , τ qr,p , τ tr,p , τ qcurrent,p } Table 3: States array at the end of these 3 algorithms Onto the paper, our criterion R to associate TOA seems to be intuitive, but maybe not instinctive.The idea of developing a new algorithm free of any more estimation than TOAs is present and now need to be consolidated.It needs to be tested with data.

Evaluation of the proposed criterion
As we tried to build a new algorithm an more exactly a new criterion, it needed to be tested.In this section the simulated data are obtained using the following values of the different parameters : • Number of sensors of the antenna: N = 5.
• Sources are arranged randomly.The coordinates of the P sources are in R 2 reported to a right-handed Cartesian coordinate system centered on the sensor 1.The only constraint is to put them with the condition y ≥ 0, dealing with the fact that we will look ahead the antenna.
Let remind the delays on sensor n : If a source is put at (x p , y p ) coordinates and a sensor at (x n , y n ), then we have : In order to simulate a mix q n on sensors, the TOA matrix is blended column by column.In the following text we will show an ordered matrix.Let resume the situation(figure 4) : In this case, we get a noiseless TOA matrix (truncated) : By applying our algorithm we must do the right TOAs association and find as same as the matrix (figure 8) : For 100 supervised tests, we got these results : • Success mean = 0.94 • Computation time mean = 0.012399 s, with plots.
For 100×1000 tests, we got : For a large number of realizations, the success rate value is high.

Estimation of R criterion
In this section we want to study the influence of propagation error.In order to do this, we introduce an "error quantity" e p,n for each τ p,n .Therefore we can write again the criterion R p qn+1 , noticed R p qn+1 + ∆R If p = q n+1 , we should have |R p p − 1 n | = 0.With errors, we get : The problem is : what happens if it exists one q n+1 which is not the right value but leading to The TOAs association concludes on an error for this case.

Conclusion
In this paper, we developed a new method for ISSP in order to build a fast and accurate TOA association algorithm without DOA estimation.We also aim to reduce numerical complexity, compared to the original ISSP algorithm.The numerous obtained results showed that the proposed method leads to good associations, even if a 100% success rate would be appreciated.Further developments would be to introduce this criterion in a Viterbi algorithm [15] and hidden Markov model.In this paper, this work was done under the far field assumption, allowing us to write a linear relation between TOA, its extension to the near-field would be interesting.

Figure 4 :
Figure 4: Sensors are colored in magenta and sources in light green (N = 5 and P = 8)

Figure 5 :
Figure 5: The TOAs association shows their linear relations.

Figure 8 :
Figure 8: R indicator when a large number of trials(615)is used.