A real-time array calibration method for underwater acoustic flexible sensor array

Applied Acoustics 97 (2015) 54–64

Contents lists available at ScienceDirect

Applied Acoustics journal homepage: www.elsevier.com/locate/apacoust

A real-time array calibration method for underwater acoustic flexible sensor array Liang An ⇑, Lijun Chen Key Laboratory of Underwater Acoustic Signal Processing of Ministry of Education, Southeast University, 210096 Nanjing, China

a r t i c l e

i n f o

Article history: Received 17 September 2014 Received in revised form 27 February 2015 Accepted 10 April 2015

Keywords: Array calibration Flexible array Passive ranging Time difference of arrival

a b s t r a c t To solve the precise passive ranging problem of underwater acoustic localization systems, we design a large aperture flexible sensor array whose sensors are deployed without rigid connections. We propose a novel real-time array calibration method to estimate the 3-D coordinates of the sensors. Spherical interpolation estimator is used to estimate the range of the acoustic source. The impact of the estimation errors in time-delay and sensor coordinates on the spherical interpolation range estimator is analyzed via computer simulations. By comparing the passive acoustic ranging results from sea trials with the GPS measurement data, we show that the proposed method can calibrate the flexible array in real time more accurately and estimate the underwater acoustic source range precisely. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction TDOA (Time Difference of Arrival) based methods have been widely used in passive localization of electromagnetic [1–3] and acoustic sources [4–7]. By using the TDOAs of underwater sound between hydrophones [4,8], they are shown to be particularly efficient in underwater acoustic passive localization to estimate source positions. The source range estimation is a key issue in underwater acoustic localization. Specially, Smith and Abel [9] presented a closedform least square approximate maximum likelihood source location method using TDOA measurements, the Spherical Interpolation (SI) estimator. This method uses range difference measurements taken from a fixed array to estimate source locations passively and was obtained from the ‘‘equation-error’’ minimization. The SI method had variance approaching the Cramèr-Rao lower bound (CRLB). However, in the SI method, the source ranging estimation accuracy depends on the time delay estimation (TDE) errors and the array aperture (the distance between two sensors). Therefore, higher precision of TDOA estimation and larger array aperture are needed to improve the ranging accuracy. The lower bounds of passive TDE were given by Weinstein and Weiss [10,11]. For wide-band signals, the TDE error variance depends on three main factors: the signal bandwidth, the observation time and the signal-to-noise ratio (SNR). The generalized cross-correlation (GCC) is the most commonly used method for ⇑ Corresponding author. Tel./fax: +86 2583794074. E-mail address: [email protected] (L. An). http://dx.doi.org/10.1016/j.apacoust.2015.04.008 0003-682X/Ó 2015 Elsevier Ltd. All rights reserved.

TDE, which, in some situations, could achieve the CRLB [12]. Certain methods were presented to improve the performance of GCC, such as multichannel cross-correlation-coefficient (MCCC) method [13,14], cross-correlation functions discrimination [15]. However, the precision of TDE is still limited by the three factors above in a basic sense. In the case of limited signal bandwidth, observation time and SNR, enlarging the array aperture is another efficient way to improve the ranging accuracy. Bottom-mounted arrays consisting of several hydrophones are one kind of large aperture arrays. However, hydrophones deployment and remote signal transmission are quite difficult for bottom-mounted arrays. A shipboard system is another choice for large aperture arrays. The hydrophone array could be deployed around the measurement ship. If the arrays with rigid connections are used like many traditional underwater acoustic arrays, the deployment problem would also occur for large arrays with an aperture of several decameters. Furthermore, the array cannot be anchored on a floating ship. In this paper, we propose a new passive ranging method using a ship-based flexible array. Instead of being mounted on a rigid frame, the receiving array is hung around the ship board randomly for the convenience of deployment in the sea. An effective method is presented to estimate the 3-D coordinates of ship-based flexible array sensors precisely. Our method extends the SI estimator from stationary array to non-stationary case. This paper was organized as follows: In Section 2, the new passive ranging method is briefly outlined and its estimation performance is evaluated by simulation results. The flexible array sensor position calibration method was also described in details

55

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

in Section 2. The performance of the proposed method is demonstrated through localization experiments in Section 3. Finally, conclusions are drawn in Section 4 along with discussion of future work.

Projector T

2. Methods 2.1. Flexible array deployment method The schematic diagram of the ship-based flexible array consisting of six sensors is shown in Fig. 1, where R0, R3 and R4 are hydrophones and T1/R1, T2/R2 and T5/R5 are composite underwater acoustic transducers. All the sensors are deployed around the ship board. Let R0 be the original of a coordinate system. Xi(xi, yi, zi) denotes the coordinate of sensor i for i = 1, . . ., 6. The coordinate of acoustic source S is Xs(xs, ys, zs). Rs = ||Xs||2 is the range between S and R0. If only Rs is needed to estimate using the SI estimator, the coordinate system could be any one of the rectangular coordinate systems which have the same original of R0 [9]. The distance between S and Ri is given by Di = ||Xi  Xs||2, while the distance from the origin to point Xi is given by Ri = ||Xi||2. The TDOA sij between sensor i and j is equal to the range difference (RD) dij divided by the sound speed c

