## Abstract

Coherent Doppler sonar is a useful tool for noninvasive measurement of ocean currents, sediment transport, and turbulence in coastal environments. Various methods have been proposed to separately address two of its inherent limitations: velocity ambiguity and measurement noise. However, in energetic turbulent flows, both factors may be present simultaneously. The presence of measurement noise complicates velocity ambiguity resolution, and conversely velocity ambiguity presents a challenge for existing noise suppression methods. A velocity estimator based on maximum a posteriori (MAP) estimation has been developed to resolve velocity ambiguity and suppress measurement noise simultaneously rather than separately. The estimator optimally combines measurements from multiple acoustic carrier frequencies and multiple transducers. Data fusion is achieved using a probabilistic approach, whereby measurements are combined numerically to derive a velocity likelihood function. The MAP velocity estimator is evaluated using a high-fidelity coherent Doppler sonar simulation of oscillating flow and data from a towing tank grid turbulence experiment where both velocity ambiguity and backscatter decorrelation were present. Time series and spectra from MAP velocity estimation are compared to those obtained with conventional Doppler signal processing. In addition to robustly resolving velocity ambiguity, the MAP velocity estimator is shown to reduce high-frequency noise in turbulence spectra.

## 1. Introduction

Pulse-to-pulse coherent Doppler sonar has been widely used to investigate oceanic transport and mixing processes. Applications include the measurement of tidal flows (Lhermitte 1983), surface wave breaking (Veron and Melville 1999), internal waves (Plueddemann 1992), sediment transport (Smyth et al. 2002), and near-bed turbulence (Hay 2008).

Backscatter autocorrelation is the fundamental measurement in coherent Doppler sonar. For a sequence {*z _{n}*} of backscatter samples, the autocorrelation at a lag of

*k*pulse-to-pulse intervals is

where *E* denotes expected value, * denotes complex conjugation, and *τ* is the time interval between successive acoustic transmissions. The pulse-to-pulse autocorrelation magnitude is given by

Radial velocity is proportional to the pulse-to-pulse phase shift *φ* = ∠*R*(*τ*) (Lohrmann et al. 1990),

where *c* is the speed of acoustic wave propagation and *f*_{0} is the carrier frequency of the transmitted signal. In (3), the phase interval [−*π*, *π*] corresponds to an ambiguity velocity

In other words, scatterer velocities of *υ _{a}* and −

*υ*are indistinguishable because phase shifts of

_{a}*π*and −

*π*are equivalent.

In practice, autocorrelation is estimated from a finite sequence of *M* pulse pairs (Zrnić 1977),

The phase of the autocorrelation estimate is determined from real and imaginary components of . Corresponding to each phase measurement, an estimate of autocorrelation magnitude may also be defined (Zedel et al. 1996) as

The coefficient provides a measure of data quality and satisfies .

Measurement noise is caused by pulse-to-pulse backscatter decorrelation from (i) scatterer advection through the sample volume, (ii) velocity shear and turbulence within the sample volume, (iii) phase distortion of the transmitted wave, and (iv) electronic noise in the receiver circuitry (Zedel et al. 1996; Hurther and Lemmin 2001). Noise suppression is possible when redundant measurements with uncorrelated measurement errors are available. An example of this approach is the use of two nonoverlapping sample domains (Garbini et al. 1982). Redundant measurements with uncorrelated errors may also be obtained either by using additional acoustic receivers (Rowe et al. 1986; Hurther and Lemmin 2001) or by simultaneously transmitting two or more acoustic carrier frequencies (Hurther and Lemmin 2008; Hay et al. 2008). The results in Hurther and Lemmin (2008) show that a carrier frequency shift of 10% is sufficient to produce statistically independent simultaneous velocity measurements.

Coherent Doppler sonar measurements are also limited by the existence of velocity ambiguity. Velocity ambiguity may be mitigated either by cross-correlating broadband pulses (Brumley et al. 1991) or by invoking the assumption that the velocity field is continuous in time and space, for example, by adding multiples of 2*π* to the autocorrelation phase to minimize discontinuities in time (Smyth and Hay 2003), space (Ray and Ziegler 1977), or both simultaneously (Franca and Lemmin 2006). Another method for velocity ambiguity resolution involves varying the pulse-to-pulse interval *τ* to obtain two or more staggered pulse repetition frequencies (Lhermitte and Serafin 1984; Lohrmann et al. 1990). Multifrequency coherent Doppler sonar is a recent development where multiple acoustic carrier frequencies are transmitted simultaneously using a wide bandwidth acoustic transducer (Nitzpon et al. 1995; Hay et al. 2008). The use of multiple frequencies gives rise to different ambiguity velocities from which the true velocity can be determined (Zedel and Hay 2010).

