Mapping coastal circulations using moving vehicle acoustic tomography

Advances in technology have resulted in the extensive use of unmanned vehicles for various ocean survey tasks, such as monitoring ocean temperatures and currents over a large area. For example, a seaglider is typically used to collect oceanographic data along a survey dive. Among various measurement techniques, ocean acoustic tomography (OAT) is probably the most efficient method to obtain a comprehensive view of temperature and current properties in the interior ocean. The tomographic technique has been applied to shallow waters and coastal seas. The applications referred to as coastal acoustic tomography (CAT) focus on mapping depth-averaged horizontal (two-dimensional) current fields of the tidal current in coastal seas. With at least a pair of acoustic tomographic sensors installed on the moving vehicles, the travel times measured between the sensors can be used to estimate the ocean temperature and currents. With more tomographic sensors, the number of data grows quadratically with the number of sensors rather than approximately linearly for point measurements.

The basic principle of using the difference between the reciprocal travel times to estimate the currents is derived considering the relative motion between the transceivers.¹ Two transceivers (stations i and j) moving with the corresponding velocities (⁠$ui$ and $uj$⁠) transmit acoustic pulses reciprocally; the travel time for acoustic pulses from station i to station j is denoted by t_ij. The measured differential travel time (DTT) (⁠$tij−tji$⁠) is affected not only by the integral of the current field v in the direction of the ray path from station i to station j (denoted by the unit vector $rij$⁠) but also by the sum of the two stations' velocities along that direction,² 

where c is the sound speed and R is the separation of the stations.

For moving vehicle tomography, the multi-Doppler matched filter processing method² was employed to reveal the acoustic arrival pattern. The response peak of the multi-Doppler matched filter is used to obtain a representative Doppler frequency shift for obtaining the acoustic arrival pattern. This estimated Doppler frequency shift can be used, in part, to correct for the effect of stations' velocities on the DTT. Note that the Doppler frequency shift $Δf$ is related to the relative velocity between two paired mobile stations, as³ 

for the acoustic signal with a carrier frequency f_c transmitted from station i received at station j. Therefore, when both stations are moving, the estimate of the velocity sum of the stations, $(ui+uj)·rij$⁠, requires independent measurements, such as a Doppler velocity log (DVL) for submerged vehicles or a GPS receiver for surface vehicles.

After obtaining the velocity sum for a pair of stations, the corresponding induced time shift can be corrected via Eq. (1), i.e., by subtracting this time shift [the last term in Eq. (1)] from the measured DTT. The corrected DTTs, denoted by d, are used to estimate current fields via the following two inversion methods.

Assuming the ocean is stationary during a selected time period of the experiment, all the corrected DTT data collected during that period can be inverted to obtain the current field. If the ocean evolves quickly during the selected time period resulting in a comparatively sparse set of DTT data, a Kalman filter method⁴ or equivalent could be used to handle the non-stationarity during the experiment. Following the tomographic inversion method described in Ref. 7, the horizontal current field $v(x,y)$ is derived by the stream function $ψ(x,y)$ expanded by the truncated Fourier series at high enough wavenumbers. On a given numerical grid, the discretized $ψ(x,y)$ is represented by the vector $ψ=Fm$⁠. The columns of F are the Fourier sine or cosine functions of different wavenumbers, while the parameter vector m contains the corresponding Fourier coefficients. Thus, the x- and y-components of the gridded current field can be expressed as

where the subscripts x and y of F denote the corresponding partial derivatives. Substituting Eq. (3) into the first term on the right-hand side of Eq. (1), the corrected DTTs for all the rays are related to a discrete number of the Fourier coefficients, $d=Gm$⁠, in the absence of noise.

A general least-squares (LS) solution is obtained by minimizing the objective function including both the residual and solution variances

where the data weight matrix $R−1$ is determined from the uncertainty in d. The model weight matrix $βS−1$ has a positive constant β for separate control and is determined by the L-curve criterion. The general LS estimation of the Fourier coefficients is given by

where $G†=(GTR−1G+βS−1)−1GTR−1$⁠. After absorbing β into the model weight matrix, $S′=S/β$⁠, the total uncertainty of the estimated Fourier coefficients is

Estimation of areal current vectors is obtained using the distributed sensing (DS) method,^2,7 which has the potential for real-time applications. To determine an areal current vector by the on-board processor of the station i using the reciprocal transmissions with the other two stations j and k, one needs the path-averaged currents along those two adjacent paths, v_ij and v_ik. Let the angle between these two paths be θ to investigate how θ affects the uncertainty in the estimated current vector. The current along the direction of bisected angle, defined as the longitudinal current, is $vL=0.5(vik+vij)/ cos (0.5 θ)$ and the transverse current $vT=0.5(vik−vij)/ sin (0.5 θ)$ are obtained from these two measured path-averaged current velocities. The uncertainty of path-averaged current velocity, σ_v, is due to the uncertainties in estimating travel time t and Doppler shift $Δf$⁠. Then, assuming no uncertainty associated with θ, the uncertainties in estimating the longitudinal and transverse currents are, respectively, as follows:

A moving-vehicle experiment was conducted using two moving vehicles (AUV and ship) and one moored station (buoy) in WangHaiXiang Bay nearby Keelung City, Taiwan, on June 27, 2017. The water depth in the study area [indicated by the background contour thin line in Fig. 1(a)] varies from 10 to 65 m. The buoy was moored at a depth of 30 m. The AUV cruised near the shore (following a zigzag line for conducting the localization experiment⁵), while the ship surveyed near the mouth of the Bay along a square trajectory with a side length of about 1 km, for providing various ray angles. Each station was equipped with a tomographic instrument, transmitting m-sequence phase-modulated signals with a carrier frequency of about 18 kHz. The m-sequence consists of 1023 digits with four cycles per digit, resulting in a signal bandwidth of about 4.5 kHz. Hence, the digit length (time resolution) is 0.22 ms, and the signal duration is 0.23 s. The acoustic signal was emitted every 20 s and the experiment duration was 40 min, yielding a total of 121 transmitted signals. The GPS receiver on the tomographic instrument recorded the time and position of the vehicle with a sampling period of 1 s. In addition to the tomographic instruments, both moving vehicles were equipped with DVLs with the Acoustic Doppler Current Profiler (ADCP) mode enabled.