sij ¼ dij =c ¼ ðDi  Dj Þ=c:

ð1Þ

D

Hydrophone R Fig. 2. Schematic diagram of the composite underwater acoustic transducer.

We also have

di0 ¼ Di  D0 :

ð4Þ

From Eqs. (3) and (4), we obtain:

Di ¼ di0 þ Rs :

ð5Þ

If the five TDOAs (excluding i = j, and counting each dij  dji) and the six sensor coordinates are known, Xs and Rs could be decided. The performance of the estimation depends on sij, Ri, and Xi. In the next section, the effect of Ri and Xi on the range estimation would be thoroughly discussed. To acquire accurate sensor coordinates, three composite underwater acoustic transducers are used to achieve the self-calibration of the flexible array. The composite underwater acoustic transducer is shown in Fig. 2, which consists of a projector T and a hydrophone R. D denotes the distance between T and R. In Fig. 1, T1/R1, T2/R2 and T5/R5 are composite transducers. A real-time array calibration method will be discussed in Section 2.3.

h i 1 T d I  SðST SÞ ST d 1 ^s ¼ h i R 2 dT I  SðST SÞ1 ST d

2.2. Spherical interpolation estimator

where

Substituting di0 + Rs for Di in Eq. (2) yields

R2i

2

 di0  2Rs di0  2X Ti X s ¼ 0;

i ¼ 1; 2; 3; 4; 5:

If the RD measurement di0 is available, the source range Rs can be estimated by solving Eq. (6). By introducing the so-called equation error and using the least square solutions in Ref. [9], the nonlinear equations in Eq. (6) could be linearized and the closed form solutions of the least square solution is given by

T

The SI estimator presented by Smith and Abel [9] is suitable for passive localization problem of arbitrary array shape. Given the definitions in Section 2.1, we have the following basic relations:

D2i ¼ kX i k22  2hX i ; X s i þ kX s k22 ¼ R2i  2X Ti X s þ R2s ;

i ¼ 1; 2; 3; 4; 5:

ð6Þ

ð7Þ

d ¼ ½d10 ; d20 ; d30 ; d40 ; d50 

ð8Þ

S ¼ ½X 1 ; X 2 ; X 3 ; X 4 ; X 5 T

ð9Þ

h iT 2 2 2 2 2 d ¼ R21  d10 ; R22  d20 ; R23  d30 ; R24  d40 ; R25  d50

ð10Þ

ð2Þ Since R0 is the original of the coordinate, D0 could be written as:

D 0 ¼ Rs :

ð3Þ

T2/R2 Sea surface T1/R1

R0

Flexible array

R3

T5 /R5

S

R4 Sea bottom

Acoustic source

Fig. 1. Schematic diagram of the underwater acoustic flexible array.

and I is the identity matrix. If the sensors coordinates are given as X1(Da, 0, 0), X2(0, Da, 0), X3(Da, 0, 0), X4(0, Da, 0), X5(0, 0, Ha), the range estimation performance could be analyzed by Monte Carlo method. Suppose the time delay estimation follows the Gaussian distribution N(si0, r2s), the sensor coordinates estimations follow the Gaussian distributions N(xi0, r2x), N(yi0, r2y) and N(zi0, r2z), respectively. si0 denotes the true value of si, so do xi0, yi0, zi0 for xi, yi, zi. rs is the root-meansquare error (RMSE) of the time delay estimation error, and rx, ry, rz are the RMSEs of the x, y, z-coordinate estimation errors. Rsh is defined as the projection of Rs in xoy plane. Figs. 3–5 give the experimental results of 2000-trial Monte-Carlo runs. The source coordinate in Fig. 3 is Xs(259.81, 150.00, 15.00) and Rsh equals 300 m. Fig. 3(a) and (b) shows the RMSE of range estimation versus different TDOA estimation errors. Fig. 3(c) and (d) shows the RMSE of range estimation versus different sensors coordinates estimation errors. In Fig. 3(a) and (c) Ha is equal to 10 m, and in Fig. 3(b) and (d) Da is equal to 20 m. The parameters in Figs. 4 and 5 are the same as those in Fig. 3, except that the source coordinate is Xs(346.41, 200.00, 15.00) and

56

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

(b)

(a) 10

18 D =15m a

16

H =5m

9

a

D =20m

H =10m

RMSE of Range Estimation (m)

RMSE of Range Estimation (m)

a

D =25m

14

a

12 10 8 6 4 2 0

a

8

H =15m a

7 6 5 4 3 2 1

0