To resolve velocity ambiguity and suppress measurement noise *simultaneously*, a velocity estimator has been developed based on the Bayesian technique of maximum a posteriori (MAP) estimation (Kay 1993). The MAP velocity estimator shares many similarities with a Kalman smoother: (i) the estimator processes a time series recursively, (ii) a model is used to provide prior statistical knowledge, (iii) the estimator produces its own performance measure, and (iv) estimator lag is eliminated via application both forward and backward in time. However, unlike a Kalman smoother, the MAP velocity estimator makes use of non-Gaussian probability density functions (PDFs) and is inherently nonlinear. The motivation for using MAP estimation is that velocity ambiguity causes the velocity likelihood function to be multimodal; thus, a Gaussian representation is inappropriate.

The MAP velocity estimator optimally combines measurements from multiple acoustic carrier frequencies and multiple transducers. Data fusion is achieved using a probabilistic approach, whereby measurements are combined numerically to derive a velocity likelihood function. While the focus of this paper is on multistatic, multifrequency coherent Doppler sonar, the results are equally applicable to single-frequency and staggered pulse repetition frequency systems with monostatic or multistatic geometries. This paper also focuses on the effectiveness of the MAP velocity estimator rather than on computational or implementation issues.

Numerical simulation of oscillating flow is used to validate the algorithm by comparing velocity estimates with the known simulated time series. Results are also presented for a towing tank grid turbulence experiment where both velocity ambiguity and backscatter decorrelation were present. Time series and spectra from MAP velocity estimation are compared to those obtained with conventional Doppler signal processing. In addition to robustly resolving velocity ambiguity, the MAP velocity estimator is shown to reduce high-frequency noise in turbulence spectra.

## 2. Theory

### a. MAP estimation

MAP estimation is one of several parameter estimation methods based on Bayesian inference. In contrast with classical statistical estimation, Bayesian methods treat a parameter to be estimated as a random variable rather than a deterministic but unknown constant (Box and Tiao 1973). The motivation for the Bayesian approach is that prior knowledge of the PDF may be used to improve estimator performance. Once observations have been made, Bayes’ theorem permits calculation of the posterior PDF, which expresses the likelihood of a particular parameter value given the measurements that have been observed. An estimate is then calculated as a function of the posterior PDF, for example, as the location of the maximum value in MAP estimation (Kay 1993).

The unknown parameter to be estimated is the radial speed *υ* (for monostatic sonar) or the velocity vector **v** (for multistatic sonar). Let ** μ** denote a vector of measurements consisting of the autocorrelation phase and autocorrelation magnitude for each carrier frequency and transducer. The MAP velocity estimate is the location of the maximum value of the posterior velocity PDF, denoted by (Kay 1993, chapter 11)

The term *p*(**v**) describes any prior knowledge of **v** that exists before measurements have been observed. The measurement PDF *p*(** μ** |

**v**) expresses the probability of observing measurements

**when velocity takes the value**

*μ***v**.

Because autocorrelation phase can only be measured modulo 2*π* radians, the velocity likelihood function for a single phase measurement is periodic with period 2*υ _{a}* determined by the ambiguity velocity. For example, a measurement with and

*υ*= 0.25 m s

_{a}^{−1}is illustrated in Fig. 1a. Each peak in Fig. 1a is Gaussian with a standard deviation of 0.06 m s

^{−1}. From only a single measurement, it is not possible to determine the true velocity unless prior knowledge can be invoked to restrict the range of possible velocity values.

In a multifrequency system, velocity ambiguity may be resolved by noting that ambiguity velocity is inversely proportional to each carrier frequency (Zedel and Hay 2010). For example, velocity likelihood functions are shown in Fig. 1b for a true velocity of 0.5 m s^{−1} and carrier frequencies with ambiguity velocities of 0.23, 0.25, and 0.27 m s^{−1}. In this example, all three likelihood functions coincide in the vicinity of the true velocity. Although each function is periodic, periods differ according to the distinct ambiguity velocities.

Independent measurements in Fig. 1b may be combined by multiplying the corresponding probability distributions, as shown in Fig. 2a. No measurement errors have been introduced; the existence of multiple peaks in Fig. 2a is solely the result of measurement uncertainty and the presence of velocity ambiguity.

Prior knowledge may be used to further constrain the velocity likelihood function. In the example posterior PDF shown in Fig. 2b, prior knowledge has been represented by enforcing the restriction |*υ*| < 0.75 m s^{−1}. In Fig. 2b, the mean and median values for the posterior PDF are biased by the presence of multiple peaks. These two methods are therefore not suitable when the posterior PDF is multimodal. In contrast, the MAP estimator correctly selects the most likely peak.

In Fig. 3a, velocity likelihood functions are shown for the case where measurement uncertainty has been reduced to 0.02 m s^{−1}. The product PDF in Fig. 3b consists of a single large amplitude peak centered around the true velocity of 0.5 m s^{−1}. Thus, Fig. 3b shows that the multiplication of probability densities automatically rules out infeasible combinations of ambiguity intervals when measurement uncertainty is low. A key feature of the MAP velocity estimator is that when measurement uncertainty is high, for example, in turbulent flow, velocity is estimated as a continuous parameter rather than as a discrete number of ambiguity intervals. Smoothing is then performed to select the most likely value using a simple model based on temporal continuity, as described in section 2d.

