Wave breaking and wave-forced flows are important to air–sea interactions and to the transport and dispersal of materials at sea. But recent measurements have shown a discrepancy in the Eulerian response to wave groups compared to scientists’ current theoretical understanding of wave–current interactions. Flow structures on scales of centimeters to meters occur underneath breaking waves, and larger-scale flows are driven by wave–current interactions (e.g., Langmuir circulation, alongshore flows). Such details of the vertically resolved flow are just beginning to be modeled, and observational guidance is needed. Here a new instrument is described that is intended to measure waves and currents over a 2D vertical plane underwater, resolving two components of velocity on this plane. Initial observations were made near the Scripps Pier (La Jolla, California), where steep waves and strong currents can be reliably found, yet logistics are not too burdensome. To get the spatial resolution desired using 200-kHz sound, ping-to-ping “coherent processing” would have be used for Doppler estimation; however, near shore the reverberations remain too strong for far too long to get any coherence, unlike the previous experience in deep water. In view of this, using much higher frequencies (>1 MHz) with “incoherent processing” is suggested; the increased attenuation at higher frequencies then would subdue the reverberation problem, but with comparable space–time resolution.
Surface waves are central to air–sea interactions (e.g., Thorpe 1982; Edson and Fairall 1994; Pattison and Belcher 1999; Veron et al. 2008; Fan et al. 2009). For example, breaking waves directly inject momentum, gas, and turbulent energy into the mixed layer of the sea and have similar dramatic effects in the air. Although the implications of the intermittent stress transfer due to breaking have historically been the subject of much discussion, they have only recently been simulated and studied in detail (Sullivan et al. 2004, 2007). These simulations indicate that vortical structures resulting from breaking events influence the motion to a depth several times greater than the wave amplitude, deeper than previously thought.
Waves also transport mass and momentum (Stokes 1847; Longuet-Higgins 1953). Longuet-Higgins and Stewart (1962, 1964, p. 530) identified and evaluated the “excess flow of momentum due to the presence of the waves” and, in analogy to optics, named it the “radiation stress.” Changes in the radiation stress (momentum flux) of the waves are compensated for by changes in the mean field, so that overall momentum is conserved. For example, Longuet-Higgins and Stewart (1960, 1961) described the generation of group-bound forced long waves. Field measurements adequate for comparison were not obtained until decades later (Herbers and Guza 1991; Herbers et al. 1994; Smith 2006a). In this last case, a discrepancy was identified: the measured response at the surface was found to be more than twice what theory suggests, as discussed further below. For more details than given below, see Smith (2006a); for a discussion of the wave groups responsible, see Smith and Brulefert (2010).
Wave–current and wave–topography interactions are important in a variety of scenarios, including wave-induced setdown and setup near shore (Longuet-Higgins and Stewart 1964; Bowen et al. 1968); generation of long-shore currents (Bowen 1969; Longuet-Higgins 1970a,b); the interaction of freely propagating long and short surface waves (Longuet-Higgins 1969; Hasselmann 1971; Garrett and Smith 1976; and many others); and the generation of Langmuir circulation (LC), a prominent form of motion in the wind-driven surface mixed layer (e.g., Langmuir 1938; Craik 1977; Leibovich 1980; Li et al. 1995; Skyllingstad and Denbo 1995; McWilliams et al. 1997; McWilliams and Sullivan 2000; Phillips 2002; Sullivan et al. 2007; Kukulka et al. 2009). A key term in the LC case involves the bending of vortex lines by the vertically varying Stokes drift of the waves, called the “vortex force” by Craik and Leibovich (1976). In its vertically integrated form, this term corresponds to a mean momentum change to compensate for the change in wave momentum as the waves are refracted (Garrett 1976; Leibovich 1980; Smith 1980), which is a part of the radiation stress (Smith 2006b), identified first (but not used) by Hasselmann (1971). Much of the work to date on wave–current interactions still depends on depth-integrated balances, as used in the original radiation stress formulation of the interaction. For depth-resolving simulations, the depth-dependent variations of the Stokes drift, including direction as well as magnitude, are needed to evaluate this vortex force and to properly distribute the remaining radiation stress terms over depth.
As mentioned above, the theory as developed to date may not be complete. Referring to comparisons of the predicted and measured surface current response to wave groups, Smith (2006a, p. 1382) states, “It is found that Eulerian counterflows occur that cancel the estimated Stokes drift variations at the surface completely. This is a stronger surface response than expected from the group-bound forced wave analysis” (mentioned above). Also, similar sized Eulerian counterflows have been observed in laboratory experiments (Kemp and Simons 1982, 1983; Jiang and Street 1991; Nepf 1992; Groeneweg and Klopman 1998; Swan et al. 2001; Monismith et al. 2007). “The changes in the Eulerian current profile from case to case (no waves, down-current waves, up-current waves) are comparable to (and oppose) the Stokes drift profile calculated from the wave parameters as given (in both magnitude and shape)” (Smith 2006a, p. 1382). While the Coriolis force would eventually cause such depth-by-depth cancellation (McWilliams and Restrepo 1999), in these experiments the cancellation is observed immediately, over times small compared to the passing of a group, far too fast for vorticity to diffuse downward from the surface. In contrast, pressure fields measured in intermediate depth water have been found consistent with second-order theory (Herbers and Guza 1991; Herbers et al. 1994). However, the expected response at intermediate depths differs markedly from that in deep water, and it is not clear how pressure measurements made close to the bottom should correspond to currents throughout the water column.
Measured velocities in both the laboratory and field show persistent differences from theory that indicate the need for further work and understanding, especially regarding the vertical structure of the forcing and response. This is particularly vital and timely because some observational guidance is needed to aid development of models that resolve the vertical structure of the flows, which is especially important near shore, where wave–current interactions can be the dominant driving force (e.g., Groeneweg and Battjes 2003; Mellor 2003; Özkan-Haller and Li 2003; McWilliams et al. 2004; Smith 2006b; Newberger and Allen 2007; Ardhuin et al. 2008; Kumar et al. 2012). The objective here is to obtain measurements of at least two velocity components over a vertical plane in order to examine the time/space characteristics of both the wave motions and the resulting response as wave groups pass and, in particular, to explore the vertical structure of the response.
Another great challenge in the study of the air–sea interface is to measure and understand natural breaking waves and the resulting subsurface flow structures. It is exceedingly hard to get quantitative measurements closer than a few wave amplitudes to the surface with actively breaking waves, although many promising techniques are being developed (e.g., Vagle and Farmer 1992; Lamarre and Melville 1994b; Farmer et al. 1998; Veron and Melville 1999; Terrill and Melville 2000; Melville et al. 2002; Wu and Nepf 2002; Gemmrich and Farmer 2004; Gemmrich et al. 2008). The objective here is to characterize the subsurface flow structures under naturally occurring breaking waves. Emphasis is not on the small-scale turbulent kinetic energy (TKE) or dissipation (cf. Agrawal et al. 1992; Terray et al. 1996; Gemmrich and Farmer 2004), but on the slightly larger-scale vortical structures (cf. Melville et al. 2002, p. 221). As noted there, “The dominant feature of the post-breaking velocity field is the coherent eddy that drifts slowly downstream at the same speed as the patch of smaller-scale turbulence.” It has also been shown that breaker shape has an effect on the flow and resulting vortex (cf. Nepf et al. 1998). These “coherent structures” should be readily measurable by the proposed technique, although we cannot expect the same level of detail as is possible in the laboratory. Observations to date can help broadly characterize the scales and frequency of breaking waves (Gemmrich and Farmer 1999; Melville and Matusov 2002), and laboratory experiments can reproduce the flow structures under a few select breaker types (cf. Melville et al. 2002). However, there remain two other problems requiring observations from nature for guidance: 1) natural wave groups generally persist longer than the “point break” subjects of laboratory studies, and in addition to breaking repeatedly, they should be accompanied by the group-forced response referred to above (Smith 2006a,b); and 2) waves in nature seldom align perfectly with the background vertical shear, as they do in the laboratory. Thus, the wave groups, breaking, and mean flow, including the response, should be studied together, both to facilitate separating their effects and to include any joint effects resulting from their (nonlinear) superposition, for example, in the mean-flow response.
Here we describe our recent efforts to develop the technology to make continuous multicomponent velocity estimates over a 2D “vertical slice” area.
The approach evolved from a prototype deployment of an up-looking coherent phased-array system done in 1999. That prototype provided roughly 16 independent beams covering a 22°-wide pie-shaped wedge, extending upward from a deployed depth of 13.5 m. The wedge was sampled up to 52 times per second, a rate suitable for examination of the rapid evolution of the flow under breaking conditions (Fig. 1). The surface was located in each frame using a multibeam technique developed to minimize the effect of dropouts caused by dense bubble clouds. The prototype data analysis showed that image correlation techniques are effective for estimating the subsurface vertical velocities associated with the bubble clouds as well as for locating the surface (Smith 2001). It was also found that, under the actively breaking crest, a phase shift occurs in the vertical motion with depth, such that the deeper motion (1–3 m below the surface) is delayed by up to a second or so relative to the surface itself (Smith 2005). With multiple beams, the shielding of the surface by dense bubble clouds is rarely complete. The bubble clouds are sufficiently nonuniform that some fraction of the beams can almost always “sneak through” to the surface and back. Because attenuation is also correlated with the bubble-cloud density, one may worry about a corresponding reduction in sound speed, which would affect estimates of both the range and velocity; however, for acoustic frequencies greater than 40 kHz or so, the sound speed changes due to bubble density are negligible (Lamarre and Melville 1994a; Terrill and Melville 1997). Thus, from the portions observed between the densest bubble clouds, it should be feasible to reconstruct the larger-scale motions such as the “primary vortex” resulting from breaking.
In addition to preliminary results on wave breaking (Smith 2001, 2005), the prototype deployment has provided many ideas for improvements, in terms of array and deployment geometry, most of which were implemented in the newer instrument described here. Chief among these are 1) place the instrument closer to the surface to reduce the lag between pings and facilitate coherent Doppler processing; 2) increase the number of channels to improve horizontal resolution and permit better digital beamforming; 3) expand to bistatic or multistatic transmit/receive geometry to get more than one component of velocity; and 4) extend the beamformed angle out to a full 180° aperture to reduce spatial aliasing.
The basic concept for acoustic measurement of velocities is straightforward, and the performance characteristics for the simple case of a monostatic narrow-beam Doppler sonar are fairly well known for both coherent (e.g., Lohrmann et al. 1990 and references therein) and incoherent (e.g., Pinkel and Smith 1992 and references therein) operation. Extension to a 2D area of measurements is also straightforward—for example, using a linear array of receivers with digital beamforming (cf. Smith 2002). As demonstrated by Smith (2002), performance along any one of many digitally formed beams is in line with that of a single simple beam. Bistatic configurations have also been used, for example, with a narrow-beam transmitter and several fan-beam receivers offset in various directions to obtain multiple velocity components along the line defined by the transmitted beam (e.g., Rolland and Lemmin 1997; Hurther and Lemmin 1998, 2008; Hay et al. 2012a,b, and references therein), and have also been used extensively in atmospheric radar applications.
3. Development of the bistatic phased-array Doppler sonar system (BiPADS)
The basic approach was to build on prior proven designs, extending the capabilities: 64 receive channels; purely resistive tuning (relaxing the restriction to a narrow frequency range); and two or more transmitters located at a distance in the view plane to provide multiple components of velocity estimated from Doppler shifts, using a bistatic configuration (see Fig. 2).
Use of a bistatic configuration with a digital multibeam receiver (noncollocated transmitters and an array of receivers) introduces two complexities to processing: 1) the location from which backscatter is received is no longer a simple linear function of the time delay since transmission; and 2) the Doppler shift corresponds to the velocity component along the direction that bisects the angle between transmitter, target location, and receiver (and hence is a function of both receive angle and range or time delay). While the geometry guiding this is straightforward (and is illustrated below, after defining the proposed array geometry), the high data rate contemplated here means that efficient algorithms are needed. In the following subsections, we review 1) the deployment geometry, 2) acoustic properties and limitations, and 3) bandwidth and processing requirements to obtain estimation errors within tolerable bounds.
There are several technological obstacles to be dealt with (aside from the electromechanical complexity of the system itself), including the effects of strong advection on the two acoustic paths and variations in the speed of sound in the presence of bubbles. For sound at frequencies well over 40 kHz, the latter is a weak effect, with variations limited to under 0.33% or so at 40 kHz and decreasing above that (Lamarre and Melville 1994a). For round-trip distances of up to 15 m, this could account for misidentification of the scatterer location by up to 5 cm, about the resolution hoped for; however, it would require the bubble cloud to fill the entire water column over the whole acoustic path. For bubble clouds on the order a few meters in size, the error in location is proportionally less and can likely be ignored (conversely, the sound speed anomaly field likely cannot be estimated from the data with any useful accuracy). Similarly, a 0.33% error in the velocity estimate would likely not be detectable.
a. Bistatic geometry
Interpretation of the backscatter data can be considered a two-step process: first, the time-angle input data has to be interpolated to a regular grid; second, the direction and errors of each velocity estimate have to be assessed.
We start by defining a grid of “output locations” G = (gx, gy); for example, a regular rectangular grid with 5-cm spacing, relative to the center of the receive array, defined to be at R = (0, 0) (see Fig. 3). The transmitter locations are defined as T = (tx, ty) on one side and S = (sx, sy) on the other. In the following, the method is described for a particular grid point G and transmitter T; extension to other grid points and to the other transmitter is straightforward.
The angle from the center of the receive array to G is
(where the two-argument version of arctan is used to eliminate the 180° ambiguity of the one-argument version). The time delay from transmitting to receiving the signal from G is just the time it takes for sound to travel the distance from T to G plus the distance from G to R: Let
so the time delay for the signal from G is written as
where c is the speed of sound. Standard processing of the receiver array data yields a 2D array of backscatter returns versus angle and delay time, so the above equation enables us to locate the returns coming from point G (etc.).
While much can be learned from just acoustic intensity maps, we would also like to estimate the Doppler shift and to know to what directional component of scatterer velocity this corresponds. The Doppler shift arises from the rate of change of the total pathlength. First, define the (radian) acoustic wavenumbers along each path: kT (along LT from T to G) and kR (along LR from G to R). Then, the Doppler shift from transmitter T, T, can be written as
(note the minus sign on kR due to its definition as going from G to R). We evaluate the vector k (kT and kR) from the center frequency , the speed of sound c, and the known path directions:
where is as given above, and (similarly)
Then we find
where we have written V(gx, gy) = (Vx, Vy) and used trigonometric angle-sum and angle-difference identities. As anticipated at the start of the section, the direction of the Doppler estimate component bisects the angle TGR (see Fig. 3). Also, the size of the Doppler shift is reduced relative to the monostatic case, approaching zero response as the angles approach opposing directions. Graphically, in Fig. 3, the Doppler shift estimates correspond to velocity components perpendicular to the ellipses of constant time delay, and where the two sets of ellipses are more nearly perpendicular to each other there is more nearly equal resolution of both components of velocity.
Naming the corresponding velocity component “VT,” it would be estimated as
Note that the paths, angles, Doppler estimates, and sound vector ks are all functions of the gridpoint location G as well as which transmit is active. This in effect yields a set of transforms to go from the beamformed angle-delay matrix for the backscatter from each transmitter to output arrays of intensity and Doppler shift (or VT and VS) versus location on a uniform Cartesian grid.
b. Combining the transmits
The two offset transmitters can provide Doppler estimates of velocity components in two different directions, generally not orthogonal but hopefully not parallel either. Since the Doppler estimates inherently contain noise (as discussed below), simple Gaussian elimination is unstable, and an estimate explicitly accounting for the noise is preferred. Here we employ the “optimal estimate” derived by Smith (2002). We define the direction of the two estimated components at each grid point, as evaluated above for the transmitter T at grid point G: let
be the direction of the estimated velocity component VT (as above), and for the “second transmit” component VS, let
Now denote the RMS Doppler velocity errors (to be discussed below) associated with these components as and , and normalize them by the (assumed isotropic) RMS true velocity variances: let
Then we write estimates for the Cartesian components of velocity in the form
Smith’s (2002) solution corresponds to
(note that the denominators are all the same). In the (unrealistic) limit of no noise, this reduces to Gaussian elimination. In the limit of large noise on one of the transmitters (S, thus ), that estimate is ignored, and the other is just reduced by its own error, consistent with minimizing the rms error of the result, for example,
4. Acoustic noise and Doppler error
There are two kinds of noise to deal with: acoustic ambient noise and Doppler estimation error. Because there is a finite amount of ambient noise in the ocean and the backscattered signal attenuates with distance through the water, the signal will eventually fade beyond detection and the estimate becomes purely noise. However, at close range, the acoustic signal-to-noise ratio is quite large, and the ambient noise can be neglected. Even so, there remains a finite amount of Doppler error.
For ping-to-ping coherent processing, returns from the same time delay for consecutive ping pairs are correlated to estimate the Doppler shift (Lhermitte and Serafin 1984; Zedel 2008). In this case, the Doppler noise depends on the inherent decorrelation time of the scatterers, the speed of advection of scatterers into and out of the sample volume, and on the reverberation characteristics of the earlier transmissions over the time period of the subsequent ones. A potential advantage of sampling a 2D planar field is that the correlations (or covariances) can be estimated from averages on a moving trajectory on that plane, reducing the advective decorrelation effect. Presuming the carrier (center) frequency is previously remixed to zero, there should be no problem with phase stability on a moving trajectory.
In deep water, as in the Research Platform (R/P) Floating Instrument Platform (FLIP) trials mentioned above, the pings decayed rapidly enough that the time between pings was mainly limited by the desired range: with 50 pings per second, ranges up to 15 m were sampled, but the correlations remaining at the 200-kHz frequencies employed were too small and erratic to encourage Doppler estimation. In shallow water, in contrast, we found that reverberation persists for too long compared to the decorrelation time scale to get useful data.
To estimate the net noise variance σ2 relative to the signal strength, we make use of the coherence:
where Coh is the ratio of the time-lagged covariance 〈Cτ〉 to the received variance 〈s2〉, based on an average, denoted by the angle brackets (〈 〉), over M sample pairs:
where the asterisk (*) denotes the complex conjugate (see appendix A).
Given the relative noise level σ2, the error associated with the covariances upon which the Doppler velocities are based can be estimated. This turns out to depend also on how the M pairs to average are arranged. One way is to average over successive samples (pings) in time:
where the time between samples (pings) is τ. In contrast, the other way involves independent sample pairs (e.g., in a coherent processing scheme, range averaging the products from a single ping pair, assuming each range bin is independent), so with r as the range index and τ as the time between pings (again):
When the true signal covariance is all real (say = 1), the Doppler error shows up in the imaginary part of the variance. For independent sample averages, the covariance error is split equally between the real and imaginary parts:
where the “1” (real) corresponds to the true covariance, and the term on the right is the net Doppler error variance (see appendix B). For sequential sample averages, the total noise variance is the same (i.e., double that term), but the repeated entries in successive sample pairs have the effect of shifting a portion of the error from the imaginary to the real part:
(see appendix B; Fig. 4). This means the Doppler error is less for sequential sample averaging than for independent averages with the same number of sample pairs, unless the real-part variance is large enough to cause many erroneous negative real-part values. For example, assuming normally distributed noise, if , then more than 1% of the real-part error values would exceed the signal strength (=1), and half of these would be on the negative side. Note that the above equation can be rewritten with the true covariance represented as , but both the math and discussion would be more cumbersome.
5. Prototype bistatic deployment
Trials with a prototype BiPADS were done near the end of the Scripps Pier (La Jolla, California; see Fig. 5) to assess the potential for estimating two in-plane components of velocity, using transmitters offset in the same plane as the receive array (see Fig. 3). The prototype BiPADS was deployed off the side of the pier, on a frame with long arms that unfold to hold the transmitters at fixed locations ~4 m to either side of the receive array. The base configuration is 200-kHz center frequency. Beamforming (to 63 beams, dropping the outermost “Nyquist angle”) was done via time delay; the angular resolution at the center angle is about 1.3°, increasing toward the edges. In these trials, we tested the instrument’s capabilities, first trying a variety of sample and ping rates, and then a variety of center frequencies. Initially. we ran it at the 12.5-kHz sample rate, but we found we could reliably increase this to a maximum ~14.7-kHz rate, limited by the bandwidth of the total throughput from all 64 channels to the host computer. Next we tested the ping interval, starting with a conservative 12.5 pings per second and noting that it took 40 ms before the returns decayed substantially. So, we increased the repeat rate to 25 pings per second (1/40 ms) to see if we could get any detectable coherence. While we got fairly satisfactory imaging of the intensity returns (see Fig. 6), 40 ms is too long to obtain usable coherence between pings: they were uniformly statistically indistinguishable from zero.
For the sample frame illustrated in Fig. 6, the pings alternated between transmitters (at 40-ms intervals), and a pair of pings was averaged, including transmits from both sides. This was mainly to establish whether we could alternate between transmitters and image the same locations (we did not even attempt ping pairs or sequences from each side). The surface is about 6 m up on average, the waves (swell) were ~1-m peak to trough, and the wind was dead calm. The surface is dimmer than usual in such images, both because of the calm conditions and because we put sound-absorbing panels above and below the transmitter bars, limiting the aperture to ~30°, reducing specular hits on the surface above the receivers as well as bottom reflections. There are discrete intensity features at the surface and in the water that move with the waves; the criterion for “satisfactory imaging” here is that we could successfully navigate most of these bright spots from the two transmitters to coincide (there will always be some “outliers” due to directional aliasing, slightly different sample volumes, and anisotropic scatterers). It is reassuring that these spots are mostly in the same place no matter which side the sound comes from; that is, that our mapping from the datastream to physical coordinates is successful: proper location of these intensity features is sensitive to the detailed geometry of deployment, and to the precision of timing. Unfortunately, while the mapping from raw data to a Cartesian grid was successful, the reverberations remained too strong for too long to be able to get any ping-to-ping coherence for Doppler estimation.
More interesting in some ways were the results of our tests of the system’s ability to run off resonance (the transducers are all designed for 200 kHz). We successfully ran at 170, 180, 190, 200, 225, 250, 275, 300, 325, 350, 400, 450, 500, 550, and even 600 (kHz). Some aspects of note: at 170 kHz, the receive array is “endfire” capable (i.e., can beamform over 180°). At 400 kHz, there was interference from an unidentified source; the ability to shift operating frequency is critical for avoiding such nuisances. At all frequencies, the signal was still strong at the maximum distances (~10 m) for the same power setting, with little change in intensity. As the frequency increases, attenuation should increase, but the decorrelation time decreases at about the same rate through this range.
To summarize the results from these deployments, we achieved most of the improvements suggested by prior work but failed at the final step:
Move closer to surface. Check: 6 m.
Increase the number of elements for better beamforming. Check: 64.
Expand to bistatic geometry. Check.
Extend beamforming to 180°. Check, at 170 kHz.
Use purely resistive tuning to broaden range of usable frequencies. Check.
Ping fast enough for coherence. Fail (at least in shallow water).
The logical approach is to return to reliable and robust incoherent Doppler processing at a much higher frequency to get the resolution desired (say, 10-cm spatial resolution and 10 cm s−1 velocity resolution). How high a frequency and what bandwidth are required to obtain this? Using repeat-sequence coding (Pinkel and Smith 1992), the single-ping range–velocity uncertainty product is estimated as (Pinkel and Smith 1992)
where c is the speed of sound (~1500 m s−1), f is the center frequency (Hz), and L is the number of “bits” (samples) in the code. The maximum number of bits, or “code length,” that can be used depends on the frequency, bandwidth (Δf), and a “wrap velocity” (Vw) chosen to exceed 99% of actual velocity values anticipated (we typically choose Vw ≈ 2 m s−1). Setting the covariance phase anomaly to π at a time lag of one code length yields the requirement
So, for example, if we assume we can achieve a bandwidth of 10% of the carrier frequency (and with the values quoted above for c and Vw), then we find L ≤ 18.75, independent of frequency. The code length must be an integer, and it must perform well. There is a good “Barker-like” code of length 19; the next shorter good performer is 13.
With the above values of c, Vw, and bandwidth ratio, we get . At 2 MHz, this yields an rms error of 0.012 66 m2 s−1, corresponding to roughly (11 cm)(11 cm s−1), close to the target resolutions. With a 200-kHz sample rate, the 13-bit code would correspond to just under 5 cm in range, so the desired range resolution would correspond to two code lengths; in contrast, the 19-bit code would correspond to just over 7 cm. However, both codes could accommodate ~14-cm range bins (with a correspondingly smaller velocity error): 3 × 13 bits or 2 × 19 bits are close.
The other limitation is the maximum attainable range: as frequency goes up, so does attenuation. The beam geometry considered here is close enough to that used in Smith (1989) that we can examine the range limit at maximum usable power (limited by nonlinear saturation) as estimated there. Unfortunately, those calculations stop just shy of 1 MHz; but above 200 kHz or so, the maximum attainable range goes very nearly like 1/f, so we can extrapolate: this yields the (surely optimistic) estimates of 80 m at 1 MHz, or 40 m at 2 MHz. In practice, there is also the matter of less efficient scatterers (no bubbles are small enough to resonate at megahertz frequencies). However, experience with commercial sonars operated near 1 MHz in incoherent mode indicates ranges of 11 m are achievable, which is in line with our desired working range (on the order of 10 m). Limiting the range to 10–15 m also means that there is time for 50 pings per second or more, enabling significant additional noise reduction by time averaging; thus, it may also be feasible to operate down to 1 MHz, if need be, to get more range. However, it looks likely that 2 MHz will be workable, maybe even optimal.
With incoherent processing, there is a third source of noise, due to “self-clutter” from the nonoverlapping portion of each sample pair. The analysis of section 4 still applies to the other sources of noise (i.e., to both ambient noise and anomalous scatterers); however, the roles of range averaging versus ping averaging are reversed: for incoherent processing, range averaging is “sequential,” while each ping is independent. This could be helpful in extending the usable range, where the ambient noise is the limiting factor: for example, there can be progressively more range averaging at farther ranges.
The early prototype field work (1999) and initial analysis were supported by ONR (Grant N00014-02-1-0855), with additional analysis and modeling supported by NSF (Grant OCE 06-23679). The new instrument development was supported by NSF (Grant OCE-06-23679). Thanks are due to the OPG engineering team of M. Goldin, M. Bui, Tyler Hughens, and A. Aja for their tireless efforts to design, construct, deploy, and operate the BiPADS, and to help develop the software to acquire and analyze the data. Thanks are also due to many others for the many interesting and educational conversations.
The Expected Value of Coh
Let two independent series be denoted and , n = 1 to M, where 1 represents the signal, and en and fn represent complex uncorrelated Gaussian noise with variance , so
(i.e., noise is uncorrelated from one sample to another). Then we can write the estimate of the lag Δt covariance C1 in the independent sample case as
which has the expected value 1, the unit amplitude signal. In forming the coherence Coh, we are interested in the interaction with the variance estimate , where and are shorthand for the two mean-squared time series used in forming C1:
and the square root is approximately
where we have taken the liberty of setting and neglect terms of order M−2 or smaller. Now we examine the interaction between the noise terms when we form the ratio Coh−1:
(again, neglecting terms of order M−2 and smaller). In contrast, for sequential sample averaging, the off-diagonal elements (via and ) change the result to
[in this case, several additional terms need to be retained; e.g., when inverting C1 for the equivalent of (A7), we carry the analysis to ]. In practice, one can hopefully average enough to neglect the O(M−1) terms.
Real and Imaginary Parts of the Covariance
To examine the real and imaginary parts separately, we change the denotation of the received samples to , where , and and are uncorrelated Gaussian random real variables with equal variance. Then
and, defining as the total noise variance, we can write
Now consider a covariance estimate formed from M + 1 sequential samples: let
which has the desired expected value, . In ping-coherent processing, for example, we can think of the samples as being at any fixed range, with n indexing the ping number in the sequence of pings used for the average.
Next consider the variance of the real part, by taking the expected value of
The first term always contributes with a value of 1 (signal), added M2 times to cancel the outer factor of 1/M2, and so yielding . But note the M cases where n = m yields
[using (B2)]. Here the signal portion accounts for the “1,” which was already counted in the “first term sum” noted above, so this leaves as noise variance the portion .
In this sequential sample average, there are also (M − 1) nonzero terms where n = m + 1, yielding
Again, the signal accounts for the “1,” which was already counted. The same amount results from the (M − 1) cases where n + 1 = m, so the net increase in the real part of the noise variance from these together is double the last term. The net expected value is
Somewhat surprisingly, sequential sample averaging actually reduces the noise variance of the imaginary part: consider the expected value of
all of which are zero unless the qs match. There are again M cases where n = m, yielding (from the first two terms) . For this sequential sampling case, there are again 2(M − 1) cases where either n = m + 1 or n + 1 = m, yielding (from the last two terms) the exact negative of the result from these 2(M − 1) cases for the real part; thus,
Since the “off diagonal” terms of the real and imaginary parts are equal and opposite, the total noise variance of the covariance estimate is the same as found above; the sequential sampling average shifts a portion of the noise from the imaginary part to the real. Heuristically, if there is an outlier sample in the midst of the sequence, then the damage it does as the second element in one sample pair is largely undone when it is the first element in the next sample pair.
In contrast, consider an average over independent sample pairs: rewrite (B3) as
where, for example, the primes denote a second ping and now the index n runs over range bins (considered independent). In this case, the M terms where n = m still contribute, but there is nothing special about n = m + 1 or n + 1 = m anymore, so they do not. Thus, the remaining noise is equally distributed between the real and imaginary parts.
Finally, this analysis also applies to “ping incoherent” processing, but with the roles of ping averaging versus range averaging reversed: we need simply to reinterpret the sequential sample estimate (B3) as samples from an individual ping, with n indexing range bins, and to reinterpret (B10) as primes denoting a second range bin, with n indexing pings to be averaged.