5

10

0

15

0

5

RMSE of Time Delay Estimation ( μs)

(c)

15

(d) 40

80

H =5m

D =15m

a

a

70

35

D =20m

H =10m a

RMSE of Range Estimation (m)

a

RMSE of Range Estimation (m)

10

RMSE of Time Delay Estimation ( μs)

D =25m a

60 50 40 30 20 10

H =15m a

30 25 20 15 10 5

0 0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0

0.16

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

RMSE of Sensors Coordinates (m)

RMSE of Sensors Coordinates (m)

Fig. 3. Range estimation errors of SI method (Rsh = 300 m).

Xs(433.01, 250.00, 15.00) respectively. The corresponding Rsh is 400 m and 500 m. The simulation results in Figs. 3–5 show that array aperture, TDOA and sensor coordinates estimations are all critical factors for SI method. Larger array aperture could lessen the effect of the TDOA and sensor coordinate errors on the range estimation accuracy. Furthermore, the horizontal aperture scale is more important than the vertical aperture scale in underwater range estimation. Comparing the results in Fig. 3 with those in Figs. 4 and 5, we can find that the range estimation errors increase with the increase of source range. However, the growth rate of errors increase decreases with the increase of array apertures. 2.3. Real-time array shape calibration method From Section 2.2, the SI estimator is very sensitive to sensor coordinates error. Since the sensors of flexible array are not stable, an effective array calibration method is required to measure the sensors’ coordinates in real-time.

Section 2.1, the source range Rs is independent of the adopted coordinate and the only limiting condition is that R0 must be regarded as the origin of a coordinate system. Therefore, calculation of the sensors’ coordinates could be simplified by choosing a proper coordinate system. Let X0, X1 and X2 define the xoy plane. X0 and X1 define the x-axis with X1 lying on the negative direction of the x-axis. The y-axis and x-axis intersect at X0 and are mutually perpendicular. The z-axis goes through X0 and is perpendicular to the xoy plane. Then we have y1 = z1 = z2 = 0. The number of unknown coordinate values reduces to twelve and only twelve equations are needed. The coordinate system is shown in Fig. 6. Define the distance between Xi and Xj as rij, i = 0, 1, 2, 3, 4, 5, j = 0, 1, 2, 3, 4, 5. There are totally fifteen distance values between any two sensors of the array. In this paper, twelve of them, i.e., r01, r12, r31, r41, r02, r32, r42, r05, r15, r25, r35 and r45, are chosen to form the sensor coordinate equations to determine the sensor coordinates Xi(xi, yi, zi) for i = 1, 2, . . ., 5, as shown in Eqs. (11–14).

x1 ¼ r 01 2.3.1. Calculation of sensors’ coordinates In Eq. (7), fifteen coordinate values of X1 to X5 must be known to calculate Rs. A sophisticated method is to establish fifteen different equations about the fifteen coordinate values. As explained in

