Position Determination of Acoustic Source Using a Discrete Probability Density Method Supported by Joint Detection and Time of Arrival Estimation

The localization of various acoustic sources in a battlefield (such as weapon rounds, mortars, rockets, mines, improvised explosive devices, vehicle-borne improvised explosive devices and airborne vehicles) nowadays has significant history. In acoustic source localization systems multiple sensors (such as microphones or microphone arrays), placed at known positions and networked, are used to detect signals emitted from the source and perform localization of the source. In this paper, Discrete Probability Density (DPD) method, as a method for position determination, has been used to estimate location of acoustic sources such as artillery weapons. Wavelet decomposition of the acoustic signal enables to emphasize N-wave in the time domain. Time of Arrival (TOA) estimation is realized by a statistical analysis of narrow time segment which hold Nwave. Based on the estimated TOA Time Difference of Arrival (TDOA) is calculated. TDOA presents relative time difference of arrival of the sound signal between pairs of sensors. Localization of the acoustic sources is performed by time difference estimation of acoustic signal arrivals using Discrete Probability Density method. A number of experiments with fire off mortar projectiles have been conducted in order to verify the performance of the proposed approach.


Introduction
HE localization of various acoustic sources in a battlefield (weapon rounds, mortars, rockets, mines, improvised explosive devices, vehicle-borne improvised explosive devices and airborne vehicles) nowadays has significant history. In the open literature there are many papers in which methods for acoustic source localization are presented [1][2][3]. Some of them are based on cross-correlation techniques and some on near-field beam-forming that can provide good localization accuracy even in noisy and reverberant environments [1][2][3][4][5].
Generally, the passive methods for determining emitter positions are divided into two major groups: one-step positioning techniques (centralized processing) and two-step positioning techniques (decentralized processing). One-step techniques or direct methods are based on direct position determination. Acquired samples from all sensors are sent to a common central unit which jointly processes the received data and estimates the emitter location in a single step. Two-step positioning techniques or indirect methods are based on the estimation of a specified parameter of signals such as the Angle of Arrival (AOA) or the Time of Arrival (TOA) at each sensor. In the second step, the estimated parameters are sent to the central sensor in order to determine the emitter positions.
All of these methods utilize the sound recordings of the sensors to calculate the most probable source location.
However, data transfer from the sensors to the processing unit requires relatively high bandwidth, which can be a serious limitation for practical implementation. This is especially case if one-step techniques are used for position determination. Another possible bottleneck in such systems is the huge amount of calculations to be performed and limiting the response time [5].
In this paper we propose localization method that utilizes a distributed two-step approach. The proposed method for acoustic source localization using Discrete Probability Density (DPD) method is based on using distributed two-step approach. The first step of the proposed method is performed at the sensors where TOA is determined (certain points of the acoustic signals are marked). These TOAs from all the sensors are recorded and sent to the central processing unit. In the second step the central processing unit calculates Time Difference of Arrival (TDOA) between sensors to obtain the position of the signal source using the DPD method.
In order to perform DPD method for position determination an important step is a reliable detection of explosion, which implies precise estimate of the moment or TOA estimation when this event occurs. Detection of the large explosions in the near field is relatively trivial task. A particular problem is detection of small explosions in presence of background noise at large distances. Background noise is mostly a consequence of various human activities. In addition, local turbulences in air have influence on transmission paths of the sound thus making the problem more complex. In many cases there is a real need to make a distinction between different types of T explosions, where a physical acoustic signature of explosion sources manifests such differences and this kind of approach introduces another problem -identification, which is not a topic of this paper.
In this paper, particular attention is dedicated to detection of N-wave, which generates artillery projectile along its trajectory after fire off. N-wave is a characteristic pressure wave or shock wave. Related problem with detection of Nwave is localization of the source of N-wave and that is, commonly, to perform by space distributed and networked acoustic sensors placed in region where N-wave will probable occur. The localization techniques, which are applicable for large ranges, utilize phase shifts of the acoustic signals. Generally, in literature many different techniques for localization of the acoustic sources are described, such as cross-correlation technique, beamforming, and many variations of the localization techniques based on the estimation of the time of arrival, or TOA.
Optimal and wavelet-based N-wave detection and estimation method is presented in the paper Sadler et al. [6]. Significant attention is focused on gradient-based detection and estimation schema, as a classical means of estimating changes of the signal amplitude. Gradient method exploits very fast rise and fall times of the N-wave edges, and therefore is effective and robust with respect to unknown environmental interference.
Detection procedures of the signals with fast changing amplitudes, which are based on changes of statistical distributions such as a step change of the mean are relatively frequently used, see Basseville and Benveniste [7]. The binary hypothesis test -likelihood ratio test of two hypotheses, H1: change of the mean, against the hypothesis H0: no change of the mean is used for detection of the step changes of the signal amplitudes.
The method for TOA estimation is based on a statistic analysis of the N-wave amplitudes of data segment, which hold N-wave. Estimation of the data segment, which holds Nwave, is done by the wavelet decomposition of the acoustic signal. Estimation of the moment, which belongs to N-wave, is performed by the analysis of an appropriate wavelet coefficients of the acoustic signal in time and amplitude domain. This method enables to estimate TOA of an N-wave with the acceptable accuracy.
In purpose to perform a verification of the proposed method for TOA estimation, precise measurement of the source position was necessary to do. That was done by a commonly positioning technique based on GPS. Such approach enables verification of the quality of proposed method for estimation TOA of N-wave by comparing positions of the mortar obtained by two different methods. The first is very accurate GPS method and the second is DPD method based on the estimated TDOAs by proposed method for TOA estimation.
This paper consists of six parts. Introduction is given in Section I. The method for joint detection and TOA estimation is presented in Section II. The concept of DPD method in combination with TDOA measurement is given in Section III. Proposed method for acoustic source localization is given in Section IV. Some results are presented in Section V. Conclusions are given in Section VI.

