Position Determination of Acoustic Source Using Higher-Order Spectral Analysis for 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 a significant history. In acoustic source localization systems of multiple sensors, networked and placed at known positions are used to detect signals emitted from the source and perform localization of the source. Time of arrival estimation of the gun fire shock waves at large distances is significant when the purpose is the gun localization. At large distances shock wave signal is significantly deformed and therefore it is relatively difficult to accurately estimate the time of arrival. In this paper methods for TOA estimation are proposed, which are based on cumulants of the acoustic signal which originate from distant gunfire events. Localization of acoustic sources is performed in two steps by time difference estimation of acoustic signal arrival and using Discrete Probability Density method for positioning. Numerous field experiments have been conducted in order to verify performance of the proposed approach.


Introduction
COUSTIC signal which originates from a gunfire event is a stochastic non-linear process which generally more or less deviates from the Gaussian process of a known mean. Non-linearity is a consequence of explosive adiabatic process, which is generated by the acoustic signal. Therefore, if the goal is to get the answer whenever explosive event occurred, it is necessary to obtain the information regarding the deviations from the Gaussian assumption and the presence of nonlinearity in the acoustic signal. The presence of characteristic nonlinearity in acoustic signal is possible to use as an indicator of the existence of explosive events. Classical means of data analysis, such as power spectrum where the process is considered as a superposition of statistically uncorrelated harmonic components and suspended phase information is of little help in this case. Additional means of statistical analysis of the acoustic data is necessary to detect possible nonlinearities in the signal.
Large explosions detection in the near field is relatively trivial task. Particular problem is the detection of small explosions in the presence of background noise at large distances. Background noise is mainly consequence of various human activities. In addition, local turbulences in the air have influence on transmission paths of the sound and thus making that problem more complex. In many cases, it is necessary to distinguish different types of explosions, where a physical acoustic signature of explosion sources manifests such differences and this kind of approach introduces another problem -identification, which is not the topic of this paper. time of arrival (TOA), of the acoustic signal of the interest, and to perform weapon localization. In this paper, it is used the localization method proposed in [3], based on the Discrete Probability Density (DPD) method.
Estimation of the TOA in situation of a BoN-wave in acoustic signal is resolved by algorithm, which is developed in the Military Technical Institute [3]. Starting point in the development of the algorithm was knowledge that the main contribution to the acoustic signal gives characteristic shock, BoN-wave and muzzle blast wave, B-wave, which occurs during shooting from the weapon. There are not acoustic sources with similar characteristics in normal circumstances. The shock and muzzle blast waves, which are generated by projectiles of different calibers, are difficult to mask and therefore it is possible to utilize these waves for localization of near placed small arms, or distant placed cannons. The first step of the algorithm is recognition of such characteristic acoustic events. This is performed by filtration of the raw acoustic data in order to remove noise and to improve signal to noise ratio (SNR). The best results were performed by wavelet decomposition, see papers [3,4].
In this paper, another approach in the algorithm for TOA estimation was made by using the fourth order cumulant in order to identify time interval when BoN-wave occurs, because BoN-wave produces nonlinearity of the acoustic signal during its occurrence. Using cumulants is a promising method because it enables the detection of BoN-waves in acoustic signal until nonlinearities exist. From the point of view extending the range of detection and localization of these events leads to the fact that TOA of these events is possible to estimate until the coherence of the atmosphere for acoustic signal is preserved and nonlinearities produced by these events exist.
The methods based on higher order statistics are very useful in problems where non-Gaussian, non-minimum phase, phase coupling or nonlinear behavior and robustness to additive noise are important. Detection and classification using higher order statistical and spectral techniques are proposed for using in many areas where pattern recognition is important. These methods have the potential to elicit better performance from sensors and sensor networks [5].
The paper is organized as follows. Introduction is given in Section I, theory of the cumulants is given in Section II, examples of the application of cumulants are given in Section III, some practical results are presented in Section IV. Conclusions are given in Section V.

The cumulants theory
The relation between moments and cumulants is shown by the following equations which exist for orders k=1,2,3,4 for the sequence of data (n).
According to the theory, the first-order cumulant, C 1x , of a stationary process is the mean, The most significant parameter for estimation nonlinearities in signal is the fourth order cumulant, 4 ( , , ), defined as follows: where 2 ( ) = { ( ) ( + )}, and equals C 2x (m), for a real-valued process, and {(k,l,m): 0 ≤ m ≤ l ≤ k ≤ ∞} are time lags. According to citied literature the first-order cumulant, 1 ( ), is the mean of the process, and the second-order cumulant, 2 ( ), is the autocovariance sequence. The zerolag cumulants have special names: C 2x (0) is the variance, but C 3x (0,0) and C 4x (0,0,0) are usually denoted by γ 3x and γ 4x . The normalized quantities are referred as the skewness and as the kurtosis. These normalized quantities are both shift and scale invariant. If x(n) is symmetric distributed, its skewness is necessarily zero (but not vice versa); if x(n) is Gaussian distributed, its kurtosis is necessarily zero (but not vice versa). The terms skewness and kurtosis are used to refer to the un-normalized quantities, γ3x and γ 4x , see literature [6]. The use of cumulants rather than moments has a number of advantages as pointed out by J. Mendel [7]. The only available data are those related to the type of weapon from which the gun was fired. Based on the analyses performed using the fourth-order cumulants and the results obtained, as well as prior knowledge of these processes, it is possible to obtain a clearer picture of the conditions and processes that occurred during the realization of these experiments. In addition, the examples given are illustrative enough to understand the potential and importance of applying the fourth order cumulants.