### b. Bistatic Doppler shift

A bistatic geometry is depicted in Fig. 4, where the scatterer, transmitter, and receiver positions are indicated by *S*, *T*, and *R*, respectively. Unit vectors directed from the scatterer to the transmitter and receiver are denoted by **r*** _{T}* and

**r**

*. The transmitter–receiver pair measures a velocity component in the direction of the bistatic baseline, that is, toward the line segment RT and along the bisector of the angle ∠RST. Let*

_{R}**r**

*denote a unit vector directed toward the baseline in the direction*

_{B}**r**

*+*

_{T}**r**

*, as shown in Fig. 4. The angle subtended by vectors*

_{R}**r**

*and*

_{T}**r**

*is denoted by 2*

_{R}*θ*.

If a sinusoidal acoustic source of frequency *f*_{0} emanates from the transmitter, then the observed frequency *f _{R}* at the receiver is given by

Because |*υ*| ≪ *c* for oceanographic applications, (8) may be approximated with a first-order Taylor series to obtain

Because **r*** _{T}* +

**r**

*= 2*

_{R}**r**

*cos*

_{B}*θ*, the Doppler frequency shift

*f*=

_{D}*f*−

_{R}*f*

_{0}is

Thus, a bistatic configuration measures the velocity component *υ _{B}* =

**v**·

**r**

*toward the bistatic baseline. When the covariance method is used to estimate Doppler frequency, autocorrelation phase and the velocity component*

_{B}*υ*are related by

_{B}### c. Velocity estimation

A symmetric two-dimensional multistatic sonar geometry is illustrated in Fig. 5. Acoustic pulses are transmitted from transducer 3 and backscatter is received by all three transducers. Velocity measurement is possible in the region where transducer beam patterns overlap, as indicated with dashed lines in Fig. 5, producing a velocity profile that extends above and below the nominal intersection point by approximately 10 cm (Hay et al. 2008).

Let *γ* denote the tilt angle of transducer 1 and 2 baselines with respect to the *x* axis in Fig. 5. For the isosceles geometry in Fig. 5, *θ* = *γ* for a scatterer at the beam intersection point. Velocity components measured by each receiver are denoted by *υ _{j}* (

*j*= 1, 2, 3) with unit vectors

**r**

*directed from the beam intersection point toward the corresponding baselines for*

_{j}*j*= 1, 2 and toward the center transducer for

*j*= 3. If the velocity vector is specified in Cartesian coordinates as [

*υ*], then velocity measurements are given by

_{x}υ_{z}For comparison with MAP estimation, a least squares inverse transformation will be used to determine Cartesian velocity components from transducer measurements,

Division by sin*γ* in (15) indicates that the transverse velocity component *υ _{x}* is more difficult to estimate than the radial component

*υ*. For a tilt angle

_{z}*γ*of 7°, the measurement error ratio

*σ*/

_{x}*σ*is approximately equal to 10.

_{z}In the least squares formulation for *υ _{x}* and

*υ*, each measurement is weighted equally. In contrast, the MAP velocity estimator uses (12)–(14) to project measurement PDFs onto a two-dimensional velocity space, thereby weighting each measurement according to its PDF.

_{z}For a multifrequency system with *N _{f}* carrier frequencies and

*N*transducers, measurements consist of autocorrelation phase and autocorrelation magnitude corresponding to each carrier frequency

_{t}*f*and transducer

_{i}*j*. The MAP velocity estimator takes the form

Numerical simulation of a Gaussian random process was performed to determine the measurement PDFs . The backscatter autocorrelation sequence *R _{k}* for simulation of a Gaussian random process is given by (Dillon et al. 2011a)

where autocorrelation in (1) has been normalized by *R*(0). Let ** ζ** = [

*ζ*

_{1}, … ,

*ζ*]

_{L}^{T}denote an ensemble of

*L*independent identically distributed samples from a complex Gaussian distribution with zero mean and unit variance, where the superscript T denotes the matrix transpose. Simulated backscatter samples

**z**= [

*z*

_{1}, … ,

*z*]

_{L}^{T}were generated using the transformation

**z**=

**U***

**, where**

*ζ***U*** is obtained from Cholesky decomposition of the desired covariance matrix (Kay 1993, chapter 15),

A histogram of phase error was calculated for 10^{6} ensembles of 10 pulse pairs with autocorrelation coefficients *ρ* ranging from 0 to 0.98 in increments of 0.02, and additional values of 0.99, 0.995, and 0.999. To evaluate , the expected phase *φ* was calculated from (3), an unbiased estimate of autocorrelation magnitude *ρ* was calculated from using the bias correction in Dillon et al. (2011a), and the histogram of phase error was interpolated linearly to obtain the corresponding probability density.

### d. Filtering and smoothing

The following presentation extends the one-dimensional MAP velocity estimator described in Dillon et al. (2011b) to higher dimensions. Superscripts *f*, *b*, and *s* denote results from the forward filter, backward filter, and smoother, respectively. Let denote the forward-time posterior PDF, that is, *p*(**v*** _{n}* |

*μ*_{1}, … ,

*μ**), describing the likelihood of*