Joint Detection and TOA estimation method
Moving projectile with supersonic speed, V>c (Mach number is M=V/c>1, quotient of projectile velocity V and speed of sound in air is c), generate acoustic monopoles at its surface which produce acoustic pressure cone of a constant phase in direction which with the projectile trajectory close the angle M θ , Eq. (1), see Pierce [8], Acoustic pressure, p, of the front of N-wave is: where ρ is density of the air, F W is Whitham function (depends on the projectile shape) and r is radial distance from trajectory of the projectile. When the shape of the projectile is deduced on a projectile diameter, d, and projectile length, l, and taking into account parameter, y, perpendicular distance from the projectile trajectory to the sensor, or miss distance, last equation is possible to write in something different form, see Sadler et al. [6], where Δp is the pressure jump at the start of N-wave and p 0 is ambient pressure of the atmosphere. Estimation of the TOA of an N-wave in acoustic signal is resolved by the algorithm supported by wavelet analysis. The problem of choosing the level of smoothing and type of the wavelet function is resolved in the research phase by a strategy, which implies employing a multiscale wavelet decomposition of the acoustic signal and improving signal to noise ratio by denoising process. Such strategy is accomplished with Daubechies wavelet coefficient of the seventh order, or "db7" wavelet, see Daubechies [9]. Ten-step schema is used for TOA estimation. Rejection of the false peaks is performed on de-noised signal afterwards by taking into account only the first peak in the series of closely separated peaks.
As an illustration how a proposed new algorithm for TOA estimation works its flow chart diagram is shown in Fig.1. Wavelet decomposition of the acoustic signal realizes denoising of the acoustic signal and increases probability to emphasize position of the N-wave in the acoustic signal. Relatively complex analysis of the amplitude peaks of "db7" wavelet coefficients enables estimation of the moment where there is a maximum probability that the narrow time segment around estimated moment hold N-wave. In the flow chart diagram, see Fig.1, is shown how the algorithm for TOA estimation works.
TOA estimation based on described method is not timeconsuming process. Processor time of standard PC computer with relatively obsolete Intel Core 2 Duo, 2.93 GHz processor, in case when in signal, whose length is 30 s and sampling rate 8333 Hz, and where only one N-wave signal exists, is of order hundred milliseconds. Characteristic shape of the acoustic signal with an N-wave is shown in Fig.2 upper diagram, and shape of the sameN-wave with more details is shown at the lower diagram of Fig.2. N-wave has characteristic shape and its appearance is accompanied with muzzle waves which have multipath spread from the source through the sensor position and therefore have certain delays but muzzle wave is still within the echo signal and cannot be clearly resolved, see the lower diagram in Fig.2. Joint detection of the acoustic signal is performed by five at terrain of interest placed and networked sensors, or microphones. Acquisition of the acoustic data is done simultaneously with synchronizing error less or equal 20 ns. One acoustic sensor was taken as referent for calculation time difference of TOA, or TDOA. As example one concrete measurement, in Fig.3 is shown estimated TOA value and Nwave obtained on the base of the acoustic signal which was measured at the first sensor's position.
The results of TOA estimation for four other sensor's positions are shown in diagrams, see Fig.4. All of these TOA values of an N-wave are obtained in an automatic manner without an operator in the loop. It is necessary to point out that the implementation of this algorithm implies the existence of a relatively good signal to noise ratio. Natural noise interferences such as thunder in near surroundings of the area of interest have no significant impact on power of the algorithm to resolve N-waves and estimate its TOA values (see Fig.5). Significant natural interference is possible to notice in diagram, see Fig.5, where few thunder events are registered between 9 th and 13 th second.

