Time-varying Array Shape Estimation by Mapping Acoustic Field Directionality Jonathan L. Odom and Jeffrey L. Krolik Dept. of Electrical and Computer Engineering Duke Univ., PO Box 90291 Durham, NC 27708 [email protected], [email protected] Abstract—This paper introduces a towed-array shape estimation technique that exploits the directional structure of the time-varying acoustic field. Unlike conventional array shape estimation methods that use discrete sources of opportunity, the proposed approach does not assume knowledge of the number of sources in the field or their estimated directions. Instead, the entire time-varying field directionality map is used. Additionally, maneuverability of the array is exploited to improve endfire resolution and left/right discrimination for a nominally linear array. The algorithm forms an approximate joint maximumlikelihood estimate of time-varying field directionality and array shape using an iterative Expectation-Maximization (EM) approach. Simulations are given to evaluate the array shape estimation error during a maneuver. In a simulated multi-source scenario, the proposed method is shown to be more robust than methods that rely on direction-of-arrival estimation when the full field around the array is considered.

I. I NTRODUCTION Passive towed-array sonar target detection and localization algorithms typically assume the array shape is known. However, during tow platform maneuvers the array necessarily deviates from its nominally linear shape. Traditionally, organic positioning systems (e.g. inertial navigation and heading sensors) are used with dynamic motion models to approximate the array shape, but these methods do not have sufficient accuracy in all cases. Small diameter towed arrays often have only a few heading sensors and are subject to frequent failure. To account for sensor failure, dynamic models are used to filter sensor outputs, but current models are unable to predict ocean currents that can dramatically alter array shape. When a tow cable is incorrectly balanced or the towing platform dramatically slows, depth of hydrophones may vary along the length of the array. Different sections of the array may experience different currents as a function of depth [1]. While position measurement systems can provide accurate shape knowledge in many cases, they are often insufficient to facilitate accurate localization and tracking of targets during maneuvers. Maneuvers provide information from a wide range of orientations and can thus potentially be exploited to increase coverage area. The location of best resolution, broadside, may sweep through a large bearing space but array shape knowledge is required to avoid significant performance loss. This paper thus considers the joint estimation of array shape and source parameters. c 2012 IEEE to appear in IEEE OCEANS 2012

Traditionally, passive acoustic-based array shape estimators model the unknown field in terms of a few discrete sources at unknown locations [2]–[4]. Typically the sources are assumed to exist over bearings from 0◦ to 180◦ or smaller. This assumption will result in array shape estimation errors due to left/right ambiguities when sources exist over a full 360◦ field. In this paper, the entire field directionality map is used to estimate array shape instead of a few discrete source locations in a restricted bearing space. The sources of opportunities exploited through the field directionality map are assumed to form a directional noise field (such as those created by distant merchant shipping lanes). Multiple distributed sources are approximated by a grid of point sources covering the entire field. This avoids the problem of correctly identifying the number of sources in the field and mitigates left/right ambiguities by extending the bearing space used for array shape estimation to the full 0◦ to 360◦ . Near-field sources are not considered here, although they may compound the problem due to the ambiguity between array shape curvature and wavefront curvature. In this paper, maximum-likelihood estimation of the entire time-varying spatial spectrum is performed using a variation of the Expectation-Maximization (EM) algorithm [5]. Since parametrization of arbitrary array shapes presents a difficult problem [6], this paper uses a piecewise linear towed array shape model in addition to a tow point indicator (TPI) [7]. The piecewise linear model parametrizes the array shape in terms of headings between rigid array segments. The absolute array orientation can be modeled using a coordinate system rotation. II. R ECEIVED S IGNAL M ODEL Consider an array of M acoustic pressure elements receiving signals from S narrowband sources. The data received by the array at time n is assumed to be given by pn = A(µn )sn + ηn , where the vector of S sources is denoted sn , the sensor noise is denoted by ηn , and the array manifold matrix,   A(µn ) = a1 (µn ) a2 (µn ) · · · aS (µn ) ,

is a function of the array shape parameters, µn . The source signals are assumed to be stochastic and circularly symmetric with complex normal distribution sn ∼ CN (0, Σn ), and the sensor noise is given by ηn ∼ CN (0, ση2 I). The

received data is thus assumed to be distributed according to pn ∼ CN (0, Rn ) with covariance matrix Rn = A(µn )Σn AH (µn ) + ση2 I. The array shape is assumed to follow a water pulley model that is valid for turns at slow speeds. The array shape is described in terms of the headings, µn , between rigid array segments connecting hydrophone elements [8]. The linear model describing the array shape is given by µn = Fµn−1 + un + ǫn