_{n}**v**

*occurring after the sequence of measurements {*

_{n}

*μ*_{1}, … ,

*μ**}. The sample-to-sample velocity increment is assumed to obey a Gaussian probability distribution. Therefore, given a posterior PDF , the predicted likelihood is the convolution of with a Gaussian distribution (Rosenthal 2000),*

_{n}where *ν* is the dimension of the velocity vector, “det” is the determinant, and denotes the covariance matrix of the velocity increment.

The forward-time MAP velocity estimator is obtained by maximizing the posterior PDF as in (17),

The forward-time prior distribution is initialized to be uniform and is updated using the convolution operation in (20).

Let denote the backward-time posterior PDF, that is, *p*(**v*** _{n}* |

*μ**, … ,*

_{n}

*μ**), describing the likelihood of*

_{N}**v**

*occurring after the measurement sequence {*

_{n}

*μ**,*

_{N}

*μ*

_{N}_{−1}, … ,

*μ**}. Analogous to the forward-time case, the backward-time estimator is given by*

_{n}The backward-time prior distribution is initialized to be uniform and is updated as in (20),

The terms and in (21) and (23) describe prior information gained from measurements {*μ*_{1}, … , *μ*_{n}_{−1}} and {*μ*_{n}_{+1}, … , *μ** _{N}*}, respectively.

MAP estimation can also be implemented by combining all of the measurements in the form of a smoother. Let denote the smoother posterior PDF *p*(**v*** _{n}* |

*μ*_{1}, … ,

*μ**) describing the likelihood of*

_{N}**v**

*occurring after all measurements in the time series have been recorded. The MAP smoother is formed by combining the results from forward and backward filtering,*

_{n}With this formulation, the smoother utilizes all of the available measurements to estimate **v*** _{n}*. The product reflects the independence of two predictions for

**v**

*derived from disjoint measurements sets {*

_{n}

*μ*_{1}, … ,

*μ*

_{n}_{−1}} and {

*μ*

_{n}_{+1}, … ,

*μ**}.*

_{N}## 3. Numerical simulation

To evaluate MAP estimation of a known velocity, simulations of oscillating flow were performed using the coherent Doppler sonar model described in Zedel (2008). The model simulates pulse-to-pulse coherent scattering from a cloud of moving particles for arbitrary multistatic sonar geometries. Physical effects such as spherical spreading, acoustic absorption, frequency-dependent beam patterns, transducer frequency response, and receiver noise are included in the model. The model supports simulation of arbitrary pulse shapes, including the use of multiple carrier frequencies.

Oscillating horizontal flow was simulated with a sinusoidal amplitude of 4.0 m s^{−1} and a period of 10 s. Thirty seconds of simulated measurements were generated, representing three periods of oscillation. In the model, the multistatic configuration shown in Fig. 5 was tilted 5° from vertical to reproduce the geometry of the grid turbulence experiment described in section 4. Parameters for the coherent Doppler sonar simulation are listed in Table 1.

In Fig. 6a, the phase time series is presented for the 2.1-MHz receiver channel of the center transducer in the oscillating flow simulation. Numerous phase discontinuities at ±*π* are present. Similar results were obtained for carrier frequencies at 1.2, 1.5, and 1.8 MHz, except that phase discontinuities occurred at a higher speed for lower carrier frequencies. When velocity is near zero, for example, at *t* = 5 and 10 s, phase measurement noise is minimal. However, when the speed is high, for example, at *t* = 2.5 and 7.5 s, phase measurement noise increases resulting from backscatter decorrelation caused by scatterer advection through the sample volume. The oscillating flow simulation therefore provides a dataset in which both velocity ambiguity and significant measurement noise are present.

The time series of autocorrelation magnitude is presented in Fig. 6b for the 2.1-MHz receiver channel of the center transducer in the oscillating flow simulation. Autocorrelation magnitude also exhibits periodic behavior corresponding to the phase noise in Fig. 6a. When velocity is near zero, autocorrelation magnitude for all frequency channels approaches one. When the speed is high, autocorrelation coefficients as low as 0.1 are observed. However, even during periods of peak velocity, high (i.e., greater than 0.7) autocorrelation magnitudes are still observed. Thus, accurate velocity estimates may be obtained if samples are combined in a way that emphasizes the higher-quality measurements.

### a. Radial velocity

The MAP velocity estimator was configured to calculate velocity likelihood functions on a one-dimensional grid extending over ±1 m s^{−1} in increments of 0.01 m s^{−1}. Gaussian interpolation was used to refine the location of the likelihood peak (Raffel et al. 2007, chapter 5). Velocity smoothing was performed as described in section 2 with in (20) and (25) and *σ* equal to 0.01 m s^{−1}. The smoothing parameter *σ* is a user-configurable parameter that represents the expected amplitude of velocity increments from sample to sample. The value of *σ* was chosen to be less than the radial ambiguity velocity and greater than the maximum sinusoidal velocity increment in one ensemble interval *Mτ*. For the oscillating flow simulation, radial velocity varies as