Discrete Probability Density method for position determination
The DPD method is a technique developed at Defense R&D Canada for numerically combining data containing uncertainty such as Electronic Warfare (EW) sensor measurements [10]. The DPD method is based on the assumption that the sensor measurement of a given parameter, in this case for example AOA or TOA, can be modeled by some probability density function (PDF) over the range of possible values. These PDFs may be represented by various distributions such as Gaussian distribution or von Mises PDF. The DPD method combines the probability density distribution of measurements directly by sampling each PDF at common intervals and calculating the joint product over the space. The DPD method is applied to determine the spatial locations by projecting the measurement PDF from each sensor, or in the case of TDOA method combination of two sensors onto a common grid of sample points. This requires that some transform function exists in order to map the measured parameter into the 2-dimensional space (2-dimensional grid in X, Y of the Cartesian coordinate system). The two-dimensional TDOA DPD array can be represented by: where F XY is the TDOA DPD array of size N×M and f(θ ij ) is Gaussian PDF for a pair of sensors i and j. In this case TDOA PDF can be represented by a Gaussian function by: where µ is the mean value of tdoa, and estimated error is the standard deviation σ.
For multiple DOA measurements, the joint DPD array is calculated over a common N×M grid

Thunder
Results obtained using DPD method based on the measurement of TDOA is presented in Figures 6 and 7. It is assumed that there are four sensors and one target at the position (7.411.656,84, 4.966.355,68). Given position is at Gauss-Krieger coordinates. In Fig.6 are shown F XY PDFs for three estimated TDOAs. Estimated target position is at that point that represents an intersection of these three estimated PDFs. In Fig.7 is shown an estimation of the target position P XY using (7) with spatial resolution of 125 meters.

Proposed method for acoustic source localization
In order to provide geo-localization of an acoustic source in the battlefield an approach has been proposed in this paper. The proposed approach is based on two-step techniques for geolocalization using sensor network. In Fig.8 the setups of three sensors for position determination of a mortar are shown.
The block diagram of the proposed approach is shown in Fig.9. After choosing sensors positions, based on the 3D model of terrain for area of interest (AOI), the first step is to calculate distances from all possible target positions to each sensor, and to save that data in the memory. 3D terrain model of AOI, with dimensions 40km × 40km, is generated using geographical information system (GIS) software. In this software, it is possible to select a desired AOI, to define dimensions of AOI and to define the spatial resolution. Based on these input data, GIS software generates corresponding 3D terrain model of AOI.
By processing 3D terrain model of AOI, useful data are extracted and as a result are obtained matrix with terrain heights and 2-dimensional grid in X, Y of the Cartesian coordinate system. Possible target positions in 3D terrain model of AOI are defined by 2-dimensional grid in X, Y of the Cartesian coordinate system. Distances from all possible target positions to each sensor are calculated based on obtained matrix with terrain heights and 2-dimensional grid in X, Y of the Cartesian coordinate system. For example, to calculate distance from one possible target position to one sensor, a terrain profile is observed in such a way that terrain heights between these two points are taken into account. All calculated distances are placed in matrix. This matrix is called the referent matrix and is used during the proposed approach to obtain matrices with distances from new possible target positions to each sensor. New possible target positions are defined with better spatial resolution which means with new 2-dimensional grid in X, Y of the Cartesian coordinate system in the same AOI (see Fig.10). Matrix modifications are done using interpolation. Interpolation is used for distances estimations, from new possible target positions to each sensor, based on distances from the referent matrix.
The number of possible target positions depends on the 2dimensional grid in X, Y of the Cartesian coordinate system used for DPD method. For example, if the dimension of AOI is 40km × 40km (in this AOI are included also positions of the sensors), depending on the processing power and accuracy, we can choose different spatial resolutions. For spatial resolution of 100m, there are 160.000 possible target positions. The distances from all possible target positions, the number of positions depend of the dimension of AOI, to all sensors are calculated. Distances are determined based on 3D model of terrain. After preliminary estimation of target position, it is possible to focus only on the part of AOI close to the preliminary estimation of the target position and to perform DPD method with better spatial resolution in order to achieve better accuracy. Proposed method for acoustic source localization is based on more iteration, K (see Fig.9).
Calculated distances are important because based on these distances and other important parameters such as temperature and wind flow, TOAs are calculated from all possible target positions to the different sensors. Based on these TOAs and chosen reference sensors TDOAs can be calculated for given sensors setup and grid of possible target positions.
The next step in the proposed approach is estimation of TOAs at each sensor in the network. TOAs can be estimated manually or automatically using the proposed method joint detection and TOA estimation method.
If there is a possibility to transfer all data from each sensor to central unit, TOAs can also be determined at the central unit.

