Underwater target tracking by means of acoustic and electromagnetic ...
Underwater target tracking by means of acoustic and electromagnetic data fusion E. Dalberg, A. Lauberts, R.K. Lennartsson, M.J. Levonen and L. Persson Swedish Defence Research Agency SE-17290, Stockholm, Sweden (eva.dalberg, andris.lauberts, ron.lennartsson, mika.levonen, leifp)@foi.se Abstract - An interesting possibility for improved surveillance capabilities in littoral waters is the integration of multisensor systems by means of data fusion. Here, we describe how data fusion can be used for localisation and tracking of targets by means of passive underwater acoustic and electric field sensors. The data was fused using a Kalman filter. The filter was applied on bearing estimates from the acoustic data and estimates of the target position from the electrode array. The results show the positions of the targets are gained by fusion between the acoustic bearing and the electrode sensors, within the limits of their sensitivity.
Keywords: Tracking, localisation, bearing estimation, Kalman filter, data fusion.
1
Introduction
Passive underwater surveillance is dominated by acoustic methods. However, in littoral environments the hydrographic and bathymetric conditions can be such that acoustic systems have a limited range of operation. Also ambient noise (e.g. heavy shipping) may further degrade the system performance. If surveillance is performed at such conditions combined with low acoustic target signatures the use of additional sensor systems is an intersting option. In this paper we investigate the possibility of fusing data from an electrode sensor system with a hydrophone sensor system. The gains from fusing data could imply increased detection ranges or improved robustness when tracking low signature targets in areas with heavy shipping, or more reliable, automated target classification. The data fusion scheme, [1], involves several stages of cooperating and overlapping parts. Localisation, detection and classification may be done independently at each sensor systems and at different times. Data fusion can be done at several different levels. The first levels are at raw data and observation level, if the sensors are measuring the same physical phenomena. The fusion can also be done with decisions or features where the data is aggregated to estimates such as target bearings and positions. Here, we use bearing estimates from acoustic data together with position estimates from electromagnetic
data in a Kalman filter for fused target tracks. We describe the used methods and point out some lessons learned. We discuss briefly the electromagnetic and acoustic estimation methods followed by a description of the Kalmar filter tracker and finally some tracking results and conclusions.
2
Target localisation using underwater electric fields
Low frequency underwater electric fields from surface vessels or underwater vehicles are caused by electric currents flowing between the hull and the propulsion system of the vessel. The currents have their origin in the use of materials with different electrochemical potential that are in contact with the sea water which acts like an electrolyte. Typical frequencies of these signatures are in the Extreme Low Frequency (ELF) range, defined as 0.3-3000 Hz. In this frequency range the propagation conditions in shallow waters are such that the fields can be detected up to a few kilometres away from the target. These fields can be measured using pairs of electrodes. The electric potential Φ, is measured between the electrodes, and the electric field is derived as E = Φ/d, where d is the distance between the electrodes. This distance varies between a few metres in systems primarily used for signature measurements, up to a few hundred metres in surveillance systems optimised for long detection ranges. A target vessel can be modelled (in the far-field) as an electric dipole source, with the dipole moment as its source strength. The propagation was modelled using a horizontally stratified model of the underwater environment, where the most important parameters are the conductivity and thickness of each layer. Such a model can be more or less refined, but it is often sufficient to use a model where the sea water constitutes one layer in addition to one or two sediment layers and the bedrock. Using such an environmental model, together with an electric dipole model of the vessel, the electric field caused by a vessel can be calculated at any point in the area surrounding the target. If the environment in an area is known, measurements of the underwater electric field may be used to localise the source of the fields by minimising the residual between
the orthogonal decomposition of the covariance matrix Rxx
the calculated and the measured fields. f (u) = ||r(u)||2 ,
r(u) = m(u, w) − d
(1)
where f (u) is the square of the L2 -norm of the residual r(u). The vector m(u, w) is the data set predicted by the model, for a parameter vector u. The parameter vector w includes all the variables that are held constant during the minimisation, and d is the measured field. The unknown parameters are the horizontal position (x, y), and the heading of the source φ. This inversion scheme was used to track a source moving through an area where electric field sensors were deployed. The electric field modelling and inversions were done using the program packages nlayer [2], and inv nlayer [3], respectively. The model parameters for the specific environment had typical values estimated from previous experiments in the area [4]. Furthermore, acoustic bearing estimates were used to limit the allowed positions to a 10 degree circle sector with a maximum range of 2000 metres.
3
Acoustic bearing estimation
Several robust high-resolution methods for directionof-arrival (DOA) estimation of impinging wave fronts are available in the literature, e.g. [5],[6], [7] and [8]. In sonar array processing it is possible to estimate both the temporal (e.g. time-delays) and spatial (e.g. transmission loss) structure of the received wave field. However, the estimation of DOA of wavefronts, i.e. bearing estimation, is the most fundamental task in sonar array signal processing. To obtain bearing estimates several sensors are used to form sensor arrays in various configurations. This is done in order to gain advantages of directivity, resolution and beamwidth, compared to single sensor systems. The resolution for conventional beamforming is bounded by the wavelength of the received signals, the sensor configuration and the aperture of the array. In many underwater situations it is important to obtain a higher resolution than achievable with conventional beamforming. During the past decades several highresolution eigenstructure-based DOA estimators have been developed and are found in the literature, e.g. [5] and [9]. Common to all methods of DOA estimation is the signal model, based on the assumption that a source signal x(t) generates a wave field at a bearing θ relative to a linear array of M sensors. A general model for the received wave field consists of the source signals attenuated with αm , and then delayed with τm (θ) at sensor m as follows, xm (t) = αm x(t − τm (θ)) + φm (t).
(2)
The noise term φm (t) represents the additive ambient noise at each sensor. In this paper, the bearing estimation was done using the minimum variance (MV) algorithm described in [5], chapter 3. We use the standard formulation of high resolution methods formed by
Rxx =
N 1 X ˆ x ] [ x(t) − µ ˆ x ]H [ x(t) − µ N t=1
(3)
where N is the number of snapshots, H is the conjuˆ x is the estimated mean (over N gate transpose and µ snapshots). The covariance matrix is divided into a signal subspace Rss and a noise subspace Rnn , Rxx = Rss + Rnn
Rss =
P X
λr Vr VrH
(4)
(5)
r=1
Rnn =
M X
λr Vr VrH .
(6)
r=P +1
The signal subspace Rss (Eq. 5) is spanned by the P largest eigenvalues λr and eigenvectors Vr of Rxx , and it is separated from the noise subspace, Rnn (Eq. 6), formed by the smallest eigenvalues and vectors. The MV method is also often called the maximum likelihood method or the minimum energy method in the literature. The MV method implementation we use is formulated as a signal subspace method by using the singular value decomposition of the covariance matrix above, Eq. 4 as: PMV (Θ) =
1 CH (Θ)R−1 SS C(Θ)
(7)
where the steering vector C(Θ) is defined as C(Θ) = [1, ejΦ(1) , ..., ejΦ(M−1) ]′ .
(8)
Therefore, MV is applied on the reduced subspace spanned by the largest signal related eigenvalues for the bearing estimation.
4
Kalman filtering
Data fusion, in the context of this paper, is used for target tracking by combining target position information from an electrode array and target bearings from an hydrophone array. The Kalman filter provides a general solution to the recursive minimised mean square linear estimation problem. The mean square error will be minimised as long as the target dynamics and the measurement noise are accurately modelled. In addition, the Kalman filter provides a convenient measure of the estimation accuracy through the covariance matrix, and the gain sequence automatically adapts to the variability of the data.
4.1
Kalman filter definition
Derivations of the Kalman filter are available in the literature, e.g. [10], [11] and [12]. Assume that the
target dynamic process can be modelled in the discretetime Markov form x(k + 1) = Φx(k) + q(k) + f(k + 1|k)
(9)
where x(k) is the n-dimensional target state vector with time index k, Φ is the assumed state transition matrix, q(k) is the zero-mean, white, Gaussian process noise with assumed known covariance Q, and f(k+1|k) is the assumed known deterministic input. For this application, the discrete-time Markov process, described by the difference equation in Eq. 9, implies that the statistical representation of the process in the future (k + 1) is completely determined by the present state (k). Assume further that the steps k are advanced by equal time steps T . Measurements are linear combinations of the system state variables, corrupted by additive, uncorrelated noise, defined as y(k) = Hx(k) + v(k)
Tracker Implementation
The tracker implementation here mainly adopts the algorithm in [13]. In the following, a 2-D planar system with state vector x = [x, y, vx , vy ] is assumed. Track initiation is done by selecting ”rule-of-thumb” parameter values for measurement noise, process noise, state information, state estimate and covariance. Since we expect slow manoeuvres, the state transition matrix Φ, the process noise covariance matrix Q, the (unrotated) noise covariance matrices R0 (one for each sensor type), and the measurement matrix H have constant, simple forms
1
ˆ (k|k) = x ˆ (k|k − 1) + · · · x
0 Φ= 0 0
0
1
0
0
0
T 0 1
0
0
(17)
0
1
0
0
1
0
0
0 0 1
σx2
0
0
σy2
1
0
0
0
0
1
0
0
(18)
(19)
(11)
H= (12)
P(k|k) = · · · [I − K(k)H]P(k|k − 1)[I − K(k)H]′ + · · · K(k)RK′ (k)
0
1
R0 =
K(k) = · · · −1
T
1
0 Q = q 0 0
P(k|k − 1)H′ [HP(k|k − 1)H′ + R]
0
(10)
where y is the M -dimensional measurement vector, H is the M × n measurement matrix, and v is zero-mean, white, Gaussian noise with covariance R. Note that in general Q and H may also vary with time, though kept constant in our application. Given the target dynamics and measurement models from Eq. 9 and Eq. 10, the Kalman filter equations become
K(k) [y(k) − Hˆ x(k|k − 1)]
4.2
(20)
where T is the time step, q is a scale factor, and the components σx2 and σy2 correspond to measurement x and y errors. The state x is initialised to x0 = [x0 , y0 , 0, 0], taking the coordinates of first accepted detection. The state covariance P is initialised to a unit matrix times a scale factor. The matrix V co-rotates R0 through a bearing θ into R
(13) R = VR0 V′
ˆ (k + 1|k) = Φˆ x x(k|k) + f(k + 1|k)
′
P(k + 1|k) = ΦP(k|k)Φ + Q
(14)
(15)
where [..]′ denotes the transpose, and Eq. 12 gives the stabilised form of the covariance update. Note that in Eq. 14 the term f(k + 1|k) is zero (no ownship motion). The state covariance matrix is defined in terms of the zero-mean Gaussian estimation error vector ˆ ] [x(k) − x ˆ ]′ }. P(k) = E{[x(k) − x
(16)
V=
cos θ
− sin θ
sin θ
cos θ
(21)
.
(22)
Typically, a new target is first observed by the nearest or the most sensitive sensor. An ellipse, derived from an eigenvalue decomposition of the R matrix, indicates the observation uncertainty region (highly elongated in the bearing direction having bearing data only). The first position estimate of the target is in the centre of this ellipse. Subsequent observations from the same sensor add more information. Eventually, other
sensors contribute with independent observations improving the tracking. An important point is the association of an position estimate to the target. Therefore, reliable observations are gated within a window of maximum size, maxgate satisfying the equation below. (y − x ˆ)′ P−1 ˆ) < maxgate, p (y − x
(23)
where Pp is a sub matrix of the state covariance containing only the position covariance. The following points outline the main steps in the Kalman filter process: 1. Initialize state vector x (position and velocity for single target) and covariance P. Take the start position as the first electromagnetic (EM) measurement and set the velocity to zero. Assume the measurement errors for EM and acoustic (Ac) sensors with fixed statistics expressed by their covariance matrices. 2. Convert measurements y to rectangular coordinates. In particular, the Ac sensor bearings are complemented by a rough ”mean distance”. Its large uncertainty is represented by a co-rotating very elongated error ellipse. 3. Look for new observations at regular times with step T. Note that the EM measurements need not coincide with the Ac measurements. 4. If the measured position falls within a fixed gate (Eq. 23) then update the Kalman filter equations and estimate a new state vector. 5. If the measurement is missing or falls outside the gate then do a new state estimate by only using the past settings of the Kalman filter. 6. Repeat steps 2-5. Note that EM and acoustic sensor measurements yield multiple model answers with different weights. Here in this study, the most probable measurement is chosen for each sensor at a given time.
5
Tracking targets using data from the sea trial
Here, we study both electric field and acoustic data from a sea trial conducted in the Baltic Sea at the southern archipelago of Stockholm. The sea trial area is typical for this part of the Baltic Sea with several small islands and a high variability of the bathymetric conditions. Further details on the experiment can be found in [14]. The electric field measurement system consisted of a linear array with eight carbon fibre electrodes. The distance between the sensors varies from 50 to 200 m. The array was placed on the sea floor in a strait where the water depth varied between 10-45 m. A synthetic dipole source was towed at a depth of 2 m along different tracks in the area. The position of the
synthetic source was continuously logged with a GPSsystem. Continuous sinusoids of a few Hz were transmitted at different dipole moment strengths. Here, we present the evaluation on estimated results from four tracks. During track 1 a low dipole strength of around 4 Am and 4 Hz was used. The signal-to-noise ratio (SNR) was high enough in four of the sensors, whereas the other four sensors were not used in the analysis of that track. On the other hand, the depth profile of that track was less varied than the other three tracks. During tracks 2, 3a, and 3b, the source strength was 28 Am, which enabled the use of all eight sensors. During the tracks 3a and 3b the target vessel was manoeuvred over relatively constant depths, whereas the bathymetry along track 2 varied more. The acoustic data was recorded using a 31-element uniform line array, with a sensor spacing of 1.5 m, placed on the seafloor at a water depth of around 42 m. Half of that array was placed in a vertical position during tracks 3a and 3b. In this case, only data from the 15 elements remaining in the horizontal part of the array has been used. The hydrophone array was placed some 100 m west of the electrode array. An acoustic source was towed after a vessel in parallel with the electromagnetic source during tracks 2, 3a and 3b, while the only acoustic source during track 1 was the vessel. The acoustic source transmitted a continuous wave at 200 Hz. The bearing estimation was estimated in a frequency band between 150 and 300 Hz, and the bearing track used in the data fusion was limited to bearing changes of less than 10 degrees, between consecutive bearing estimates.
5.1
Track 1
Acoustic (Ac) and electromagnetic (EM) data from track 1 have been analysed resulting in coordinates for the target. The acoustic array delivered bearings and the EM array gave positions of a signal projector towed after a vessel. A GPS on board the vessel provided reference positions showing the true track. The acoustic signals were segmented into 1 s windows and the EM signals were segmented into 5 s windows suggesting a 0.5 s time step in the Kalman filter. For each data point we have assumed that the acoustic array gives the target bearing within ±5 degrees about the true direction. Consequently the EM data are mainly used to get the target distance. The electromagnetic measurement covariance has been scaled with the local mean square residual between model and data, influencing the Kalman tracker. Occasional large residuals thus make a negligible contribution to the track estimate. Figure 1 shows a track estimate using data from track 1. The track estimate starts in the upper left corner of the figure, and is estimated from the first EM measurement. The initially large error ellipses indicate an uncertain starting point, followed by some fairly quick steps past the hydrophone array centre. The cloud of EM points surrounding position (100,-800) lie at the limit of EM sensitivity so the estimates rely on acoustic bearings only. However, not surprisingly this
1000
0
EM center
EM center GPS
500 −500 y [m]
y [m]
0 Ac center GPS
−500
Estimate
−1000
Ac center
GPS
path −1000 Estimate −1500 0
500 x [m]
1000
−1500 −1500
1500
Figure 1: Track 1. Estimated target track (red dots each time step) using bearings from an acoustic array (Ac, big blue cross) and position data from an electrode array (EM, big black ring). Black dots denote calculated EM positions, green crosses denote GPS positions. State error ellipses are shown in blue every 100 steps. Relatively high weight is put on the acoustic data compared with the EM data.
−500
0 x [m]
500
1000
1500
Figure 2: Track 2. Estimated target track (red dots each time step) using bearings from an acoustic array and position data from an electrode array. See Figure 1 for further details.
500
GPS
EM
GPS
0 y [m]
later part of the estimated track deviates quite a lot from the track given by the GPS readings. The final track estimate is running freely, out of reach for both sensor systems.
−1000
−500 Estimate
5.2
Track 2
Figure 2 shows a track estimate using data from track 2. Note how the state error ellipses steadily grow after passing the swarm of EM positions just right of the picture centre. A bit surprising is that using EM data only the resulting track turns out to be about the same, indicating the weak influence of the accompanying acoustic array. The large discrepancy between EM and GPS positions in this part of the track evidently calls for new interpretation of the EM data or extra supporting sensors. Using the GPS positions a second hydrophone array positioned at (1400,-1100) was simulated, see Figure 3. The simulated acoustic array (Ac2) outputs bearings at the same rate and error characteristics as the real one (Ac1). This time, thanks to the extra information from the second acoustic array, the track estimate follows the true track (GPS readings) much better, completely avoiding the central cloud of false EM positions. It should be noted that the second array must be carefully positioned in order to optimally support the bearings from the first array.
5.3
Tracks 3a and 3b
Figure 4 shows a track estimate using data from track 3a, using positions from the EM array and bearings from the acoustic array. Except from an initial discrepancy between the track estimate and GPS positions, the main remaining part of the track estimate
Ac1 −1000 Ac2 −1500 −1000
−500
0 x [m]
500
1000
1500
Figure 3: Track 2. Estimated target track (red dots) using bearings from two acoustic arrays, Ac1 and the simulated array, Ac2, together with position data from an electrode array, EM. See Figure 1 for further details. follows the GPS positions rather well. The differently sized error ellipses reflect the changing state estimate errors, getting smaller towards the EM centre. The actual position error, however, has a local maximum at this point. This is because the used EM model does not take into account near-field effects. The acoustic array, on the other hand, has the best performance when the resulting track is broadside to the array. In this case, the target passes the end of the array. This means that the geometry of the track and the arrays are such that the acoustic bearings are of little help here. The returning part of the track, running from right to left in Figure 5, uses EM and acoustic data from track 3b. Except from a uncertain start, the remaining track estimate lies somewhat closer to the GPS readings than in case 3a. Again there is a maximal deviation between the electric field positions and acoustic
500 path
Estimate
0
y [m]
−500 GPS
−1000
−1500
−1000
−500
0
500
1000
1500
x [m]
Figure 4: Track 3a. EM positions fused with acoustic bearings. The track goes from left to right. See Figure 1 for further details. 500 Estimate
0
y [m]
−500 path
GPS
−1000
−1500
−1000
−500
0
500
1000
1500
x [m]
Figure 5: Track 3b. EM positions fused with acoustic bearings. The track goes from right to left. See Figure 1 for further details. bearings, which could be attributed to the particular layout of these arrays, i.e. near endfire position for the acoustic array and near-field effects for the electrode array.
6
Conclusions
Tracking was performed on real data from a sea trial in the Baltic Sea, where two passive underwater sensor systems were used. One of the systems was an acoustic uniform line array placed on the seabed, and the other one was a long base line electrode array for underwater electric field measurements. The sea trial took place in shallow water, where acoustic and electric field sources were towed from a ship, transmitting signals simulating an underwater target. Bearing information from the acoustic array has been fused with target positions derived from the electric field using an inversion algorithm. The bearing- and positioning information was
fused using a Kalman filter to estimate the track of the target. This technique can be useful for applications where accurate target tracking is needed within a limited area. Acoustic beamforming is a well established technique, and there are numerous algorithms that have been tuned for different applications. Tracking using electric field sensors is a fairly new discipline. This work shows that it is possible to use methods for localisation of underwater targets based on electric field together with acoustic bearing information to form a fused target track. Compared to the acoustic sensors, the electro magnetic sensors have relatively short detection ranges. In addition to that, target ranges are underestimated in this analysis, even when the signal-to-noise ratio is high. This is probably caused by a mismatch between the true underwater environment and the used environment model. Filtering the EM data, constrained by acoustic bearings, helps to some extent. A real improvement can be achieved by adding carefully positioned extra sensors. Achieving the best possible performance gain from multisensor data fusion requires complete control over the entire chain from array design to the fusion algorithm. Our focus was kept on passive acoustic and electric field underwater surveillance, which leaves many possible approaches of data fusion unexamined. Another possibility would be to use combined matched field processing. It would also be of interest to evaluate the relative benefits of using one acoustic and one electric field array, compared with two acoustic arrays.
References [1] D.L. Hall and J. Llinas, An introduction to multisensor data fusion, Proc. of the IEEE, vol 85, no 1, pp. 6-23, 1997. [2] L. Abrahamsson and B-L. Andersson, User’s Guide to NLAYER 2.0, FOA-R–01310-409–SE, 1999. [3] P. Krylstedt and J. Mattsson, INV NLAYSCA, FOA-R–00-01548-409–SE, 2000. [4] P. Krylstedt and J. Mattsson, Environment assessment for underwater electric sensors, Proc. UDT Europe, Malm¨o, 2003. [5] S. Haykin (ed.), Array Signal Prentice-Hall, New Jersey, 1985.
Processing,
[6] R. O. Schmidt, Multiple emitter location and signal parameter estimation, in Proc. RADC Spectral Estimat. Workshop, pp. 243-258, 1979. [7] T. J. Shah and T. Kailath, Adaptive beamforming for coherent signals and interference, IEEE Trans. Acoust., Speech, Signal Processing, ASSP-33, pp. 527-536, 1985.
[8] E. Dudgeon, Fundamentals of digital array processing, Proc. of the IEEE, vol 65, no 6, pp. 898904, 1977. [9] P.S. Naidu, Sensor Array Signal Processing, CRC Press LLC, 2001. [10] S. Blackman and R. Popoli, Design and Analysis of Modern Tracking Systems, Boston, MA, Artech, pp. 157-160, 1999. [11] A. S. Gelb, Applied Optimal Estimation, Cambridge, MA, MIT Press, 1974. [12] Y. Bar-Shalom and T. E. Fortmann, Tracking and Data Association, Orlando, FL, Academic Press, 1988. [13] K. Murphy, Homepage for Matlab Kalman filter toolbox, http://www.cs.ubc.ca/ murphyk/Software/Kalman/kalman.html. [14] L. Crona, E. Dalberg, R. Lennartsson, B. Lundqvist, Per M´ oren, L. Persson, P. S¨ oderberg, An experiment on passive multi-sensor underwater surveillance, FOI-R–1336–SE, 2004.
Keywords: Tracking, localisation, bearing estimation, Kalman filter, data fusion.
1
Introduction
Passive underwater surveillance is dominated by acoustic methods. However, in littoral environments the hydrographic and bathymetric conditions can be such that acoustic systems have a limited range of operation. Also ambient noise (e.g. heavy shipping) may further degrade the system performance. If surveillance is performed at such conditions combined with low acoustic target signatures the use of additional sensor systems is an intersting option. In this paper we investigate the possibility of fusing data from an electrode sensor system with a hydrophone sensor system. The gains from fusing data could imply increased detection ranges or improved robustness when tracking low signature targets in areas with heavy shipping, or more reliable, automated target classification. The data fusion scheme, [1], involves several stages of cooperating and overlapping parts. Localisation, detection and classification may be done independently at each sensor systems and at different times. Data fusion can be done at several different levels. The first levels are at raw data and observation level, if the sensors are measuring the same physical phenomena. The fusion can also be done with decisions or features where the data is aggregated to estimates such as target bearings and positions. Here, we use bearing estimates from acoustic data together with position estimates from electromagnetic
data in a Kalman filter for fused target tracks. We describe the used methods and point out some lessons learned. We discuss briefly the electromagnetic and acoustic estimation methods followed by a description of the Kalmar filter tracker and finally some tracking results and conclusions.
2
Target localisation using underwater electric fields
Low frequency underwater electric fields from surface vessels or underwater vehicles are caused by electric currents flowing between the hull and the propulsion system of the vessel. The currents have their origin in the use of materials with different electrochemical potential that are in contact with the sea water which acts like an electrolyte. Typical frequencies of these signatures are in the Extreme Low Frequency (ELF) range, defined as 0.3-3000 Hz. In this frequency range the propagation conditions in shallow waters are such that the fields can be detected up to a few kilometres away from the target. These fields can be measured using pairs of electrodes. The electric potential Φ, is measured between the electrodes, and the electric field is derived as E = Φ/d, where d is the distance between the electrodes. This distance varies between a few metres in systems primarily used for signature measurements, up to a few hundred metres in surveillance systems optimised for long detection ranges. A target vessel can be modelled (in the far-field) as an electric dipole source, with the dipole moment as its source strength. The propagation was modelled using a horizontally stratified model of the underwater environment, where the most important parameters are the conductivity and thickness of each layer. Such a model can be more or less refined, but it is often sufficient to use a model where the sea water constitutes one layer in addition to one or two sediment layers and the bedrock. Using such an environmental model, together with an electric dipole model of the vessel, the electric field caused by a vessel can be calculated at any point in the area surrounding the target. If the environment in an area is known, measurements of the underwater electric field may be used to localise the source of the fields by minimising the residual between
the orthogonal decomposition of the covariance matrix Rxx
the calculated and the measured fields. f (u) = ||r(u)||2 ,
r(u) = m(u, w) − d
(1)
where f (u) is the square of the L2 -norm of the residual r(u). The vector m(u, w) is the data set predicted by the model, for a parameter vector u. The parameter vector w includes all the variables that are held constant during the minimisation, and d is the measured field. The unknown parameters are the horizontal position (x, y), and the heading of the source φ. This inversion scheme was used to track a source moving through an area where electric field sensors were deployed. The electric field modelling and inversions were done using the program packages nlayer [2], and inv nlayer [3], respectively. The model parameters for the specific environment had typical values estimated from previous experiments in the area [4]. Furthermore, acoustic bearing estimates were used to limit the allowed positions to a 10 degree circle sector with a maximum range of 2000 metres.
3
Acoustic bearing estimation
Several robust high-resolution methods for directionof-arrival (DOA) estimation of impinging wave fronts are available in the literature, e.g. [5],[6], [7] and [8]. In sonar array processing it is possible to estimate both the temporal (e.g. time-delays) and spatial (e.g. transmission loss) structure of the received wave field. However, the estimation of DOA of wavefronts, i.e. bearing estimation, is the most fundamental task in sonar array signal processing. To obtain bearing estimates several sensors are used to form sensor arrays in various configurations. This is done in order to gain advantages of directivity, resolution and beamwidth, compared to single sensor systems. The resolution for conventional beamforming is bounded by the wavelength of the received signals, the sensor configuration and the aperture of the array. In many underwater situations it is important to obtain a higher resolution than achievable with conventional beamforming. During the past decades several highresolution eigenstructure-based DOA estimators have been developed and are found in the literature, e.g. [5] and [9]. Common to all methods of DOA estimation is the signal model, based on the assumption that a source signal x(t) generates a wave field at a bearing θ relative to a linear array of M sensors. A general model for the received wave field consists of the source signals attenuated with αm , and then delayed with τm (θ) at sensor m as follows, xm (t) = αm x(t − τm (θ)) + φm (t).
(2)
The noise term φm (t) represents the additive ambient noise at each sensor. In this paper, the bearing estimation was done using the minimum variance (MV) algorithm described in [5], chapter 3. We use the standard formulation of high resolution methods formed by
Rxx =
N 1 X ˆ x ] [ x(t) − µ ˆ x ]H [ x(t) − µ N t=1
(3)
where N is the number of snapshots, H is the conjuˆ x is the estimated mean (over N gate transpose and µ snapshots). The covariance matrix is divided into a signal subspace Rss and a noise subspace Rnn , Rxx = Rss + Rnn
Rss =
P X
λr Vr VrH
(4)
(5)
r=1
Rnn =
M X
λr Vr VrH .
(6)
r=P +1
The signal subspace Rss (Eq. 5) is spanned by the P largest eigenvalues λr and eigenvectors Vr of Rxx , and it is separated from the noise subspace, Rnn (Eq. 6), formed by the smallest eigenvalues and vectors. The MV method is also often called the maximum likelihood method or the minimum energy method in the literature. The MV method implementation we use is formulated as a signal subspace method by using the singular value decomposition of the covariance matrix above, Eq. 4 as: PMV (Θ) =
1 CH (Θ)R−1 SS C(Θ)
(7)
where the steering vector C(Θ) is defined as C(Θ) = [1, ejΦ(1) , ..., ejΦ(M−1) ]′ .
(8)
Therefore, MV is applied on the reduced subspace spanned by the largest signal related eigenvalues for the bearing estimation.
4
Kalman filtering
Data fusion, in the context of this paper, is used for target tracking by combining target position information from an electrode array and target bearings from an hydrophone array. The Kalman filter provides a general solution to the recursive minimised mean square linear estimation problem. The mean square error will be minimised as long as the target dynamics and the measurement noise are accurately modelled. In addition, the Kalman filter provides a convenient measure of the estimation accuracy through the covariance matrix, and the gain sequence automatically adapts to the variability of the data.
4.1
Kalman filter definition
Derivations of the Kalman filter are available in the literature, e.g. [10], [11] and [12]. Assume that the
target dynamic process can be modelled in the discretetime Markov form x(k + 1) = Φx(k) + q(k) + f(k + 1|k)
(9)
where x(k) is the n-dimensional target state vector with time index k, Φ is the assumed state transition matrix, q(k) is the zero-mean, white, Gaussian process noise with assumed known covariance Q, and f(k+1|k) is the assumed known deterministic input. For this application, the discrete-time Markov process, described by the difference equation in Eq. 9, implies that the statistical representation of the process in the future (k + 1) is completely determined by the present state (k). Assume further that the steps k are advanced by equal time steps T . Measurements are linear combinations of the system state variables, corrupted by additive, uncorrelated noise, defined as y(k) = Hx(k) + v(k)
Tracker Implementation
The tracker implementation here mainly adopts the algorithm in [13]. In the following, a 2-D planar system with state vector x = [x, y, vx , vy ] is assumed. Track initiation is done by selecting ”rule-of-thumb” parameter values for measurement noise, process noise, state information, state estimate and covariance. Since we expect slow manoeuvres, the state transition matrix Φ, the process noise covariance matrix Q, the (unrotated) noise covariance matrices R0 (one for each sensor type), and the measurement matrix H have constant, simple forms
1
ˆ (k|k) = x ˆ (k|k − 1) + · · · x
0 Φ= 0 0
0
1
0
0
0
T 0 1
0
0
(17)
0
1
0
0
1
0
0
0 0 1
σx2
0
0
σy2
1
0
0
0
0
1
0
0
(18)
(19)
(11)
H= (12)
P(k|k) = · · · [I − K(k)H]P(k|k − 1)[I − K(k)H]′ + · · · K(k)RK′ (k)
0
1
R0 =
K(k) = · · · −1
T
1
0 Q = q 0 0
P(k|k − 1)H′ [HP(k|k − 1)H′ + R]
0
(10)
where y is the M -dimensional measurement vector, H is the M × n measurement matrix, and v is zero-mean, white, Gaussian noise with covariance R. Note that in general Q and H may also vary with time, though kept constant in our application. Given the target dynamics and measurement models from Eq. 9 and Eq. 10, the Kalman filter equations become
K(k) [y(k) − Hˆ x(k|k − 1)]
4.2
(20)
where T is the time step, q is a scale factor, and the components σx2 and σy2 correspond to measurement x and y errors. The state x is initialised to x0 = [x0 , y0 , 0, 0], taking the coordinates of first accepted detection. The state covariance P is initialised to a unit matrix times a scale factor. The matrix V co-rotates R0 through a bearing θ into R
(13) R = VR0 V′
ˆ (k + 1|k) = Φˆ x x(k|k) + f(k + 1|k)
′
P(k + 1|k) = ΦP(k|k)Φ + Q
(14)
(15)
where [..]′ denotes the transpose, and Eq. 12 gives the stabilised form of the covariance update. Note that in Eq. 14 the term f(k + 1|k) is zero (no ownship motion). The state covariance matrix is defined in terms of the zero-mean Gaussian estimation error vector ˆ ] [x(k) − x ˆ ]′ }. P(k) = E{[x(k) − x
(16)
V=
cos θ
− sin θ
sin θ
cos θ
(21)
.
(22)
Typically, a new target is first observed by the nearest or the most sensitive sensor. An ellipse, derived from an eigenvalue decomposition of the R matrix, indicates the observation uncertainty region (highly elongated in the bearing direction having bearing data only). The first position estimate of the target is in the centre of this ellipse. Subsequent observations from the same sensor add more information. Eventually, other
sensors contribute with independent observations improving the tracking. An important point is the association of an position estimate to the target. Therefore, reliable observations are gated within a window of maximum size, maxgate satisfying the equation below. (y − x ˆ)′ P−1 ˆ) < maxgate, p (y − x
(23)
where Pp is a sub matrix of the state covariance containing only the position covariance. The following points outline the main steps in the Kalman filter process: 1. Initialize state vector x (position and velocity for single target) and covariance P. Take the start position as the first electromagnetic (EM) measurement and set the velocity to zero. Assume the measurement errors for EM and acoustic (Ac) sensors with fixed statistics expressed by their covariance matrices. 2. Convert measurements y to rectangular coordinates. In particular, the Ac sensor bearings are complemented by a rough ”mean distance”. Its large uncertainty is represented by a co-rotating very elongated error ellipse. 3. Look for new observations at regular times with step T. Note that the EM measurements need not coincide with the Ac measurements. 4. If the measured position falls within a fixed gate (Eq. 23) then update the Kalman filter equations and estimate a new state vector. 5. If the measurement is missing or falls outside the gate then do a new state estimate by only using the past settings of the Kalman filter. 6. Repeat steps 2-5. Note that EM and acoustic sensor measurements yield multiple model answers with different weights. Here in this study, the most probable measurement is chosen for each sensor at a given time.
5
Tracking targets using data from the sea trial
Here, we study both electric field and acoustic data from a sea trial conducted in the Baltic Sea at the southern archipelago of Stockholm. The sea trial area is typical for this part of the Baltic Sea with several small islands and a high variability of the bathymetric conditions. Further details on the experiment can be found in [14]. The electric field measurement system consisted of a linear array with eight carbon fibre electrodes. The distance between the sensors varies from 50 to 200 m. The array was placed on the sea floor in a strait where the water depth varied between 10-45 m. A synthetic dipole source was towed at a depth of 2 m along different tracks in the area. The position of the
synthetic source was continuously logged with a GPSsystem. Continuous sinusoids of a few Hz were transmitted at different dipole moment strengths. Here, we present the evaluation on estimated results from four tracks. During track 1 a low dipole strength of around 4 Am and 4 Hz was used. The signal-to-noise ratio (SNR) was high enough in four of the sensors, whereas the other four sensors were not used in the analysis of that track. On the other hand, the depth profile of that track was less varied than the other three tracks. During tracks 2, 3a, and 3b, the source strength was 28 Am, which enabled the use of all eight sensors. During the tracks 3a and 3b the target vessel was manoeuvred over relatively constant depths, whereas the bathymetry along track 2 varied more. The acoustic data was recorded using a 31-element uniform line array, with a sensor spacing of 1.5 m, placed on the seafloor at a water depth of around 42 m. Half of that array was placed in a vertical position during tracks 3a and 3b. In this case, only data from the 15 elements remaining in the horizontal part of the array has been used. The hydrophone array was placed some 100 m west of the electrode array. An acoustic source was towed after a vessel in parallel with the electromagnetic source during tracks 2, 3a and 3b, while the only acoustic source during track 1 was the vessel. The acoustic source transmitted a continuous wave at 200 Hz. The bearing estimation was estimated in a frequency band between 150 and 300 Hz, and the bearing track used in the data fusion was limited to bearing changes of less than 10 degrees, between consecutive bearing estimates.
5.1
Track 1
Acoustic (Ac) and electromagnetic (EM) data from track 1 have been analysed resulting in coordinates for the target. The acoustic array delivered bearings and the EM array gave positions of a signal projector towed after a vessel. A GPS on board the vessel provided reference positions showing the true track. The acoustic signals were segmented into 1 s windows and the EM signals were segmented into 5 s windows suggesting a 0.5 s time step in the Kalman filter. For each data point we have assumed that the acoustic array gives the target bearing within ±5 degrees about the true direction. Consequently the EM data are mainly used to get the target distance. The electromagnetic measurement covariance has been scaled with the local mean square residual between model and data, influencing the Kalman tracker. Occasional large residuals thus make a negligible contribution to the track estimate. Figure 1 shows a track estimate using data from track 1. The track estimate starts in the upper left corner of the figure, and is estimated from the first EM measurement. The initially large error ellipses indicate an uncertain starting point, followed by some fairly quick steps past the hydrophone array centre. The cloud of EM points surrounding position (100,-800) lie at the limit of EM sensitivity so the estimates rely on acoustic bearings only. However, not surprisingly this
1000
0
EM center
EM center GPS
500 −500 y [m]
y [m]
0 Ac center GPS
−500
Estimate
−1000
Ac center
GPS
path −1000 Estimate −1500 0
500 x [m]
1000
−1500 −1500
1500
Figure 1: Track 1. Estimated target track (red dots each time step) using bearings from an acoustic array (Ac, big blue cross) and position data from an electrode array (EM, big black ring). Black dots denote calculated EM positions, green crosses denote GPS positions. State error ellipses are shown in blue every 100 steps. Relatively high weight is put on the acoustic data compared with the EM data.
−500
0 x [m]
500
1000
1500
Figure 2: Track 2. Estimated target track (red dots each time step) using bearings from an acoustic array and position data from an electrode array. See Figure 1 for further details.
500
GPS
EM
GPS
0 y [m]
later part of the estimated track deviates quite a lot from the track given by the GPS readings. The final track estimate is running freely, out of reach for both sensor systems.
−1000
−500 Estimate
5.2
Track 2
Figure 2 shows a track estimate using data from track 2. Note how the state error ellipses steadily grow after passing the swarm of EM positions just right of the picture centre. A bit surprising is that using EM data only the resulting track turns out to be about the same, indicating the weak influence of the accompanying acoustic array. The large discrepancy between EM and GPS positions in this part of the track evidently calls for new interpretation of the EM data or extra supporting sensors. Using the GPS positions a second hydrophone array positioned at (1400,-1100) was simulated, see Figure 3. The simulated acoustic array (Ac2) outputs bearings at the same rate and error characteristics as the real one (Ac1). This time, thanks to the extra information from the second acoustic array, the track estimate follows the true track (GPS readings) much better, completely avoiding the central cloud of false EM positions. It should be noted that the second array must be carefully positioned in order to optimally support the bearings from the first array.
5.3
Tracks 3a and 3b
Figure 4 shows a track estimate using data from track 3a, using positions from the EM array and bearings from the acoustic array. Except from an initial discrepancy between the track estimate and GPS positions, the main remaining part of the track estimate
Ac1 −1000 Ac2 −1500 −1000
−500
0 x [m]
500
1000
1500
Figure 3: Track 2. Estimated target track (red dots) using bearings from two acoustic arrays, Ac1 and the simulated array, Ac2, together with position data from an electrode array, EM. See Figure 1 for further details. follows the GPS positions rather well. The differently sized error ellipses reflect the changing state estimate errors, getting smaller towards the EM centre. The actual position error, however, has a local maximum at this point. This is because the used EM model does not take into account near-field effects. The acoustic array, on the other hand, has the best performance when the resulting track is broadside to the array. In this case, the target passes the end of the array. This means that the geometry of the track and the arrays are such that the acoustic bearings are of little help here. The returning part of the track, running from right to left in Figure 5, uses EM and acoustic data from track 3b. Except from a uncertain start, the remaining track estimate lies somewhat closer to the GPS readings than in case 3a. Again there is a maximal deviation between the electric field positions and acoustic
500 path
Estimate
0
y [m]
−500 GPS
−1000
−1500
−1000
−500
0
500
1000
1500
x [m]
Figure 4: Track 3a. EM positions fused with acoustic bearings. The track goes from left to right. See Figure 1 for further details. 500 Estimate
0
y [m]
−500 path
GPS
−1000
−1500
−1000
−500
0
500
1000
1500
x [m]
Figure 5: Track 3b. EM positions fused with acoustic bearings. The track goes from right to left. See Figure 1 for further details. bearings, which could be attributed to the particular layout of these arrays, i.e. near endfire position for the acoustic array and near-field effects for the electrode array.
6
Conclusions
Tracking was performed on real data from a sea trial in the Baltic Sea, where two passive underwater sensor systems were used. One of the systems was an acoustic uniform line array placed on the seabed, and the other one was a long base line electrode array for underwater electric field measurements. The sea trial took place in shallow water, where acoustic and electric field sources were towed from a ship, transmitting signals simulating an underwater target. Bearing information from the acoustic array has been fused with target positions derived from the electric field using an inversion algorithm. The bearing- and positioning information was
fused using a Kalman filter to estimate the track of the target. This technique can be useful for applications where accurate target tracking is needed within a limited area. Acoustic beamforming is a well established technique, and there are numerous algorithms that have been tuned for different applications. Tracking using electric field sensors is a fairly new discipline. This work shows that it is possible to use methods for localisation of underwater targets based on electric field together with acoustic bearing information to form a fused target track. Compared to the acoustic sensors, the electro magnetic sensors have relatively short detection ranges. In addition to that, target ranges are underestimated in this analysis, even when the signal-to-noise ratio is high. This is probably caused by a mismatch between the true underwater environment and the used environment model. Filtering the EM data, constrained by acoustic bearings, helps to some extent. A real improvement can be achieved by adding carefully positioned extra sensors. Achieving the best possible performance gain from multisensor data fusion requires complete control over the entire chain from array design to the fusion algorithm. Our focus was kept on passive acoustic and electric field underwater surveillance, which leaves many possible approaches of data fusion unexamined. Another possibility would be to use combined matched field processing. It would also be of interest to evaluate the relative benefits of using one acoustic and one electric field array, compared with two acoustic arrays.
References [1] D.L. Hall and J. Llinas, An introduction to multisensor data fusion, Proc. of the IEEE, vol 85, no 1, pp. 6-23, 1997. [2] L. Abrahamsson and B-L. Andersson, User’s Guide to NLAYER 2.0, FOA-R–01310-409–SE, 1999. [3] P. Krylstedt and J. Mattsson, INV NLAYSCA, FOA-R–00-01548-409–SE, 2000. [4] P. Krylstedt and J. Mattsson, Environment assessment for underwater electric sensors, Proc. UDT Europe, Malm¨o, 2003. [5] S. Haykin (ed.), Array Signal Prentice-Hall, New Jersey, 1985.
Processing,
[6] R. O. Schmidt, Multiple emitter location and signal parameter estimation, in Proc. RADC Spectral Estimat. Workshop, pp. 243-258, 1979. [7] T. J. Shah and T. Kailath, Adaptive beamforming for coherent signals and interference, IEEE Trans. Acoust., Speech, Signal Processing, ASSP-33, pp. 527-536, 1985.
[8] E. Dudgeon, Fundamentals of digital array processing, Proc. of the IEEE, vol 65, no 6, pp. 898904, 1977. [9] P.S. Naidu, Sensor Array Signal Processing, CRC Press LLC, 2001. [10] S. Blackman and R. Popoli, Design and Analysis of Modern Tracking Systems, Boston, MA, Artech, pp. 157-160, 1999. [11] A. S. Gelb, Applied Optimal Estimation, Cambridge, MA, MIT Press, 1974. [12] Y. Bar-Shalom and T. E. Fortmann, Tracking and Data Association, Orlando, FL, Academic Press, 1988. [13] K. Murphy, Homepage for Matlab Kalman filter toolbox, http://www.cs.ubc.ca/ murphyk/Software/Kalman/kalman.html. [14] L. Crona, E. Dalberg, R. Lennartsson, B. Lundqvist, Per M´ oren, L. Persson, P. S¨ oderberg, An experiment on passive multi-sensor underwater surveillance, FOI-R–1336–SE, 2004.