Example 1
Acoustic data used in the first example contain recorded 9 mm gunshot, with sampling frequency of 44100 Hz. The fourth order cumulant of the signal is shown in Fig. 1. The arrow marked with number 1 indicates the occurrence of a shock wave, and the arrow marked with number 2 indicates the wave that occurred when ejecting the products of combustion of gunpowder from the mouth of the barrel of the weapon Gunshot and its dynamic are clearly illustrated with the fourth order cumulant lags in the time domain. It is possible to notice that gunshot signal consists of two components. The first component of the cumulant is positioned at the beginning of the process, see arrow 1, and with relatively pronounced amplitudes because its origin is shock wave with direct transmission. The second component, see arrow 2, comes somewhat later and the source of this nonlinearity in signal is occurrence of delayed muzzle wave. The amplitudes of cumulant lags which belongs to the muzzle waves, are significantly larger than amplitudes that originate from shock wave. The reason for this is that the energy of the muzzle wave is much higher than the energy of the shock wave. The energy of the shock wave depends on the speed and shape of the bullet or projectile. In the time domain cumulant peaks are somewhat expanded in the time domain due to transmission conditions of these waves. All other smaller cumulant peaks that occur, see diagram, are the result of the propagation of the shock and muzzle waves along different and longer propagation paths due to reflection on obstacles in the environment. It can be said that a bullet whose caliber is 9 mm produces relatively strong BoN-waves.

Example 2
Two short bursts of three bullets each fired from an AK 47 assault rifle were analyzed in this example. The diagram, see Fig. 2, clearly shows the acoustic disturbances resulting from these firings. When analyzing the structure of these disorders in detail, it can be concluded that they are very similar to those analyzed in detail in the previous example. The contribution of the shock wave to the total signal is less pronounced if the recording is made from a greater distance since the amplitude of the shock wave with distance decreases significantly. The dominant content of the acoustic signal is the muzzle wave, which becomes more distorted with increasing distance. Also, at greater distances, the effect on the recorded signal is greatly enhanced by the propagation along multiple wave paths due to numerous reflections. Different methods for their recognition and extraction are possible to use. The use of the data obtained by this analysis relates to the problem of solving the place from which these bullets were fired. In this sense, it is necessary to determine as accurately as possible the moment when this process began, that is, the first occurrence of the shock wave or the muzzle wave, if the shock wave is very weak and cannot be detected at that distance. In real conditions, the frequency of occurrence of these events can be very high, so it is of practical importance to realize systems and algorithms to automatically recognize and analyze them without a human being in the loop. Neural networks and signal analysis algorithms based on them can be of practical use in this case.

Examples 3 and 4
Examples 3 and 4 illustrate the application of the fourth order cumulants in the analysis of underwater explosions. Unlike acoustic phenomena in the air environment where sound recording is made with a microphone in the water environment, it is done with a hydrophone. The shock waves generated in the water when explosions occur can be lifethreatening for the diver and therefore need to be behind solid obstacles to avoid contact with the direct shock wave front. In this case, the diver or the hydrophone itself may be exposed to shock waves reflected from the water / air boundary. However, such reflected shock waves have the property of reducing the pressure at the observed location, ie the pressure amplitude is negative. During the recording of underwater sound, the position of the hydrophone in the water relative to the position of the underwater explosion may be such that the hydrophone can be reached by: a) a direct shock wave, a shock wave reflected from the water surface, and a shock wave component reflected from the bottom; b) the reflected wave from the bottom and the reflected wave from the water surface; or c) only the reflected wave from the water surface. Considering that the analysis showed that almost all the cumulative lags were negative, it can be concluded that the position of the hydrophone with respect to the position of the underwater explosion was such that it was only possible to receive shock waves reflected from the water / air boundary surface, that is, they had a negative pressure producing in the hydrophone position dilution. Reflection of the shock waves from the water/air boundary is described by a well known principle of images. The effect of the shock wave reflected from the water / air boundary surface is effectively viewed as having no boundary surface, that it has the same intensity and time change is a negative sign. These principle is applicable whenever the pressure amplitudes are small enough for the laws of ordinary acoustics to apply [9]. As in previous example, in Fig. 4. the fourth order cumulant lags are shown in the time domain. In this case it is possible to notice strong reverberation and multi path transmission effect and five reflected components.