where *V* is the peak horizontal velocity and *T _{o}* is the period of oscillation. Therefore, the maximum velocity increment is given by

which is equal to 0.003 m s^{−1} for the parameters of the oscillating flow simulation. When choosing the number of pulse pairs *M* per ensemble, there is a trade-off between maximizing the sample rate (small *M*) and reducing the variance of individual measurements (large *M*). Equation (29) implies that the sample rate 1/*Mτ* should be sufficiently fast to ensure that the sample-to-sample velocity increment is small compared to the ambiguity velocity.

In Fig. 7, radial velocity from MAP velocity estimation is compared with the average of the multifrequency measurements from the center transducer. Velocity ambiguity has been removed from the multifrequency curve in Fig. 7a by adding multiples of 2*π* to each phase measurement to minimize the difference between the measured and actual radial velocities. The MAP estimator, however, operates on raw phase measurements constrained to the interval [−*π*, *π*].

The time series of velocity error in Fig. 7b shows that measurement error for the multifrequency measurements increases when the speed is high. In addition to resolving velocity ambiguity, the MAP velocity estimator produces an error time series with a standard deviation of 0.005 compared to 0.011 m s^{−1} for the average of the multifrequency measurements. Also, the multifrequency error time series in Fig. 7b is clearly periodic. Although the quality of raw measurements varies in time, reduced RMS error is possible using the MAP velocity estimator because each measurement is weighted according to its probability density function.

In practice, velocity ambiguity cannot be removed from individual measurements as in Fig. 7a because the actual velocity would not be known. In Fig. 8, radial velocity from the MAP velocity estimator is compared with two other methods for velocity ambiguity resolution. All three algorithms operate on the same set of raw measurements. The multifrequency phase resolution algorithm processes measurements from each sample in time to derive a linear relationship between phase and frequency based on an estimate of the slope *dφ*/*df*, as described in Zedel and Hay (2010). For the temporal continuity algorithm, offsets of 2*π* were added to the phase time series for each frequency channel to minimize phase differences between successive samples. Occasional spikes occur in the multifrequency phase resolution time series because of the difficulty of estimating *dφ*/*df* as Δ*φ*/Δ*f* when the phase noise amplitudes are high. However, spikes could be removed with low-pass filtering or a despiking algorithm.

The temporal continuity method shows a much more serious error occurring at *t* = 7.7 s. A series of poor quality measurements causes an ambiguity resolution error in the *f*_{1} frequency channel. The error persists until *t* = 22.3 s when another error happens to cancel the *f*_{1} error. However, an error occurs in the *f*_{4} channel at *t* = 22.5 s, followed by an error in the *f*_{3} channel at *t* = 27.5 s. The temporal continuity method is therefore not robust for bursts of low-quality measurements, such as those occurring at the peak horizontal velocity of 4.0 m s^{−1}. Ambiguity errors for both the multifrequency phase resolution method and the temporal continuity method disappear when the complex autocorrelation sequence is low-pass filtered with a 10-Hz Butterworth filter before calculating velocity measurements. However, low-pass filtering leads to a reduction in bandwidth, which may be unacceptable for applications in turbulent flow.

### b. Transverse velocity

The MAP velocity estimator was also configured to calculate velocity likelihood functions on a two-dimensional grid extending over ±5 m s^{−1} in the horizontal direction and ±1 m s^{−1} in the vertical direction, with a grid spacing of 0.02 m s^{−1}. The rationale for choosing a smaller range for vertical velocity is that, in practice, one-dimensional MAP velocity estimation may be used with center transducer data to quickly estimate the range of expected radial velocities. Velocity smoothing was performed using and *σ* = 0.02 m s^{−1}.

In Fig. 9, the transverse velocity *V _{X}* (i.e., perpendicular to the center transducer axis) from the MAP velocity estimation is compared with the transformation given by (15) in section 2. As for the radial case, velocity ambiguity has been removed from the multifrequency curve in Fig. 9a by adding multiples of 2

*π*to autocorrelation phase to minimize velocity error.

As before, the MAP estimator operates on raw phase measurements constrained to the interval [−*π*, *π*]. The time series of velocity error in Fig. 9b also shows a periodic increase in measurement error for multifrequency measurements when speed is high. Also, transverse error is greater than radial error from Fig. 7b, because transverse velocity is more difficult to estimate than radial velocity using a multistatic geometry with a small baseline tilt angle, as discussed in section 2. In addition to resolving velocity ambiguity, the MAP velocity estimator produces an error time series in Fig. 9b with a standard deviation of 0.018 compared to 0.063 m s^{−1} for the average of the multifrequency measurements. As in Fig. 7b, reduced RMS error is possible because the MAP velocity estimator weights each measurement according to its probability density function.

## 4. Experimental instrumentation