The experiment was conducted during the ebb tide [the time window indicated by the shaded area in Fig. 1(b)]. During the 40-min period, the depth-averaged current velocity measured by the ADCP on the ship [white arrows in Fig. 1(a)] shows the flow moving in an east-southeast direction. The higher current velocity with a magnitude of 1.34 ± 0.28 m/s was observed in the deeper area than the weaker magnitude of 0.17 ± 0.17 m/s in the shallower area (near the shore).

Obtaining the path-averaged currents between the moving vehicles requires not only the DTTs but also the vehicle velocities. First, the relative velocity between each pair of stations was obtained using the representative Doppler shift estimated from the response peak of the multi-Doppler matched filter of the acoustic signals. Figures 2(a)–2(c) show that the Doppler-estimated relative velocities between all three paired stations (circles) agree with the relative velocities derived from the GPS estimates (crosses). The discontinuities in the estimated relative velocities are due to the changing course of the moving vehicles. Better agreements are observed in the estimated relative velocities between the mobile and moored stations [Figs. 2(a) and 2(b)], while more fluctuations are observed between the two mobile stations due to low SNRs [Fig. 2(c)]. Second, for the case of the two mobile stations, the velocity sum is determined using the GPS measurements to obtain the ground velocity for one of the mobile vehicles. Figure 2(d) shows the projected ground velocities of the AUV and the ship on the direction from the AUV to the ship. To obtain the velocity sum, the projected ground velocity of the ship is used to remove its contribution from the relative velocity, Eq. (2). Finally, the time shift induced by the moving vehicles is removed from the DTT.

Validation of the DTT-estimated currents using the DS method is carried out by comparing with the point measurements using the ADCP on the ship [arrows in Fig. 1(a)]. Near the ship, the DTT-estimated current vector is synthesized using the path-averaged current velocities along the paths of ship-buoy and ship-AUV. Then, the ADCP-measured current is obtained by averaging the data over the depth with its corresponding uncertainty indicating the variability of the currents in depth. Since the angle between these two paths, θ, varies between 12.0° and 61.4° [see Fig. 3(c)], the uncertainty (measured in one STD) in the longitudinal current component [error bars with solid circles in Fig. 3(a)] is always smaller than that in the transverse component (error bars with open circles). Especially when $θ<20°$ (the period from 15:52 to 16:06), a relatively large uncertainty is observed in the transverse current component.⁷ Figs. 3(b-1) and 3(b-2) show the currents expressed in the eastward and northward components, respectively. The DTT estimates (solid circles with error bars) are in general consistent with the ADCP measurements (open circles with error bars), except when $θ<20°$⁠.

The spatial distributions of the current field within a triangle are obtained using the DS method for each transmission. When one of the adjacent path lengths, R, is less than 500 m or the angle between these two paths $θ<20°$⁠, the derived current field is inaccurate and imprecise, e.g., at the indices #2 and #37 in Fig. 3(d). As the ship moved to the deeper-water area and $θ>20°$⁠, an anticyclonic circulation within the Bay is observed, e.g., at the indices #77 and #103. However, the estimated centers (indicated by the ellipses) are located about 500-m apart with the mean uncertainty of 10 and 11 cm/s, respectively, in the derived current field. Due to the averaging nature of DTT measurements, the circulation centers are uncertain using the DS method.

To improve the estimate of the vortical structure, the DTT data with an SNR higher than 15 dB are used in an inversion employing the general LS method assuming the field is frozen during the collection period. Figure 3(e) shows the estimated field (arrows) along with the uncertainty (background color). The presence of the vortical structure is consistent with the simulation result of an ocean model based on the Princeton Ocean Model. The smallest uncertainty is located where there are concentrated crossing rays, while the largest uncertainty is observed near the tomographic boundary. The vortical structure is resolved with an uncertainty of 7 cm/s near the center. Outside the mouth of the Bay, the estimated currents are different from the ADCP-measured currents since the general LS method with the Fourier representation tends to have better wavenumber resolution, this results in a poor estimation of a uniform current field, with a large root-mean-square difference (RMSD) of 0.62 m/s in both eastward and northward components along legs #2 and #3. Further reduction of the RMSD is expected when combining the ADCP-measured currents into the inversion process.

This study has developed the technique of MVT to extend the survey area covered by autonomous vehicles and obtain a comprehensive map of the ocean current distribution. Via the DS method, with the measurements of the DTT and the vehicles' velocities projected along the transmission path, the path-averaged currents can be computed for estimation of the areal current velocity. The estimated areal currents near the ship show good agreements with the ADCP point measurements when the angle between two acoustic transmission paths is larger than $20°$⁠. Via the general LS method, a map of an anticyclonic circulation with a relatively small uncertainty is obtained using all the high-SNR rays collected during the experiment. Future work includes combining both integral data and point measurements for reducing the uncertainty of estimated ocean temperature and currents, and using the real-time estimation of these ocean properties for path-planning of a vehicle fleet.

This study was supported by grants from the Ministry of Science and Technology of Taiwan (MOST 108-2638-M-002-002-MY2 and MOST 108-2221-E-002-061). The bathymetric data are provided by the National Museum of Marine Science and Technology, Taiwan.