Search Space Reduction for Localization and Tracking of an Acoustic Source

: Experimental data from the SACLANTCEN 1993 Mediterranean Experiment are reviewed to assess the reduction of the search space for the localization and tracking of an acoustic source in a three-dimensional environment. Key to this goal is the availability of an initial estimate of source range and depth (called the 2D initial guess); an ambiguous estimate of source bearing can be obtained from the 2D initial guess through Environmental Signal Processing, and the ambiguity can be removed by searching for the source only in the range/bearing regions where bearing estimates are higher. This search provides a new estimate of source range and a single bearing, which together with the estimate for source depth constitute the center of the reduced search space for source localization and tracking. The suggested approach is tested on experimental data from the SACLANTCEN experiment considering different frequencies, as well as a stationary and a moving source.


Introduction
Source localization is at the core of Matched Field Processing (MFP) [1], a fundamental technique of Underwater Acoustics that generalizes the concept of plane-wave beamforming [2]. In MFP, a set of candidate predictions of the acoustic field (called replicas) are to be compared with a given set of signal observations acquired at a single array; the replica with the better fit to the observations can be used to infer the position of the source. A critical condition for MFP to be successful is to have a proper knowledge of the sound speed profile, bottom properties, and local bathymetry and to take advantage of a reliable underwater acoustic model. Generally speaking, there are two different types of source localization: the first one is two-dimensional (2D) source localization, which can be seen as the problem of finding the source coordinates (r s , z s ) within the search space defined by [r min , r max ] × [z min , z max ]; the second case is three-dimensional (3D) source localization, which corresponds to the problem of finding the source coordinates (r s , z s , φ s ) within the search space defined by [r min , r max ] × [z min , z max ] × [φ min , φ max ], where φ stands for source bearing, and in general φ max − φ min = 360 • (the choice of cylindrical coordinates is a typical convention). In both cases, the calculation of replicas takes place among all possible combinations of coordinates within the search space; for obvious reasons, the number of combinations is much higher in 3D source localization than in 2D source localization. Strictly speaking, 2D environments do not exist because ocean bathymetries are always 3D; however, there are often configurations of the source and the array in which the sea bottom transect can be approximated as flat. In such cases, one can replace 3D source localization with 2D source localization and rely on a normal mode model to provide an efficient calculation of replicas over a given interval of source depths; this replacement is attractive in order to reduce the size of the search space but is achieved at the cost of losing the source bearing because the calculation of replicas depends on range and depth only. Two-dimensional MFP-based localization had been discussed in detail considering simulations and experimental data for both the narrowband and broadband cases [3][4][5][6][7][8]; additionally, non-MFP 2D localization with a single hydrophone has been shown to be possible as well [9,10]. Environmental Signal Processing (ESP) was proposed to overcome the bearing limitations of 2D source localization in [11] by suggesting to take advantage of the bathymetry. To this end, one can consider two different approaches: the first one is the so-called N × 2D modeling, in which replicas are to be calculated along 2D transects over N bearings; in the second approach, replicas are to be calculated taking into account all the variability of the 3D environment. Environmental Signal Processing based on the N × 2D approach using the KRAKEN normal mode model [12] was discussed in detail for experimental data in [13], but only for a subset of all possible bearings; under such conditions, the particular bounds on the source bearing certainly allow us to reduce the size of the search space, but this is a specific type of a priori information that cannot be expected to be always available. Current computational resources certainly allow us to handle 3D source localization successfully for all combinations of parameters of the search space. However, by proceeding in this manner, the method seems to be missing one important fact, namely that the calculation of replicas over the entire 3D search space is undesirable because of the corresponding computational burden; besides, such calculation is also inefficient because a large percentage of the replicas will provide a low fit to the observations. The approach described in this discussion looks forward to the reduction of the search space for 3D localization through the identification of reliable bearings; to this end, source localization starts with a "2D initial guess", which is nothing more than the result of 2D localization obtained under the assumption that the waveguide is flat. The 2D initial guess is then combined with ESP to identify regions of reliable bearings, and these reliable bearings are further reduced to a single bearing by searching over the bearing/range space. This search provides a unique bearing and a new range, which, together with the source depth of the 2D initial guess, represent the center of a reduced neighborhood in which the 3D search can take place. Certainly, there are cases, like cross-slope propagation in a wedge environment [14], for which it is very unlikely to obtain a reliable 2D initial guess. However, the diversity of scenarios for which the 2D initial guess can be obtained is too broad to be considered here; it would be necessary to compare vertical versus horizontal arrays, narrowband versus broadband processing, mild versus rough bathymetries, single versus multiple hydrophone arrays, etc. This complex issue deserves to be addressed in future independent discussions. The proposed approach is exemplified here using the exceptional data from the SACLANTCEN 1993 Mediterranean Experiment, which are available online. (http://spib.linse.ufsc.br/sonar.html accessed on: 14 June 2023) [15]. The experiment itself is described in brief in Section 2, together with the corresponding experimental data and the 2D initial guess. Localization and tracking are discussed in Section 3, while conclusions and future work are presented in Section 4.

The Experiment
The SACLANTCEN 1993 Mediterranean Experiment is described in detail in the literature [16]; thus, a compact description is presented in this section. Transmissions took place in a shallow water area north of the island of Elba in an environment with a mild bathymetry and over a two day period on 26 and 27 October 1993. A vertical array (VLA) with 48 hydrophones spanning the water column from 18.7 to 112.7 m was deployed on the morning of the 26th, and on the same day, a source was deployed at a stationary position 5.8 km north of the VLA; on the 27th starting at 14:00 h, a support ship moving northward towed a source with similar characteristics to the one deployed at the stationary position. The stationary source was attached to a surface buoy, which was itself tethered to a ballast lying on the bottom; the length of the tether cable was such that the position of the stationary source was only known to be within a circle of radius 200 m with its center at the ballast. Depending on wind conditions, the stationary source could be closer to or further from the VLA within that circle. GPS coordinates for the VLA and the source ballast are indicated in Table 1; the corresponding accuracy for each position was around 100 m. Two signals, called RM2 and RM5, were transmitted at different moments by the stationary source, while the moving source transmitted only the RM5 signal; both RM2 and RM5 were transmissions of pseudorandom noise with bit shift registers for which there is no available information. According to the literature, center frequencies for RM2 and RM5 corresponded to 335 and 170 Hz, respectively [16]; however, frequency-based optimization indicated that a better fit between replicas and data could be achieved using slightly different values of frequency (see Section 2.2).

The Experimental Data
The received signals are distributed in MATLAB mat files, each containing around 1 min of transmissions with a sampling rate of 1 kHz; every file stores the time series in matrix form, with every row containing the data of one of the 48 hydrophones. The data were recorded on 26 and 27 October 1993 and after additional processing were organized in the following way: • Five files for the stationary source, RM2 signal, recorded on 26 October. • Ten files for the stationary source, RM5 signal, recorded on 26 October. • Ten files for the moving source, RM5 signal, recorded on 27 October.
Timestamps are not available for the files; navigation coordinates for the moving source are also not available.

The 2D Initial Guess
Personal archival data with the site bathymetry in geographical coordinates were converted to UTM coordinates, as well as the VLA and stationary source GPS positions. To simplify the calculation of replicas, the obtained cartesian coordinates were further shifted in order to place the VLA at the origin of the bathymetry (see Figure 1); the expected range between the stationary source and the VLA range was calculated from these coordinates, providing the value r s = 5836 m. The source bearing relative to the VLA was found to be φ s = 93.3 • . Combined errors of GPS positions and wind-induced source drifting provide an error in the source range of approximately 400 m and a source bearing between 91 • and 95 • . From the description of experimental conditions, bearing estimates for the stationary and moving sources were expected to lie within this interval. The map shown in Figure 1 reveals a bathymetry in which depth increases mildly from right to left, with almost vertical isobaths; in particular, the isobath passing at the source slightly misses the VLA and gently deviates to the left in the region with Y < 0.
The bottom model adopted considered a single layer over a half space (see Figure 2); this model was adopted from the literature (see [5,17]). The corresponding properties are shown in Table 2.  Several sound speed profiles (SSPs) are provided with the acoustic data; however, different trials of replica generation demonstrated that the only SSP that can provide consistent estimates of source position for the 2D initial guess is the one shown in Figure 3. This is somehow surprising in terms of the stationarity of sound speed over two days of transmissions, and this is an issue that deserves future discussions. The SSP has unusual features: above 60 m depth, it is almost constant, while in the interval from 60 to 80 m, it exhibits a strong thermocline, with a negative gradient; below 80 m depth, it decreases almost linearly. In terms of ray propagation, the SSP divides the water column into two different waveguides: for a source placed above 60 m, rays will bounce within the surface and the bottom, no matter what the launching angle is, while for a source placed below 60 m, rays with small launching angles will bounce on the bottom repeatedly, and steeper rays will "escape" the thermocline and bounce within the two boundaries. However, ray field predictions were not used in source localization due to the fine-scale irregularities of the SSP and because ray models cannot account for layer properties; replicas were instead calculated with the KRAKEN normal mode model, which provides rather smooth predictions of the modes (see Figure 4)    During calculations, it was noticed through optimization than an additional tuning of water depth and signal frequency allowed us to improve the fit between replicas and observations, relative to the values indicated in the literature; thus, the 2D initial guess was generated using the tuned values shown in Table 3. In all cases, the fit between the replicas and the observations was calculated considering the narrowband Bartlett estimator, which can be written as [18] B(r, z, φ) = e * R(r s , z s , φ s )e , where e(r, z, φ) is the normalized replica calculated for a fixed frequency, * denotes the conjugate transpose, and R stands for the covariance matrix, which was calculated from received signals with MATLAB code provided with the acoustic data. When the dependence on the bearing is weak,B(r, z, φ) =B(r, z) and 2D localization can be achieved from the condition B(r s , z s ) = max (r,z)B (r, z) .
The set with 2D initial guesses for the acoustic data was obtained considering frequencies that maximized the Bartlett estimator, but ignoring any wind influence; the corresponding results are shown in Table 4. The Bartlett estimator is shown in two particular cases in Figure 5. The 2D initial guess for the stationary source at 331 Hz is rather consistent for all files, but underestimates the source-VLA range by almost 500 m; the corresponding source position is also different from the estimate for the stationary source at 171 Hz, which is itself very consistent and lies within the range of error. The guess for the moving source at 170 Hz clearly indicates a source moving away from the VLA. The exact nature of the differences in position for the stationary sources is unclear; it seems reasonable to consider that they reflect changes in wind direction and ocean currents along transmissions.

Three-Dimensional Localization and Tracking
Calculations of the 2D initial guess with the KRAKEN model require a reduced amount of computational time because bottom properties and sound speed are known; the next step for 3D localization and tracking consists in generating estimates of source bearing. To this end, one can rely on ESP for the generation of replicas; the source at depth z s is idealized as rotating around the VLA along a circle with a radius r s and N different bearings, and each bearing defines a different source-array configuration relative to the bathymetry and corresponds to a different set of replicas, to be calculated either through 3D or N × 2D modeling. Thus, bearing estimation can take place in the search space (r s , To reduce calculations as much as possible, the VLA is considered to be always at the origin; the bathymetry is also corrected, so the bottom depth at the VLA corresponds to the tuned value indicated in Table 3. As shown in Figure 6a, this approach produces a butterfly-like shape when the Bartlett estimator is displayed in polar coordinates. The butterfly's "wings" point in the directions where bearing estimation is more reliable, although the alignment of isobaths makes the estimation ambiguous. In the case of Figure 6a, the upper wing is aligned along φ up ≈ 89 • , while the lower wing aligns at a sort of antipode at φ down ≈ 250 • , which in fact has a higher Bartlett power. Thus, blind acceptance at this stage of the highest Bartlett would provide a source bearing, which deviates almost 180 • from the true bearing. Similar results can be seen in Figure 6b,c; for Figure 6b, the wings tend to be wider, and the Bartlett curves practically overlap each other, while in Figure 6c, the wings are somehow between the two cases of the stationary sources but with a Bartlett power that increases as the source-array distance increases. Results for the moving source are remarkable because of the constant wing alignment, a feature of Figure 6c that is consistent with the description of source motion presented in Section 2.
While the Bartlett estimators shown in Figure 6 do not resolve source bearing, they drastically reduce the search space from (r s , Since the error in the estimation of source depth can be expected to be much smaller than the error in source range, the strategy used to remove bearing ambiguity consisted in calculating additional replicas in the search space where δr is to be chosen to be sufficiently "small"; a faulty choice of δr can be detected when the localization result is found at one extreme of the search space. This strategy was applied to all files, and in each case it produced a "wing motion", in which the estimator wings extend or contract as range changes (see Figure 7). Bearing ambiguity was removed in all cases by noticing that the region φ up ± δφ exhibits the highest Bartlett power at a range r * ∈ [r s ± δr] and a bearing φ * ∈ [φ up ± δφ]. The wing motion provides a new estimate of source position, given by (r * , z s , φ * ); the adoption of this result as final depends on the desired degree of accuracy. However, further refinement can proceed inside the search space given by [r * ± δr] × [z s ± δz] × [φ * ± δφ] with a new set of small parameters (δr, δz, δφ), with the advantage that the calculation of replicas inside this search space is not computationally demanding; this provides the final 3D estimate of source position (r , z , φ ). For instance, for the stationary source at 331 Hz and file 1, it was found that (r , z ) = (5435, 65) m and φ = 93 • (see Figure 8); this 3D estimate corrects r s by almost 60 m, has no effect on source depth, and provides a value for the source bearing close to the expected one. A summary of localization results at 331 Hz for all files is presented in Table 5; one can notice that all range and bearing estimates are now within the expected boundaries of range and bearing errors. Corresponding results for the other two cases are presented in Tables 6 and 7; source positions over the (X, Y) plane are shown in Figure 9. Discretization values for range, depth, and bearing corresponded in all cases to 5 m, 1 m, and 1 • , respectively. Unlike the results for the stationary source at 331 Hz, there are a few cases with bearing estimates outside the expected boundaries of error; this issue indicates perhaps that the assumption of a fixed VLA needs additional analysis.  2  5570  76  93  3  5575  76  94  4  5575  76  94  5  5565  75  91  6  5555  75  90  7  5555  75  90  8  5560  75  91  9  5580  76  95  10 5580 76 94

Conclusions and Future Work
An approach to reduce the search space for the 3D localization and tracking of an acoustic source was proposed and tested on experimental data, with encouraging results. Critical to the success of this approach was the availability of a 2D initial guess, which was explored to determine a particular interval of source bearings; from this, it followed that the guess significantly reduced the search space of parameter combinations required to obtain the position of the source in a three-dimensional environment. Certainly, this guess cannot be expected to be available in all types of scenarios; this is an issue that will require further analysis considering additional sets of experimental data, particularly those in which out-of-plane propagation is known to be relevant. The broadband processing of acoustic receptions also needs to be considered, as well as replica calculations for a three-dimensional field of sound speed. Finally, the lack of navigational data does not allow us to determine the error of source localization and tracking; this issue is expected to be reviewed considering observations for which the geographical coordinates of the source are known. An association of the suggested approach with parallel processing and/or neural network architectures can be expected to speed up significantly the processing of acoustic data, which is required to locate an underwater source.
Author Contributions: Conceptualization, methodology and validation: O.C.R., L.Z. and X.C.; software: O.C.R. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.