A grid turbulence experiment was performed in the Marine Craft Model Towing Tank at Dalhousie University. The tank has horizontal dimensions of 30 m × 1 m and a depth of 1 m. An instrumented carriage is propelled by an electric motor along rails mounted above the tank. Carriage speed is computer controlled and programmable over a range of 0–3.0 m s^{−1}. Constant speed is sustained over a rail length of approximately 25 m. The towing carriage, rectangular grid, and instrumentation are shown schematically in Fig. 10.

Velocity was measured using the multifrequency coherent Doppler sonar (MFDop) described in Hay et al. (2008). The sonar employs a symmetric multistatic geometry, as shown in Fig. 5. In Hay et al. (2008), a longer-range transducer head with an angle of 16.6° between transducers and a beam intersection distance of 70 cm was used. In this experiment, a shorter-range transducer head was used, featuring an angle of 14° between transducers and a beam intersection point located 40 cm from the face of transducer 3.

As in Hay et al. (2008), each transducer has a diameter of 2 cm, a nominal center frequency of 1.7 MHz, and a bandwidth of approximately 1 MHz. Carrier frequencies, profiling range, range resolution, pulse length, pulse-to-pulse interval, and ensemble length are configurable in software. Nominally, each sample point has a diameter of 2 cm and a height of 3 mm. The parameters in Table 1 also apply for the MFDop. The sonar was rotated to point 5° aft (i.e., counterclockwise in Fig. 10) to avoid receiving multiple reflections from the tank bottom and water surface. Also, the surface-piercing support for the MFDop transducer assembly was fared to minimize bow wave and wake effects. The sonar was located on the tank center line with the center transducer 56 cm above the bottom.

Measurements were also collected with a 10-MHz Nortek Vectrino acoustic velocimeter to independently measure carriage speed. The center transducer of the Vectrino was located 45 cm above the bottom, 16 cm forward of the MFDop center transducer, and 19 cm to the left side of the tank center line. Vectrino measurements were recorded at 100 Hz using Vectrino^{+} firmware. Accuracy is quoted by the manufacturer to be ±0.5% of the measured value.

A rectangular grid was mounted 79 cm forward of the Vectrino and 95 cm forward of the MFDop, as indicated in Fig. 10. The grid consisted of 11 × 11 aluminum bars with a 1 cm^{2} square cross section. Each of the openings between the grid elements was square with size *D* = 4 cm. Therefore, the Vectrino was located approximately 20*D* downstream of the grid, corresponding to the beginning of the homogeneous initial period of decay for grid turbulence (Batchelor and Townsend 1948).

Water in the tank was seeded with agricultural lime. Prior to each run, approximately 0.5 kg of lime was added to replace scatterers lost to settling. An order of magnitude estimate of sediment concentration was 1 g L^{−1}.

Carriage speed was varied from 0.1 to 1.0 m s^{−1} by programming the desired speed into the control system. Results are presented in section 5 for a nominal speed of 1 m s^{−1}. Because of the additional drag from the grid, actual carriage speed was less than the requested value by approximately 10%, as determined from velocity measurements with the Vectrino. The control system software automatically calculated an acceleration and deceleration profile to maximize the time at constant speed subject to the tank’s length constraint.

## 5. Experimental results

Time series of autocorrelation phase and magnitude are shown in Figs. 11a,b. MFDop results are presented for transducer 1 at a range of 40 cm for the 2.1-MHz carrier frequency. The towing carriage accelerated to a nominal speed of 1 m s^{−1}, as shown by the time series of horizontal velocity measured by the Vectrino in Fig. 11c. A phase discontinuity resulting from velocity ambiguity is evident at the 1-s mark in Fig. 11a. Similar results were obtained for carrier frequencies at 1.2, 1.5, and 1.8 MHz, except that the phase discontinuity occurred later (i.e., at a higher carriage speed) for lower carrier frequencies. The time series of autocorrelation magnitude in Fig. 11b shows a decreasing trend during carriage acceleration. Furthermore, measurement noise increased abruptly as the MFDop entered the turbulent wake of the grid at approximately the 2-s mark.

Three methods for velocity ambiguity resolution are compared in Fig. 12a for measurements from a single MFDop receiver in a bistatic configuration. For a sonar tilt angle of 5°, the velocity component *V*_{1} measured by transducer 1 is inclined 12° to the vertical. Therefore, for a carriage speed of 1 m s^{−1}, the *V*_{1} component is nominally 0.2 m s^{−1}. In Fig. 12b, horizontal and vertical Vectrino measurements were rotated by 12° to obtain the *V*_{1} velocity component. Results from one-dimensional MAP velocity estimation in Fig. 12a were obtained by evaluating velocity likelihood functions over a range of ±1 m s^{−1} with a gridpoint spacing of 0.01 m s^{−1} and a smoothing parameter *σ* of 0.08 m s^{−1}. The multifrequency time series is the result of applying the method in Zedel and Hay (2010) for velocity ambiguity resolution. For the temporal continuity method, multiples of 2*π* radians were added to the autocorrelation phase of each carrier frequency to minimize discontinuities in time before averaging the results from all carrier frequencies.