Example 5
The diagram given in Fig. 5 shows the fourth order tsunami wave cumulant. A tsunami wave is known to occur as a result of strong seismic earthquakes at the bottom of the sea or ocean, and as a consequence a powerful surface wave of tremendous destructive power develops. The diagram shows an arrow marked with the number 1 showing the cumulants of the sound signal originating from an earthquake that arrived directly at the hydrophone without significant previous reflections from sea surface. The cumulants in the diagram indicated by arrow 2 correspond to the sound generated by the tsunami wave on the coast. Considering that the cumulants are generally negative, we can conclude, as in the two previous examples, that the sound produced by impact of the tsunami wave into the coast was transmitted to the hydrophone by being repeatedly reflected from the water surface. Based on the examples shown, it can be concluded that the introduction of the fourth-order cumulants in the analysis of characteristic acoustic, hydroacoustic and even seismic phenomena (tsunami wave) enables their better separation from ambient noise, as well as a clear recognition of their characteristic features. The analysis presented is of particular importance in the case where our goal is to determine the position of the source of these disturbances based on the recorded signals at different points in the field. In order to perform this task, it is necessary to determine with great precision the moment when these disturbances arrive at our sensors -microphones, hydrophones or geophones. In other words, the time of arrival of the wave or (ToA) should be specified. It is understood that in the process there is every single best possible extraction of these phenomena from the total ambient noise that we can consider to have a Gaussian distribution. For the automatic determination of ToA, it is necessary to develop special algorithms that take into account the specificities of these phenomena.