The stiffness of the array is modeled by ρ, which is the fraction of a spatial segment that the tow motion travels down the array during a single time discretization. The linear water pulley model assumes that the own ship forcing function has a wavelength that is longer than the array length. The number of array parameters is reduced using principle component analysis, as suggested by Varadarajan and Krolik [9] such that the headings can be written in terms of a smaller basis set, µ = Mγ. A total of 2 parameters, γ, are used to describe the array shape by first simulating a maneuver and then preforming PCA on the resulting headings. These parameters approximately correspond to angle and curvature as described by Gertstoft et al. [10]. The principle components of the headings for a 30 element array with 1 m inter-element spacing preforming a ±20◦ sinusoidal maneuver are shown in Figure 1. The shapes are physically intuitive as they are similar to the modes of a string. The small scale of the y-axis in Figure 1 is due to the unit length normalization of singular vectors. 0.04 0.02

y position, m

0 −0.02 −0.04

−0.08 −0.1 0

Fig. 1.

mode 1 mode 2 mode 3 mode 4 mode 5 mode 6 mode 7 5

10

15 20 x position, m

25

The array shape and field directionality map are jointly estimated using a maximum-likelihood (ML) approach. This is accomplished by first estimating the array shape using the previous field directionality map and the current acoustic data, as shown in Figure 2. The field directionality map is

(1)

where the own ship forcing function arriving at the first element of the array is denoted, un = ρ [µTPI , 0, · · · , 0], process white noise is denoted, ǫn , and displacement travels down the array according to   0 0 0 0  1 0 0 0    F = (1 − ρ)IM−1 + ρ  .  0 ... 0 0  0 0 1 0

−0.06

III. J OINT A RRAY S HAPE AND F IELD D IRECTIONALITY M AP E STIMATE

30

35

Principle components of sensor headings as Cartesian positions

Fig. 2.

Block diagram of joint array and source parameter estimation

then updated using the new array shape estimate. A single iteration occurs with each data measurement. The estimate can be replaced by a straight line assumption or heading sensor estimates when the maneuver is complete. The field directionality map estimates power as a function of bearing, formed by assuming a grid of Q sources uniformly spaced in bearing. For uncorrelated sources, the field directionality map provides an estimate of the diagonal signal covariance matrix, Σn . The number of grid points is much larger than the number of true sources such that Q > S. The steering matrix corresponding to the Q grid points is denoted by V. Using the algorithm first described by Rogers and Krolik [5], the ML estimate of the field directionality map is given by ! N X 1 old new old ˆ ˆn = Σ ˆ − Σ ˆ old , Σ G(n) Σ (2) N k=n−N +1