Velocity spectra for the multifrequency and MAP velocity *V*_{1} time series are shown in Fig. 13 for the steady-state portion of the towing carriage motion. Spikes in the multifrequency time series were removed prior to calculation of the spectrum by shifting phase measurements by multiples of 2*π* to best match the MAP velocity time series. Because there is no high-frequency noise floor in Fig. 13, the automatic tuning algorithm in Dillon et al. (2011b) does not apply. Instead, the smoothing parameter of *σ* = 0.08 m s^{−1} was obtained by decreasing *σ* as much as possible without introducing high-frequency attenuation in the velocity spectrum. Results similar to Fig. 13 were obtained when MAP velocity smoothing was applied to measurements from transducers 2 and 3.

In Fig. 14, the *V _{X}* velocity component (i.e., transverse to the MFDop, as shown in Fig. 5) is estimated from measurements from all three MFDop transducers. Results are shown for the range bin at 40 cm. The center transducer is inclined 5° to the vertical. In Fig. 14b, horizontal and vertical Vectrino measurements were rotated by 5° to obtain the

*V*velocity component. Results from two-dimensional MAP velocity estimation in Fig. 14a were obtained by evaluating the velocity likelihood functions over a horizontal range of 0–2 m s

_{X}^{−1}and a vertical range of ±1 m s

^{−1}with a gridpoint spacing of 0.01 m s

^{−1}and a smoothing parameter

*σ*of 0.08 m s

^{−1}, where . For multifrequency and temporal continuity time series, the respective methods were applied to measurements from each transducer before estimating

*V*using (15).

_{X}Velocity spectra for the multifrequency and MAP velocity *V _{X}* time series are shown in Fig. 15 for the steady-state portion of the towing carriage motion. As for the one-dimensional case, spikes in the multifrequency time series were removed prior to calculation of the spectrum. Also shown in Fig. 15 is a spectral slope of corresponding to the inertial subrange of homogeneous isotropic turbulence (Pope 2000). The Nyquist frequency occurs at 33.3 Hz for a pulse-to-pulse interval of 1.5 ms and ensembles of 10 pulse pairs. The velocity spectrum for the

*V*component was similar to Fig. 13.

_{Z}## 6. Discussion

An illustrative example was presented in section 2 to show how multimodal velocity likelihood functions may occur when velocity ambiguity is present. It was shown that estimators based on the posterior mean and median can be biased when the velocity likelihood function is multimodal. When velocity ambiguity is present, estimators that assume Gaussian statistics (e.g., Kalman filters) are also inappropriate and will lead to biased velocity estimates. In contrast with classical statistical estimation, MAP estimation makes use of prior knowledge to identify the most likely peak. Two types of prior knowledge result if it is known that (i) velocity lies in some predefined range, or (ii) velocity is continuous in time and/or space.

The MAP velocity estimator was developed in section 2. Although the discussion focused on the two-dimensional multistatic geometry shown in Fig. 5, estimator equations were presented in a way that applies for one-, two-, or three-dimensional velocity estimation. The MAP velocity smoother combines results from forward and backward filtering to eliminate estimator lag and to reduce measurement uncertainty. Prior knowledge is derived from a model where the velocity increment from sample to sample is assumed to be discrete-time white Gaussian noise. As the smoothing parameter *σ* decreases, the likelihood function is increasingly shaped by the assumption of temporal continuity. However, the resulting velocity estimate is not constrained to obey the model precisely because the estimate is ultimately determined by available measurements.

In the prediction step, convolution is implemented by filtering the posterior PDF of the previous sample with a finite impulse response (FIR) filter having a Gaussian impulse response. MAP velocity estimation could also be applied in the spatial domain using a model based on the assumption of spatial continuity. The smoothing parameter would represent a spatial diffusion of velocity likelihood, and convolution would be implemented by filtering the posterior PDF of adjacent samples on a spatial grid.

A simulation of multifrequency coherent Doppler sonar measuring an oscillating flow was used to evaluate the MAP velocity estimator via a comparison with the known simulated velocity. Measurement noise was caused by randomly positioned scatterers advecting through the sample volume. Numerous phase wraps at ±*π* were also present because of velocity ambiguity. The one- and two-dimensional MAP velocity estimators were shown to resolve velocity ambiguity while suppressing measurement noise during periods of peak flow speed. The key observation was that although autocorrelation coefficients as low as 0.1 were observed, low-quality measurements tended not to occur simultaneously on all frequency channels and for all transducers. Thus, by combining measurements according to their corresponding probability density functions, it was possible to suppress spurious measurements, resulting in an amplitude error that was approximately constant in time despite the varying quality of raw measurements.