(

x22 þ y22 ¼ r 202 ðx2  x1 Þ2 þ y22 ¼ r 212

ð11Þ ð12Þ

57

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

(a)

(b)

30

18

D =15m a

a

H =10m

a

RMSE of Range Estimation (m)

RMSE of Range Estimation (m)

H =5m

16

D =20m

25

D =25m a

20

15

10

a

14

H =15m a

12 10 8 6 4

5 2

0

0

0

5

10

15

0

5

10

(c)

(d)

140

70

D =15m

H =5m

a

a

60

D =20m a

RMSE of Range Estimation (m)

RMSE of Range Estimation (m)

120

D =25m a

100

80

60

40

H =10m a

H =15m a

50

40

30

20

10

20

0

15

RMSE of Time Delay Estimation ( μs)

RMSE of Time Delay Estimation ( μs)

0

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

RMSE of Sensors Corrdinates (m)

RMSE of Sensors Corrdinates (m)

Fig. 4. Range estimation errors of SI method (Rsh = 400 m).

8 2 x þ y25 þ z25 ¼ r205 > > < 5 ðx5  x1 Þ2 þ y25 þ z25 ¼ r 215 > > : ðx5  x2 Þ2 þ ðy5  y2 Þ2 þ z25 ¼ r 225 8 ðxk  x5 Þ2 þ ðyk  y5 Þ2 þ ðzk  z5 Þ2 ¼ r 2k5 > > < ; ðxk  x1 Þ2 þ y2k þ z2k ¼ r2k1 > > : 2 2 ðxk  x2 Þ þ ðyk  y2 Þ þ z2k ¼ r 2k2

where

ð13Þ

k ¼ 3; 4

ð14Þ

y1 ¼ 0;

x2 ¼ r 02 cos a;

y ¼ r 02 sin a;

r 205  r 215 þ x21 ; 2x1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼  r 205  x25  y25

x5 ¼ 

xk ¼

Gk þ

z1 ¼ 0

y5 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G2k  4F k Hk 2F k

z2 ¼ 0

r 205  r 225 þ x22 þ y22 x5 x2  ; 2y2 y2

yk ¼ Ak  xk Bk ;

F k ¼ 1 þ B2k þ D2k ; x22

A2k

y1  y5 ; z5

Ek ¼

2Ak ðy1  y5 Þ  C k 2z5

Gk ¼ 2ðy2 Bk  Ak Bk þ Dk Ek  x2 Þ;

 2Ak y2 þ y22 þ E2k  r 22 :

Hk ¼

ð16Þ

2.3.2. Measurement of distances between sensors In underwater acoustic signal processing, distance values are often measured indirectly by TOA (Time of Arrival) method. In this paper, we assume that the velocity of sound is constant or changes slowly with depth and range. Given the velocity of sound in water c, if the signal propagation time s from the source to the receiver is measured, the distance between the source and the receiver could be calculated as

z5

zk ¼ Dk xk þ Ek

r 2k1  r2k2  x21 þ x22 þ y22 ; 2y2

ð15Þ

ð17Þ

;

Ak ¼

C k ¼ r 2k5  r 2k1  x25  y25  z25 þ x21 þ y21 ; Dk ¼ x1  x5  Bk

By solving Eqs. (11–14), we are able to obtain the sensor coordinates Xi(xi, yi, zi) for i = 1, 2, . . ., 5 in the following equations:

x1 ¼ r 01 ;

r 201 þ r202  r 221 ; 2r 01 r 02 x1 þ x2 Bk ¼ y2

a ¼ arccos

ð18Þ

r ¼ cs

þ

ð19Þ

58

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

(a)

(b)

50

30

D =15m

H =5m

a

45

a

D =20m 40

RMSE of Range Estimation (m)

RMSE of Range Estimation (m)

a

D =25m a

35 30 25 20 15 10

H =10m

25

a

H =15m a

20

15

10

5

5 0

0

5 10 RMSE of Time Delay Estimation ( μs)

0

15

0

5

10

(c)

(d)

300

120

D =15m

H =5m

a a

H =10m

100

RMSE of Range Estimation (m)

RMSE of Range Estimation (m)

a

D =20m

250

D =25m a

200

150

100

50

0

15

RMSE of Time Delay Estimation ( μ s)

a

H =15m a

80

60

40

20

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

RMSE of Sensors Corrdinates (m)

RMSE of Sensors Corrdinates (m)

Fig. 5. Range estimation errors of SI method (Rsh = 500 m).

y T2/R2 T1/R1

z

R0 R3 T5/R5

x

R4 Fig. 6. Coordinate system of sensors array.

In this paper, the distances between array sensors were measured with the composite underwater acoustic transducers described in Section 2.1. As shown in Fig. 2, the projector T emits special array calibration signal and the hydrophone R receives the calibration signal. Based on this, the three projectors in the three composite transducers are used to emit three calibration signals respectively. Two measures are adopted to distinguish the three calibration signals and reduce the mutual interference. First, three Direct Sequence Spread Spectrum–Differential Phase Shift Keying (DSSS–DPSK) pulse signals are used as the calibration

signals to reduce their correlations. Since the calibration signal is emitted while the acoustic source radiating noise, the frequency band of calibration signal should not overlap that of radiated noise. In this paper, three m sequences [16], which are widely used in spread spectrum communication systems, are adopted as the modulating signals, called m1, m2 and m3, generated by 7th polynomials. Each m sequence has 128 binary codes including one reference code. m1 and m2 are preferred pair [17], so do m3 and m2. The length of one binary code is 0.1 ms and the length of the m sequence is 12.8 ms accordingly. The carrier signal s0 is a 60 kHz sinusoidal signal. Its phase is modulated by the three m sequences respectively with DPSK mode and the output signals are named s1, s2 and s3. The spectrum of s1 is shown in Fig. 7. The normalized auto-correlation function of s1 and the cross-correlation functions between s1, s2 and s3 are shown in Figs. 8–10. Second, the three projectors emit calibration signals orderly. In Fig. 2, the hydrophone is close to the projector. The received calibration signal intensity of a hydrophone projected by the projector just above itself would be higher than those projected by other two projectors. Therefore, the three calibration signals must be projected in a time-sharing mode to avoid signals superimposing in time domain. A simulated flexible array schematic diagram is shown in Fig. 11. Suppose that the coordinates of the sensors in Fig. 11 are X0(0, 0, 0), X1(Da, 0, 0), X2(0, Da, 0), X3(Da, 0, 0), X4(0, Da, 0) and X5(0, 0, Ha). The projectors emit s1, s2 and s3 successively. The time

59

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64 0.25

-10

0.2

Normalized Auto-correlation Function

0

Relative Amplitude (dB)

-20 -30 -40 -50 -60 -70 -80 -90

0

10

20

30

40

50

60

70

80

90

100

0.05 0 -0.05 -0.1

-0.005

0

0.005

0.01

0.015

The envelope sequencing of projected and received calibration signals are shown in Fig. 12. In Fig. 12(b), the time interval between the falling edge of s1 and the rising edge of s2 is 3.2 ms (the area in the circle), which is the minimum of the time intervals between two envelopes. Therefore, the sequencing of the projected signals can guarantee that all the calibration signals do not superimpose each other in time domain.

0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6

-0.01

-0.005

0

0.005

0.01

0.015

0.01

0.015

Time (s) Fig. 8. The auto-correlation function of s1.

0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2 -0.015

-0.01

Fig. 10. The cross-correlation function of s3 and s2.

1

-0.8 -0.015

-0.2 -0.015

Time (s)

Fig. 7. The spectrum of array calibration signal s1.

Normalized Auto-correlation Function

0.1

-0.15

Frequency (kHz)

Normalized Auto-correlation Function

0.15

-0.01

-0.005

0

0.005

Time (s) Fig. 9. The cross-correlation function of s1 and s2.

interval between two calibration signals is 35 ms, when Da equals to 20 m and Ha equals to 10 m.

2.3.3. Performance of the array calibration method From Eqs. (15–18), the accuracy of sensors coordinates estimation depends on the distances estimation between sensors. According to Eq. (19), the estimation accuracy of distances is proportional to that of the time-delay estimation. In other words, the time-delay estimation performance determines the sensors’ coordinates estimation performance. Since it is difficult to analyze the sensors coordinates estimation performance analytically, Monte Carlo method is also used to simulate the RMSE of sensors coordinates’ estimation. If the RMSE of time delay estimation is rt, the RMSE of sensors distance estimation is rd = crt, while c is equal to 1530 m/s and the time-delay estimation error follows the Gaussian distribution N(0, r2t). Fig. 13 gives the experimental results of 2000-trial Monte Carlo runs with different Da and a fixed Ha of 10 m. From the results, we can see that Da has no significant effects on sensors coordinates’ estimation errors. First, we consider the situation when source range is 300 m. It can be seen from Fig. 3(c) that the RMSE of range estimation (Column 2) would be less than 9 m (about 3% of the source range) if the RMSE of sensors coordinates is less than 0.02 m, 0.035 m and 0.055 m when Da equals to 15 m, 20 m and 25 m respectively. Correspondingly, the RMSE of time-delay estimation should be less than 2.5 ls, 4.5 ls and 6.0 ls respectively given the results in Fig. 13. However, the RMSE of range estimation would be less than 18 m (about 6% of the source range) if the RMSE of sensors coordinates is less than 0.04 m, 0.075 m and 0.115 m when Da equals to 15 m, 20 m and 25 m respectively. And the corresponding RMSE of time-delay estimation is 5.0 ls, 9.0 ls and 12.5 ls from Fig. 13(a). Similarly, when the source range is 400 m and 500 m, the relations between the RMSE of sensors coordinates, the RMSE of timedelay estimation and the RMSE of range estimation can be found from Figs. 4(c), 5(c) and 13. The results are shown in Tables 1–3, in which the RMSE of sensors coordinates (Column 3) and RMSE of time-delay estimation (Column 4) are the upper bound directly

60

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

20m

T2/R2

20m 20m

R3

R0

T1/R1

R0

20m

R3

20m

10m

20m

T1/R1

T5/R5 R4

(a) Top view of the flexible array

(b) Sideview ofthe flexible array

Fig. 11. Simulated schematic diagram of the flexible array.

Calibration signals received by R0

Calibration signal projected by T

Nomalized Amplitude

Nomalized Amplitude

1

Envelope of s

1

1

0.5 0 0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

1

0

0.01

0.02

Nomalized Amplitude

Nomalized Amplitude

Envelope of s 2

0.5 0 0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

Envelope of s

0.5

0.04

0.05

0.01

0.02

0.06

0.03

Nomalized Amplitude

0.07

0.08

0.09

1

0.1

0 0

0.01

0.02

0.03

Nomalized Amplitude

Envelope of s1

Envelope of s3

Envelope of s2

0 0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

0.09

0.1

Nomalized Amplitude

Calibration signals received by R4 1 Envelope of s1

Envelope of s2

Envelope of s3

0 0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

Calibration signals received by R5 1 Envelope of s1

0.5

Envelope of s3

Envelope of s2

0 0

0.1

0.09

0.1

0.09

0.1

Envelope of s

3

0.05

0.06

0.07

0.08

0.04

0.05

0.06

0.07

0.08

(b) The envelop of received signals by R 0 , R 1 and R 2

1

0.5

0.09

Envelope of s 2 Envelope of s 3

1

Calibration signals received by R3

0

0.08

TIme (s)

(a) The envelop of projected signals

0.5

0.04

Envelope of s

0.5

Time (s)

Nomalized Amplitude

Nomalized Amplitude

0 0.03

0.07

Calibration signals received by R 2

0.5 0.02

0.06

Envelope of s2

1

0

3

0.01

0.05

0

0.1

Envelope of s3

0

0.04

1

Calibration signal projected by T 1

0.03

Calibration signals received by R1

2

0

Envelope of s3

0

0.1

Calibration signal projected by T 1

Envelope of s2

Envelope of s1

0.5

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

Time (s)

(c) The envelop of received signals by R 3, R 4 and R5 Fig. 12. The envelope of projected and received calibration signals.

61

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

Table 1 The relations between range, coordinates and time-delay estimation errors (Da = 15 m).

0.12 R

RMSE of Sensors Corrdinates (m)

1

R

2

0.1

R

3

R 0.08

R

5

0.06

RMSE of range estimation (m)

RMSE of sensors coordinates (m)

RMSE of time-delay estimation (ls)

300 400 500 300 400 500

8.42 11.07 11.50 16.69 22.36 28.69

0.02 0.015 0.01 0.04 0.03 0.025

2.5 1.5 1.0 5.0 4.0 3.0

0.04

0.02

0

Table 2 The relations between range, coordinates and time-delay estimation errors (Da = 20 m). 0

5

10

15

RMSE of Time Delay Estimation ( μs)

(a) Da = 15m 0.14 R

1

RMSE of Sensors Coordinates (m)

Source range (m)

4

0.12

R

Source range (m)

RMSE of range estimation (m)

RMSE of sensors coordinates (m)

RMSE of time-delay estimation (ls)

300 400 500 300 400 500

8.12 10.21 13.05 17.86 22.96 29.58

0.035 0.025 0.02 0.075 0.055 0.045

4.0 3.0 2.0 9.0 6.5 5.0

2

R

3

0.1

R

4

R

Table 3 The relations between range, coordinates and time-delay estimation errors (Da = 25 m).

5

0.08

0.06

Source range (m)

RMSE of range estimation (m)

RMSE of sensors coordinates (m)

RMSE of time-delay estimation (ls)

0.04

300 400 500 300 400 500

8.41 11.84 14.43 17.44 22.43 28.31

0.055 0.045 0.035 0.115 0.085 0.07

6.0 4.5 3.5 12.5 9.0 7.5

0.02

0

0

5

10

15

RMSE of Time Delay Estimation ( μs)

(b) Da = 20m

3. Experimental results

0.14 R

RMSE of Sensors Corrdinates (m)

1

0.12

R

2

R

3

0.1

R

4

R

5

0.08

0.06

0.04

0.02

0

0

5

10

15

RMSE of Time Delay Estimation ( μs)

(c) Da = 25m Fig. 13. Sensors coordinates estimation errors varying with time-delay estimation errors of different Da.

read from Fig. 13 to meet the error requirements (Column 2). It can be seen from the three tables that the larger array is more robust than the smaller ones to the time-delay estimation error.

To verify the method in this paper, the SI ranging method and flexible array calibration method discussed in the above sections were realized in a real-time underwater acoustic passive localization system consisting of 6 TigerSHARC 201 Digital Signal Processors. It takes 0.2 s for the localization system to calibrate the flexible array and estimate the source range each time. The system flow chart is shown in Fig. 14, where fL is the lower-frequency limit, fH is the higher-frequency limit of the source radiated noise, fc is the carrier signal frequency of calibration signal, and Df is half of the main lobe width of the calibration signal power spectrum. The real-time underwater acoustic passive localization systems were tested in a sea trial taken on South China Sea. As shown in Fig. 15, the flexible array consisted of six sensors, where T1/R1, T2/R2 and T5/R5 were composite underwater acoustic transducers, R0, R3 and R4 were hydrophones. R0, T1/R1, T2/R2, R3 and T5/R5 were hung around the measurement ship. R4 was hung on an inflatable floating body made of PVC (Polymer of Vinyl Chloride) and was 16 m away from R0. A 7 kg lead block was hung under every flexible array sensor to maintain its stability. All the sensors depths were about 25 m except that of T5/R5, which was about 35 m. The acoustic source S was hung on the target ship with a depth of 40 m and radiated wide band noise from 4 kHz to 12 kHz. During the sea trial, the measurement ship anchored, while the target ship shut down its power devices and drifted downstream. Projectors T1, T2 and T5 projected array calibration signal as described in Section 2.3.2. Hydrophones R0, R1, R2, R3, R4 and R5

62

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

Calibration Signal

Digital Band Pass Filter fc-Δ f ~f c+Δ f

Differential Coherent Demodulation

Sensors Coordinates

Sensors Coordinates Estimation

Time-delay Estimation

SI Source Range Estimator

Received Signal

Source Radiated Noise

Digital Band Pass Filter fL ~ fH

Source Range

5 TDOAs

TDOA Estimation

Fig. 14. Signal processing flow chart of flexible passive ranging system.

GPS 2

Target Ship

PVC Floating Body (16m)

Acoustic Source S

R4

GPS1 R3

T5/R5

R0

20m

16m

20m

T1/R1

Measurement Ship T2/R2 Fig. 15. Schematic diagram of the transducers deployment in sea trial.

received the source noise and array calibration signals simultaneously. Two differential GPS (DGPS) equipments are fixed above R0 and S to record the movement locus. The filtered array calibration signal received by R3 is shown in Fig. 16.There are crosstalk signals of the calibration signal as well as the direct-path signals in the received signals because the projection subsystem and the receiving subsystem were installed in

the same electronic case. However, the amplitudes of crosstalk signals are less than that of the direct-path signals and the time-delays of crosstalk signals relative to the projection time are equal to zeroes approximately. It is convenient to discriminate the true peaks in cross-correlation functions when using cross-correlation method to estimate time-delays of the direct-path signals. The SNR of s1 in frequency band 50–70 kHz is about 24 dB, while the SNR of s2 and s3 is about 28 dB. Fig. 17 gives the twelve distances between flexible array sensors in real time measurement around 10 min. Since the flexible array sensors were not fixed, the distances between sensors fluctuated within the limit of about 2 m. The standard deviations of sensor distances estimations are shown in Table 4, the maximum of which is 9.499 mm and the minimum is about 1.650 mm. This means that maximum and minimum of the standard deviations of time delay estimations are 6.21 ls and 1.08 ls respectively. It is important to point out that the sensors fluctuations are inevitable in sea trial. Therefore, the slowly varying components of the distances were first estimated by PCA (Principal Component Analysis) method, then the residual errors between distance values and the slowly varying components are used to estimate the standard deviations, which could reflect the sensor distances estimation performance approximately. Fig. 18 gives the measurement results of the array sensor positions varying during about 10 min. It’s worth noting that R4 oscillated in the vertical direction more greatly than other sensors. The vertical motion range is about 2 m. This is because R4 was hung on a floating body, whose vertical motion with waves leaded to the oscillation of R4. The measurement results agree well with the real situation and verify the efficiency of the array calibration method. The range estimation results can further verify the efficiency of the array calibration method. Fig. 19 gives the real time range

0.2 Direct-path signal of s

0.15

Direct-path signal of s

2

1

Crosstalk signal of s

3

Amplitude (V)

0.1 Crosstalk signal of s 1 0.05

0 -0.05

-0.1

-0.15

-0.2

Crosstalk signal of s

0

0.01

0.02

0.03

0.04

2

0.05

Direct-path signal of s

0.06

0.07

Time (s) Fig. 16. The signal waveform received by R3.

0.08

0.09

3

0.1

63

r

02

(m) 31

r

0

100

200

300

400

500

600

0

100

200

300

400

500

600

21 20 19

(m) 35

r

200

300

400

500

600

r

41

(m)

100

100

200

300

400

500

600 42

r

0

100

200

300

400

500

0

100

200

300

400

600

500

600

20 19 18 22 21 20 32.5 32 31.5

(m)

0

(m)

(m) 12

r

(m) 15

r (m)

25

r

23 22 21

0

45

r

8.5 26.5 26 25.5

r

05

(m)

9

37 36 35 22.5 22 21.5

(m)

(m)

17.5 17 16.5

32

21 20 19

r

r

01

(m)

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

20 18 16

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

Time (s)

Time (s) Fig. 17. Sensors distances estimations.

Table 4 Standard deviations of the sensor distance estimations.

900

Sensor distance

Standard deviation (m)

Sensor distance

Standard deviation (m)

r01 r02 r05 r12 r15 r25

0.003974 0.003294 0.001650 0.006310 0.007863 0.004570

r31 r32 r35 r41 r42 r45

0.007878 0.004107 0.004375 0.004712 0.006226 0.009499

800

results using calibrated sensor coordinates results using DGPS reuslts using fixed sensor coordinates

Distance (m)

700 600 500 400 300 200

R

z (m)

5

3

R

100 0

300

400

500

Fig. 19. Comparison of the source range estimation with GPS measurement.

T/R1

T/R2

-5

200

Time (s)

R0

0

100

4

9

T/R

-10 20

5

8

-20 0

-10

7

10 -40

y (m)

20

Fig. 18. The array calibration results.

estimation results between acoustic source S and R0. The solid line represents the range estimation results using the calibrated sensor coordinates in Fig. 18. The dash line represents the range measurement results using the DGPS devices. The dash-dot line represents the range estimation results using fixe sensor coordinates (sensor coordinates without calibration). The comparison of the results shows that the range estimation values with calibrated sensor coordinates are stable and match the DGPS measurement results well, while the range estimation values using fixed sensor coordinates deviate significantly from the DGPS measurement results. Fig. 20 gives the range estimation error ratio using calibrated sensor coordinates. The range measured by GPS is used as the true value of source range. The error ratio is defined as the ratio between TDOA estimation error and the true value of source range.

6

Error ratio (%)

x (m)

0

-20

5 4 3 2 1 0 150

200

250

300

350

400

450

True value of source range (m) Fig. 20. The source range estimation error ratio using calibrated sensor coordinates.

When the source range is less than 300 m, the estimation error is less than 3% of the true value of the source range. The errors are ascending when the source range increases. If the source range is

64

L. An, L. Chen / Applied Acoustics 97 (2015) 54–64

less than 400 m, the estimation error is less than 6% of the true value of the source range. When the source range is larger than 400 m, the increase of TDOA estimation error of the source radiated noise signal causes the increase of the source range estimation error. The experiment results agree with the simulation results in Section 2. 4. Conclusion and discussions A large aperture flexible sensor array is used in passive underwater acoustic ranging system in this paper. To reduce the uncertainty of the sensor coordinates, a new real time array selfcalibration method is presented with composite underwater acoustic transducers. The array consists of three composite acoustic transducers and three hydrophones. Sensor coordinates are calculated with the TOAs of DPSK pulses. These coordinates are then adopted to estimate the source range. Ranging experiment results on the sea indicate that the method in this paper could give accurate real time source range. However, it could be seen from the sea trial ranging results that there exist big range estimation deviations when the source range is larger than 400 m with an array aperture of about 40 m. The first reason is the increase of TDOA estimation error of the source radiated noise, and the other important reason is the array calibration error. The acoustic centers of the projector and the hydrophone of a composite transducer are not coincident. The errors of noncoincidence are omitted in this paper, which may lead to large error in range estimation when the source is far away. More research work is needed to reduce this error. Acknowledgements The authors appreciate Professor Xiang Gao for helpful discussions and suggestions. This work was supported by National Science Fund of China under Grant Nos. 11104029 and 11104141.

References [1] Martin RK, Yan CP, Fan HH. Bounds on distributed TDOA-based localization of OFDM sources. In: IEEE int. conf. acoust. speech signal process (ICASSP 2009), Taipei; April 2009. [2] Bocquet M, Loyez C, Benlarbi-Delaï A. Using enhanced-TDOA measurement for indoor positioning. IEEE Microwave Wirel Compon Lett 2005;15(10):612–4. [3] Ho KC, Chan YT. Solution and performance analysis of geolocation by TDOA. IEEE Trans Aerospace Electron Syst 1993;29(4):1311–22. [4] Giraudet P, Glotin H. Real-time 3D tracking of whales by echo-robust precise TDOA estimates with a widely-spaced hydrophone array. Appl Acoust 2006;67:1106–17. [5] Zhao YL, Chen X, Wang B. Real-time sound source localization using hybrid framework. Appl Acoust 2013;74:1367–73. [6] Dvorkinda TG, Gannotb S. Time difference of arrival estimation of speech source in a noisy and reverberant environment. Signal Process 2005;85:177–204. [7] Lombard A, Zheng YH, Buchner H, et al. TDOA estimation for multiple sound sources in noisy and reverberant environments using broadband independent component analysis. IEEE Trans Audio Speech Lang Process 2011;19(6):1490–503. [8] Xerri B, Cavassilas JF, Borloz B. Passive tracking in underwater acoustic. Signal Process 2002;82:1067–85. [9] Smith JO, Abel JS. The spherical interpolation method of source localization. IEEE J Ocean Eng 1987;12(1):246–52. [10] Weiss AJ, Weinstein E. Fundamental limitations in passive time-delay estimation-part II: narrow-band systems. IEEE Trans Acoust Speech Signal Process 1983;31(2):472–86. [11] Weinstein E, Weiss AJ. Fundamental limitations in passive time-delay estimation-part II: wide-band systems. IEEE Trans Acoust Speech Signal Process 1984;32(5):1064–78. [12] Knapp CH, Carter GC. The generalized cross-correlation method for estimation of time delay. IEEE Trans Acoust Speech Signal Process 1976;24(4):320–7. [13] Chen J, Benesty J, Huang Y. Robust time delay estimation exploiting redundancy among multiple microphones. IEEE Trans Speech Audio Process 2003;11:549–57. [14] Benesty J, Chen J, Huang Y. Time-delay estimation via linear interpolation and cross-correlation. IEEE Trans Speech Audio Process 2004;12:509–19. [15] Wan XW, Wu ZY. Sound source localization based on discrimination of crosscorrelation functions. Appl Acoust 2013;74:28–37. [16] Proakis JG. Digital communications. 4th ed. McGraw-Hill Education; 2001. p. 766–8. [17] Dinan EH, Jabbari B. Spread codes for direct sequence CDMA and wideband CDMA cellular networks. IEEE Commun Mag 1998:48–54.