The field experiment results
More examples of acoustic signals, which deal with different types of explosive events, are analyzed with the fourth order cumulant. Validation of the proposed method for acoustic source localization is done during two filed experiment in different time period. Reference position of the acoustic source was determined by GPS and distances from sensors to acoustic source position were precisely measured. Meteorological data necessary for sound speed estimation, such as: temperature, humidity, wind velocity, and wind direction have been acquired during the field experiments by appropriate sensors. Joint acquisition of the acoustic signals from all sensors (microphones) was synchronized with the accuracy of 20 ns. For all acquired acoustic signals, detection and TOA estimation is performed using the proposed method. First two field experiments are performed in order to investigate behavior of the proposed method in case of near distance between sensors and acoustic source.
In the first experiment, which was performed in a field, the source of the BoN-waves was firing off 120 mm mortar shell. 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.
In the second filed experiment the source of BoN-waves was the fire off the 105 mm howitzer projectiles. Five networked acoustic sensors were placed at different locations in the area of the field experiment. Measured distances between mortar position and six sensors were 2596, 1529, 1072, 470 and 755 meters, respectively.
In this example signal from mortar projectiles was analyzed. The source of BoN-waves was the fire off the 120 mm mortar projectile. The distance between mortar position and sensor is D=2090 m. Samples of the signal recorded at position of the sensor with sampling frequency of 8333 Hz are shown in Fig. 6a. Recorded signal is analyzed with the fourth order cumulant and the results are shown in Fig. 6b In Table 1. are shown the results of TOAs estimation and determination of position for the first field experiment of one firing off. TOAs are estimated for different number of samples in the time series of 10, 100 and 250 and for 10 times lags. Based on the shown results, it can be concluded that increasing the number of samples in the time series causes worse accuracy of TOA estimation and that directly reflects the accuracy of acoustic source position determination.
Based on shown results, it can be concluded that using fourth order cumulant analysis it is possible to determine TOA of the acoustic signals. In order to calculate magnitude of the fourth order cumulant it is necessary to provide the significant usage of the memory space and the processing time. Processing time for calculation depends on the number of the samples per frame or time series. Decreasing the number of samples per frame causes increasing of the processing time. However, decreasing the number of samples per frame causes better accuracy of the TOA estimation.
For determination of TOA based on calculation of the fourth order cumulant it is necessary to locate a maximum of magnitude 4 ( , , ). The maximum is determined using well known method of interpolation or curve fitting. The usage of described first and second field experiments estimation of TOA is performed by both algorithms: algorithm for joint detection and TOA estimation method proposed in [3] (Alg.1) and the algorithm based on the fourth order cumulant (Alg.2). In Table 2. are shown obtained results for one firing off. Based on the shown results, it can be concluded that obtained results for position estimation are in both case acceptable. Time series of 10 samples and 10 time lags are used for TOAs estimation.
Accuracy of the howitzer position, which is obtained by DPD method using the fourth order cumulate analysis for TOA estimation, is also validated by GPS preliminary estimated position of the mortar. Positions of the sensors and estimated positions for 9 fires off and for TOAs estimation are used time series of 100 samples and 10 time lags are shown in Fig.8. The estimation of the position using the proposed algorithm in [3] is shown at the same Figure. RMS error is calculated for joint x-y howitzer position . Joint RMS error for 9 fire off in the first field experiment is 3.31m for the fourth order cumulant and 4.10m for algorithm proposed in [3] and for x-coordinates and y-coordinates mean RMS errors are 2.11 m and 2.55 m, respectively in the case of Alg.2 and 2.94m and 2.76m in the case of Alg.1.
In accordance to referent position of the sensor system, RMS error of azimuth is 2.74 mils and RMS error of distance is 1.40 m in the case of the Alg.2. In the case of Alg.1 RMS error of azimuth is 3.56 mils and RMS error of distance is 1.08 m. Generally, two main sources have an influence on RMS errors of the artillery weapons position estimated by DPD method. The first one is quality of TOAs, or TDOAs, estimation, and the second is a consequence of the uncertainty in estimation 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 whole sound trajectory. Networked acoustic sensors cover relatively great area, but meteorological parameters are measured only in one position, in this case.
Accuracy of the mortar position, which is obtained by DPD method using the fourth order cumulate analysis for TOA estimation, is also validated by GPS preliminary estimated position of the mortar. For TOAs estimation time series of 100 samples and 10 time lags are used. Sensors positions and the estimated positions for 34 fires off are shown in Fig.9. The same Figure shows the position estimation using the proposed algorithm in [3]. Joint RMS error for 34 firings off in the second field experiment is 2.23m in the case of fourth order cumulant and 1.92m in the case of algorithm proposed in [3] and for xcoordinates and y-coordinates mean RMS errors are 1.66 m and 1.49 m, respectively in the case of Alg.2 and 1.64m and 1.02m in the case of Alg.1.
In accordance to referent position of the sensor system RMS error of azimuth is 1.10 mils and RMS error of distance is 1.83 m in the case of Alg.2. In the case of Alg.1 RMS error of azimuth is 1.07 mils and RMS error of distance is 1.51 m.
The next two field experiments are performed in order to investigate behavior of the proposed method in the case of far distance between sensors and acoustic source [10]. In the next field trial the source of BON-waves was fire of the 130 mm cannon projectiles. The distance between mortar position and sensor is more than D=4200 m in the first case and more than D=8100 m in the second case. In Fig. 10a. time samples of the signal recorded at the position of the sensor for the distance of D=4200 m are shown, and in Fig.  11a samples of the signal recorded at the position of the sensor for the distance of D=8100 m. The fourth order cumulate analysis has been performed in both cases. The results are shown in Fig. 10a and Fig. 11b Accuracy of the cannon position, which is obtained by DPD method using the fourth order cumulate analysis for TOA estimation, is also validated by GPS preliminary estimated position of the mortar. Positions of the sensors and estimated positions for all 10 fires off for the second case at the distance of more than 8100 m are shown in Fig. 12. Joint RMS error for all fire off in the first field experiment in the case of joint (x-y) estimate is 95.05 m, and for x-coordinates and y-coordinates mean RMS errors are 75.9 m and 57 m, respectively. Errors in coordinates estimation for the all 9 fires off are shown in Fig. 13. In accordance to the referent position of the sensor system, RMS error of azimuth is 2.79 mils and RMS error of distance is 101.83 m. The presented experimental results show that the time of arrival estimation of shock wave, based on cumulants, can be used in the process of acoustic source localization. Using Higher-Order Spectral Analysis for the estimation of TOA on each sensor/node in combination with Discrete Probability Density (DPD) method for localization is possible to achieve good performance of the system for acoustic source localization. Especially in the case that shock wave is distorted due to transmission through fluctuating atmosphere at large distances.

Conclusion
In this paper, an approach for time of arrival estimation using High Order Spectral Analysis in the process of position determination of the acoustic sources was proposed. The proposed approach is based on two-step techniques for localization. The calculation of the fourth order cumulant in order to estimate TOA of acoustic signal is implemented in the first step in localization techniques.
Numerous field experiments with fire off mortar, howitzer and cannon projectiles have been conducted in order to verify the performance of proposed approach. Based on the results from these field experiments, it can be concluded that proposed approach for acoustic source localization can be applied and implemented into the system for acoustic source localization. The proposed approach exhibits almost similar characteristics as algorithm proposed in [3] especially in the case when acoustic sensor is near to acoustic source and there is no distortion due to transmission through fluctuating atmosphere. However, this approach outperforms algorithm proposed in [3] in case when there is a distortion due to transmission through fluctuating atmosphere at large distances.