Results of the field experiment
Validation of the proposed method for acoustic source localization is done during three field experiments in different time period. Reference position of the acoustic source was determined by GPS and distances from sensors to acoustic source position were precise measured. Meteorological data necessary for sound speed estimation, such as: temperature, humidity, wind velocity and wind direction have been gathered during the field experiments by the appropriate sensors. Joint acquisition of the acoustic signals from all microphones was synchronized with the accuracy of 20 ns. For all gathered acoustic signals detection and TOA estimation is performed using the proposed method.
At the beginning of the process of acoustic source localization by DPD method it was necessary to establish so called AOI with some linear dimensions. For an implementation of DPD method, 3D terrain model of AOI is used. In this case, AOI with dimensions of 20×20km has been set. DPD method for source localization is implemented as an iterative algorithm and seven iterations were necessary to perform localization of the mortar. During the first step, AOI was represented as grid 40×40, with spatial resolution of 500 m. After the initial -rough estimation of the mortar position, new AOI is formed around previously estimated position of the acoustic source with same grid dimensions of 40×40, but with improved spatial resolution. The final spatial resolution was 11 m during the implementation of this localization step.
In the first field experiment the source of N-waves was fire off 120 mm mortar projectiles. Six networked acoustic sensors were placed at different locations at the area of the field experiment. Measured distances between mortar position and six sensors were 472, 828, 1547, 2090, 2310 and 285 meters, respectively.
The first sensor is used as a referent 1 sensor for TDOAs estimation. There were 34 recorded acoustic signals during this experiment. In Figures 11 and 12    In the third field experiment the source of N-waves was fire off 82 mm mortar projectiles. Four networked acoustic sensors were placed at different locations at the area of the field experiment. Measured distances between mortar position and four sensors were 996, 486, 760 and 1585 meters, respectively. The second sensor is used as a referent sensor for TDOAs estimation. There were 29 recorded acoustic signals during this experiment.
In Figures 15 and 16     Accuracy of the mortar position, which is obtained by DPD method, is also validated by DGPS preliminary estimated position of the mortar. In Fig.17 RMS errors of the localization estimation for each of the three field experiments are shown. RMS errors are calculated for x-coordinate, y-coordinate of the mortar position and joint x-y position of the artillery weapons. Generally, on RMS errors of the artillery weapons position estimate by DPD method two main sources have influence. The first one is a quality of TOAs, or TDOAs, estimation, and the second is a consequence of the uncertainty in the estimate of the sound speed and transmission path at AOI, due to the fact that atmospheric meteorological parameters, which have influence on sound speed and sound transmission, are not exactly known in sufficient number of points through the whole sound trajectory. Networked acoustic sensors cover relatively great area, but meteorological parameters are measured only in one position, in this case. Joint RMS errors for all fire off in the first field experiment in case of joint (x-y) estimate is 7.96 m, and for x-coordinates and y-coordinates mean RMS errors are 4.87 m and 6.30 m, respectively. Joint RMS errors for all fire off in the second field experiment in case of joint (x-y) estimate is 9.06 m, and for x-coordinates and y-coordinates mean RMS errors are 4.66 m and 7.77 m, respectively. Joint RMS errors for all fire off in the third field experiment in case of joint (x-y) estimate is 7.08 m, and for x-coordinates and y-coordinates mean RMS errors are 4.96 m and 5.04 m, respectively.

Conclusion
In this paper, we proposed an approach for the acoustic source localization using joint detection and time of arrival estimation combining with Discrete Probability Density (DPD) method, as a method for position determination. The proposed approach is based on an estimation of time of arrival of acoustic signal at each sensor and in the second step applying DPD source localization is determined. This proposed method does not require high communication bandwidth for data transfer between sensors.
The proposed approach has been used to estimate location of the acoustic sources such as artillery weapons. A number of experiments with fire off mortar projectiles have been conducted in order to verify the performance of the proposed approach. Based on the results from these field experiments it can be concluded that the proposed approach for acoustic source localization can be applied and implemented in system for acoustic source localization. According [11] CEP of a few meters has been observed under very favorable conditions. CEP is more typically of 100 m under 8 km range, and 1%-2% of range at longer range. Based on that, it can be concluded that achieved CEP by the proposed method in these three field experiments satisfied.