The grid turbulence experiment also provided a dataset where both velocity ambiguity and backscatter decorrelation were present. Time series of autocorrelation phase and magnitude in Fig. 11 indicate that there was an increase in measurement noise during carriage acceleration, with additional measurement noise occurring as the MFDop entered the turbulent wake of the grid. The Vectrino and MFDop velocity time series in Figs. 12 and 14 are not directly comparable because of the spatial separation of the sensors, as indicated in Fig. 10. However, the Vectrino provides an independent measurement of carriage speed, which is useful for evaluating the performance of algorithms for velocity ambiguity resolution.

Figures 12a and 14a show that backscatter decorrelation leads to velocity ambiguity resolution errors for the multifrequency and temporal continuity methods. Because the multifrequency method processes each point of the time series independently, occasional errors appear as velocity spikes. The algorithm has no “memory,” and therefore spikes can be removed with a despiking algorithm. The temporal continuity method, on the other hand, has infinite memory, in the sense that ambiguity resolution errors persist indefinitely unless two errors happen to cancel. MAP velocity smoothing can be interpreted as a compromise between these two extremes where the best features of each approach are retained. The estimator incorporates temporal information to help choose the correct peak in a multimodal velocity likelihood function. However, temporal continuity is not enforced in a strict sense, and the estimator is able to recover if an error occurs. The two-dimensional MAP velocity estimator is more tolerant of measurement noise than the one-dimensional methods because probability distributions from all three transducers are combined before estimating velocity.

While large amplitude spikes may be removed from the velocity time series, lower amplitude measurement noise is difficult to separate from real turbulent fluctuations. Averaging and low-pass filtering are frequently applied to reduce measurement errors; however, there is a corresponding reduction in effective sample rate and bandwidth. Ideally, noise suppression methods would improve velocity estimates while preserving the sample rate and bandwidth of the original data. Figures 12 and 13 show an example where one-dimensional MAP velocity estimation robustly resolved velocity ambiguity without introducing high-frequency attenuation in the velocity spectrum. If a noise floor had been present in the velocity spectrum, the one-dimensional estimator could also be used to suppress noise, as demonstrated in Dillon et al. (2011b). Because the MAP velocity estimator operates in the time domain, a smoothed time series at the full sonar sample rate is produced (e.g., Figs. 12a and 14a) in addition to the smoothed spectra shown in Figs. 13 and 15.

For transverse velocity measurement, the component *V _{X}* is more difficult to estimate than radial velocity using a multistatic geometry with a small baseline tilt angle. Figures 14 and 15 show that two-dimensional MAP velocity estimation suppresses high-frequency noise inherent in transverse velocity measurement while resolving velocity ambiguity and retaining high-frequency fluctuations that are characteristic of turbulent flow. Comparison of Figs. 13 and 15 illustrates how two-dimensional MAP velocity estimation differs from the one-dimensional case. The two-dimensional estimator projects measurements from each transducer onto the space of Cartesian components [

*υ*,

_{x}*υ*] so that each measurement is weighted according to the geometry of the sonar in addition to the corresponding autocorrelation coefficient. Smoothing is then applied to [

_{z}*υ*,

_{x}*υ*] rather than to the time series from each transducer.

_{z}The focus of this paper has been on the effectiveness of the MAP velocity estimator rather than computational or implementation issues. Velocity likelihood functions were evaluated on a discrete set of uniformly spaced grid points to estimate the maximum value of the posterior PDF. Although the grid search is not computationally efficient, there are no convergence issues and the search produces the global maximum when the grid spacing is small compared to the width of peaks in the posterior PDF. Velocity estimates were refined by fitting a Gaussian peak to a 3 × 3 grid of points in the vicinity of the maximum value, as is typically performed for subpixel estimation in particle image velocimetry (Raffel et al. 2007, chapter 5). Using the peak fitting algorithm, velocity results were not sensitive to the choice of gridpoint spacing. Future work will investigate the use of parallel processing to evaluate velocity likelihood functions, and the use of an adaptive mesh and more sophisticated optimization algorithms to improve computational efficiency.

## 7. Conclusions

In this paper, a velocity estimator for coherent Doppler sonar has been described. The estimator employs MAP estimation to optimally combine multifrequency and multitransducer measurements. The estimation framework accommodates monostatic, bistatic, and multistatic sonar configurations. The algorithm operates on a time series of measurements from a single point in space. Spatial smoothing is also possible if the flow can be assumed to be spatially continuous. Each phase measurement is weighted according to a PDF determined by a corresponding autocorrelation coefficient. A user-configurable smoothing parameter controls how measurements are combined in the time domain. The MAP velocity estimator was evaluated using a high-fidelity coherent Doppler sonar simulation of oscillating flow. Results were also presented from a towing tank grid turbulence experiment where both velocity ambiguity and backscatter decorrelation were present. Time series and spectra from MAP velocity estimation were compared to those obtained with conventional Doppler signal processing. In addition to robustly resolving velocity ambiguity, the MAP velocity estimator was shown to reduce high-frequency noise in turbulence spectra.

## Acknowledgments

We thank Richard Cheel for experimental support and Robert Craig for development of data acquisition software. Financial support for J. Dillon was provided by the Link Foundation and the Natural Sciences and Engineering Research Council of Canada.