where the snapshot dependent terms are   ˆ −1 )V(ˆ G(n) = VH (ˆ γn )(K−1 − K−1 RK γn ) ˆ old VH (ˆ K(n) = V(ˆ γn )Σ γn ) + ση2 I,

and the previous NPdata samples form the covariance matrix H ˆ = 1 n estimate, R k=n−N +1 xk xk . The field directionality N map estimate is written in terms of the array shape estimate. Previous work assumed the array shape was known during maneuvers [5]. In this paper, the array shape is also estimated using an ML approach. The two-parameter basis set reduces the computational complexity of numerical optimization. The loglikelihood of the array shape parameters are given by   ˆ n−1 VH (γn ) + ση2 I L(γn ) = −N ln det V(γn )Σ o n ˆ n−1 VH (γn ) + ση2 I)−1 R ˆn (3) − tr (V(γn )Σ

where the true signal covariance matrix is approximated by the ˆ n−1 . When the tow field directionality map estimate, Σn ≈ Σ platform heading is known, the array shape estimate can be constrained [9]. For a single tow point indicator, the constraint ˆ = µTPI , where m1 is a column vector of is defined by mT1 γ

the first row of the basis matrix M. Using the likelihood from (3), the estimate is defined as ˆn = arg max L(γn ) sub. to mT1 γn = c γ (4)

γn

ˆ n = Mˆ µ γn . The estimate finds the best array shape that fits the acoustic data from N snapshots and the previous field directionality map. Using an eigenvector expansion of the covariance matrix estimate containing the previous field directionality map and array shape, let the qth eigenvector be denoted eq (µn ) and the qth eigenvalue be denoted λq . For geometric insight into the estimate, consider the second term of (3) written as 1 N

Q X

n X

H 2 λ−1 q |xk eq (γ)| .

(5)

k=n−N +1 q=1

The eigenvalue term places higher weight on noise eigenvectors than source eigenvectors. In this way, the second term of (3) matches the estimate to the acoustic data. Unlike adaptive beamforming techniques that place nulls in the direction of interferers, the minimization of (5) results in pushing the eigenvectors to locations that are orthogonal to the noise eigenvectors. For half-wavelength spaced linear arrays with fewer sources than elements, the number of source eigenvectors is equal to the number of sources. Additionally, the eigenvectors, eq (γ), orthogonal to the noise subspace are parallel to the signal subspace. As both the noise and signal subspaces are considered, the array shape estimate fits both the nulls and peaks in the field directionality map. It should be noted that this method is not derived as a subspace-based algorithm. The extension for broadband data, assuming uncorrelated and independent frequency bins, results in a field directionality map as described in [11] and a broadband extension to (3). The frequency dependence of the narrowband likelihood is denoted by expressing it as a probability density function, L(γn ) = ˆ n−1 (ωb )), where the recieved data matrix of f (Pn (ωb )|γn , Σ the previous N samples is given by,  Pn = pn

pn−1

···

 pn−N +1 .

The broadband estimate results in a sum of narrowband loglikelihoods with estimate, ˆn = arg max γ γn

B X

ˆ n−1 (ωb )). f (Pn (ωb )|γn , Σ

b=1

An important feature of both the FDM and array shape estimates is that no inversion of the sample covariance matrix is required. As an array preforms a maneuver, the impinging wavesfronts may change rapidly and the sample covariance matrix will vary over a large of number of snapshots. For long arrays, the number of snapshots that the received signal covariance matrix is constant may be less than the total number of elements.

IV. R ESULTS A simulation is used to compare the performance of the proposed method, referred to as ML-FDM, to the traditional array shape estimate based on direction of arrival estimation, referred to as ML-DOA [2]. In order to demonstrate the ability to distinguish sources across the entire 360◦ around the platform, consider the scenario where a source travels through endfire (105◦ to 84◦ relative to North) while another source transitions from 60◦ to 15◦ relative to North. The broadside of the array is initially pointed towards North/South. The signal-to-noise ratio for each source before array gain is 10 dB, and the noise of the array shape in the water pulley model is -40 dB, ǫ ∼ CN (0, σǫ2 I). The shape of the 30 element array with half-wavelength spacing (1 m assuming sound speed of 1500 m/s) is assumed known for the first 40 seconds then estimated for the remaining time. The tow platform performs a ±20◦ sinusoidal maneuver throughout the entire simulation assuming constant forward motion of 2 m/s. The maneuver is first simulated with a 10 Hz update rate to determine the array shape, then the headings of each array segment between hydrophones are used to form the basis set. The final simulation uses the basis set to calculate array shapes and is calculated with a 2 Hz update rate. The array shape is assumed to be constant between shape updates. A single snapshot of acoustic data collected for each array shape and a sliding window of 10 snapshots with 90% overlap is used to calculate estimates. Adaptive beamforming is snapshot limited during maneuvers, which is a significant issue [10]. The traditional technique, referred to as ML-DOA, used for array shape calibration is to treat the source power and direction of arrivals as nuisance parameters and solve for the array shape parameters. However, ambiguities due to the linear array geometry cause incorrect direction of arrival estimates and degrade array shape performance for the ML-DOA technique. This can be seen in the example shown in Figure 3, where the ML-DOA method fails between 150 and 250 seconds. The ML-FDM method is able to maintain an unambiguous field directionality map during the entire maneuver. A Monte Carlo simulation of 100 runs with different signal realizations is used to demonstrate the estimate sensitivity to the received signals with limited snapshot support. The performance is computed by calculating root mean square error of the hydrophones positions and averaging across realizations. The ML-DOA method does not always fail at the same time instance, but the average error of the array shape is higher using ML-DOA as compared to ML-FDM over the entire maneuver, shown in Figure 4. The main reason for the performance difference is that the MLDOA algorithm is sensitive to source localization accuracy, which results in large variation in error performance. The standard deviation of the error is high for ML-DOA, as shown in Figure 5. The overall error in both methods will rise over time as estimates stray from the local ML solutions. The MLDOA method provides an estimate with lower error when the direction of arrival estimates are highly accurate for both

9 8 7 Std. dev., m

sources. However, the large number of times the estimates are incorrect results in poor array shape estimation. This is shown in the empirical cumulative distribution function (CDF) of error across time and realizations combined. The the probability of low error for ML-DOA is higher than MLFDM below the median, where F(x)=.5 in Figure 6. However, above the median the ML-FDM shows a significantly higher probability over ML-DOA for the same error threshold. For example, the probability that the RMSE is at most 5 m is 0.9 for ML-FDM but only 0.6 for ML-DOA.

6 5 4 3 2

20

1

0

0 0

ML−DOA ML−FDM 50

100

Heading, degrees

−20

150 200 Time, sec

250

300

−40

Fig. 5. shape

True −60

ML−FDM

Standard deviation of Root Mean Square Error (RMSE) of array

ML−DOA −80

1 −100

0.9 −120

0.8 50

Fig. 3.

100

150 Time, sec

200

250

0.7

300

0.6 F(x)

−140 0

Array heading estimate for sensor 15

0.5 0.4 0.3 0.2

12 10

ML−DOA ML−FDM

ML−DOA ML−FDM

0.1 0 0

5

10

15 20 x, RMSE m

25

30

Mean, m

8 Fig. 6. Empirical cumlative distribution function of Root Mean Square Error (RMSE) of array shape

6 4

the estimate will be extended to consider severe turns and physical shape models that do not require prior heading knowledge.

2 0 0

Fig. 4.

50

100

150 200 Time, sec

250

300

Average of Root Mean Square Error (RMSE) of array shape

ACKNOWLEDGMENT This work was supported by ONR under grant N0001412-1-0053. Bruce Newhall at JHU/APL contributed insightful discussion to this work. R EFERENCES

V. C ONCLUSION The ML-FDM method in this paper has been used to preform towed array shape estimation directly from the field directionality estimate without source direction finding. Additionally, left/right ambiguities were mitigated to increase performance by exploiting array maneuverability. In the future,

[1] I. Schurman, “Reverberation rejection with a dual-line towed array,” IEEE J. of Ocean. Eng., vol. 21, no. 2, pp. 193–204, 1996. [2] A. Weiss and B. Friedlander, “Array shape calibration using sources in unknown locations-a maximum likelihood approach,” IEEE Trans. Acous. Speech, Signal Process., vol. 37, no. 12, pp. 1958–1966, 1989. [3] J. Goldberg, “Joint direction-of-arrival and array-shape tracking for multiple moving targets,” IEEE J. of Ocean. Eng., vol. 23, no. 2, pp. 118–126, 1998.

[4] B. Flanagan and K. Bell, “Array self-calibration with large sensor position errors,” Signal Processing, vol. 81, no. 10, pp. 2201–2214, 2001. [5] J. Rogers and J. Krolik, “Time-varying spatial spectrum estimation with a maneuverable towed array,” J. Acoust. Soc. Am., vol. 128, no. 6, pp. 3543–3553, 2010. [6] Y. Rockah and P. Schultheiss, “Array shape calibration using sources in unknown locations–Part I: Far-field sources,” IEEE Trans. Acous., Speech, Signal Process., vol. 35, no. 3, pp. 286–299, 1987. [7] J. Smith, Y. Leung, and A. Cantoni, “The Cram´er-Rao lower bound for towed array shape estimation with a single source,” IEEE Trans. Signal Process., vol. 44, no. 4, pp. 1033–1036, 1996. [8] D. Gray, B. Anderson, and R. Bitmead, “Towed array shape estimation using Kalman filters-theoretical models,” IEEE J. of Ocean. Eng., vol. 18, no. 4, pp. 543–556, 1993. [9] V. Varadarajan and J. Krolik, “Array shape estimation tracking using active sonar reverberation,” IEEE Trans. Aerosp. Electron. Syst., vol. 40, no. 3, pp. 1073–1086, 2004. [10] P. Gerstoft, W. Hodgkiss, W. Kuperman, H. Song, M. Siderius, and P. Nielsen, “Adaptive beamforming of a towed array during a turn,” IEEE J. of Ocean. Eng., vol. 28, no. 1, pp. 44–54, 2003. [11] J. Odom and J. Krolik, “Broadband field directionality mapping with spatially-aliased arrays,” in IEEE Sensor and Multichannel Signal Proc. Workshop, Jun. 2012, pp. 449–452.