Estimating relative channel impulse responses from ships of opportunity in a shallow water environment

The uncertainty of estimating relative channel impulse responses (CIRs) obtained using the radiated signature from a ship of opportunity is investigated. The ship observations were taken during a 1.4 km (11 min) transect in a shallow water environment during the Noise Correlation 2009 (NC09) experiment. Beamforming on the angle associated with the direct ray-path yields an estimate of the ship signature, subsequently used in a matched filter. Relative CIRs are estimated every 2.5 s independently at three vertical line arrays (VLAs). The relative arrival-time uncertainty is inversely proportional to source bandwidth and CIR signal-to-noise ratio, and reached a minimum standard deviation of 5 ls (equivalent to approximately 1 cm spatial displacement). Time-series of directpath relative arrival-times are constructed for each VLA element across the 11 min observation interval. The overall structure of these time-series compares favorably with that predicted from an array element localization model. The short-term standard deviations calculated on the direct-path (7 ls) and bottom-reflected-path (17 ls) time-series are in agreement with the predicted arrival-time accuracies. The implications of these observed arrival-time accuracies in the context of estimating sound speed perturbations and bottom-depth are discussed. VC 2018 Acoustical Society of America. https://doi.org/10.1121/1.5052259


I. INTRODUCTION
Controlled sources enable many underwater acoustic systems such as active sonars, undersea acoustic communication, and various imaging and inversion methods.Their role is the emission of known, accurately timed signals allowing their separation or deconvolution from the channel impulse response (CIR).The pervasive presence of ocean ambient noise along with the environmental and cost issues of using active acoustic sources motivates considering how noise sources can be used in many of these same applications.
2][3] Further, lack of uniformity in the spatial distribution of sources can make this approach impractical.
In this paper, we investigate the feasibility of using the radiated signature from ships in place of controlled sources.The single most challenging aspect of this approach is the timing uncertainty of the estimated CIRs.For example, given a known source position and transmit time, absolute travel times measured over several kilometers require accuracy of order milliseconds for estimating ocean sound speed perturbations within 1 m/s. 4 In contrast, differential (relative) arrival-times between multiple hydrophones or CIR multipaths may be used in place of absolute travel times if the source timing or position have a much larger uncertainty than the travel time variations of interest.Using just the path differences cancels the common error but estimating sound speed perturbations of the same order (i.e., 1 m/s) increases the requirements on relative arrival-time accuracy.
We show that this accuracy is attainable from ships of opportunity together with ship position information provided by the Automatic Identification System (AIS) used as a constraint.An enabler of this method is that moving ships provide the equivalent of many high amplitude source transmissions.As an example of demonstrating this accuracy, we show that array element localization on the order of centimeters can be obtained.In the absence of significant systemic uncertainty, the CIR estimate contains information on the geoacoustic characteristics of the waveguide including water column sound speed structure, seafloor depth, and subbottom characteristics.
This paper is organized as follows.In Sec.II, we discuss the processing steps to separate the source signal from the CIR and define arrival-time uncertainty.Processing of the Noise Correlation 2009 (NC09) experiment data set and the development of an array element localization (AEL) model is presented in Sec.III.A discussion of the results is given in Sec.IV, followed by conclusions in Sec.V.For the interested reader, the appendixes provide additional detail on beamforming/sparse Bayesian learning and figures illustrating signal-to-noise ratio (SNR), integration time, and uncertainty of the direct-path and the bottom-reflected-path arrival-times.

II. SHIPS OF OPPORTUNITY AS SOURCES
Using a moving and controlled source for the purpose of inference has been of interest for decades. 5,6Transiting a) Electronic mail: gemba@ucsd.edu2][13][14][15] In contrast to controlled sources, using the radiated signature from ships of opportunity requires additional processing because the broadband signal emitted by the anisotropic radiator 16 is unknown (note that approaches exist for modeling and classification [17][18][19] ).The task is to separate the unknown source signature from the unknown CIR.
Blind deconvolution attempts to separate the source signal from a set of diverse CIRs by use of an optimization.The terminology "blind" is used because both the signal and CIRs are unknown, which an algorithm must learn simultaneously. 20urther, the source signal is not used explicitly in the system (i.e., CIR) identification process. 21To render the problem tractable, statistical models of the input signal or the CIRs may be used as constraints.However, the presence of environmental noise and uncertainty (e.g., sensor position) can make this approach problematic from a practical standpoint.
In the context of ships as random radiators, one approach to constraining the problem is to use AIS information to identify and model the direct ray-path angle-of-arrival at a vertical line array (VLA).Beamforming on this angle yields an estimate of the source radiated signature. 22This information reduces the complexity of the problem to a known signal in noise approach.The locally estimated signal can be used in a matched filter in order to estimate the CIR relative to the unknown travel time from the ship to the reference element of the array (termed relative CIR).
A. Estimating the source signal and relative CIR Assume a vessel emits a random signal s(t) at position rs .Its Fourier transform over a segment f of a larger observation interval is Sðf Þ ¼ jSðf Þje iU s ðf Þ , where U s (f) denotes the phase of the source at frequency f.The duration of the segment will be referred to as the integration time in units of seconds.To simplify notation, dependence on f will be carried only for relative CIRs.The signal s(t) is convolved with a CIR and embedded in additive environmental noise N yielding the signal p n (t) observed by the nth VLA element located at rn .In the frequency domain, p n (t) is represented as where G n ðr n ;r s ; f Þ is the Green's function (Fourier transform of the CIR) between the source and the nth receiver.We drop the dependence on N n ðf Þ since the ship signature substantially dominates the received signal.An estimate of S(f) is obtained by beamforming on the vertical angle h ¼ h D associated with the well-resolved, direct ray-path, The beamforming phase shifts are denoted by a Ã n ðh; f Þ [see Eq. (A1)] and the asterisk denotes complex conjugate.T denotes the travel time from the source to the reference element of the array.Hence Ŝðf Þ is a time delayed version of S(f).
To separate the source signature from the CIR, a matched-filtering technique is used, The phase [22][23][24] from Eq. ( 2) argð Ŝðf phone records received at the VLA in Eq. ( 1).The phase rotation removes the dependency of P n (f) on the source phase and provides an estimate of the Green's function relative to the unknown travel time T.Although T evolves with time as the ship moves, its influence is removed in Eq. ( 3) resulting in the direct ray-path arrival used to estimate S(f) in Eq. ( 2) becoming the reference point for the Green's function estimate Ĝn ðr n ;r s ; f Þ in Eq. ( 3).This estimate is shaded by the magnitude spectrum of the source, which the denominator in Eq. ( 3) attempts to remove. 22,25he output of the matched filter in Eq. ( 3) subsequently is transformed into the time domain using an inverse Fourier transform yielding an estimate of the relative CIR ĝn ðt; s n ; fÞ.The nth relative CIR depends on time t and relative arrival-time s n .In this work, s n is referenced to the direct-path arrival at the VLA reference hydrophone.
The dependence on integration time f is illustrated using an example.Consider beamforming on a direct ray-path angle h D for a non-transiting source.f controls the length of the estimated source signal [i.e., the time domain equivalent of Eq. ( 2)] used in a matched filter [Eq.( 3)].For a transiting source, however, actual CIRs change with source position.Hence, the choice of f influences both the SNR of the estimated relative CIR and smearing of the relative CIR due to source motion.This approach is comparable to increasing the duration of a controlled signal transmitted through a peak-power limited system.
Note that Eq. ( 2) provides a local estimate of the source signal at the reference element of the VLA and includes the effects of a combination of parameters such as Doppler, Lloyd's mirror interference pattern, the anisotropic radiation pattern of the vessel, and signal distortion due to propagation.In this work, Eq. ( 2) is used separately at each VLA.This approach acknowledges the dependence of Ŝðf Þ on those parameters.As a result, the CIR relative arrival-times between arrays are not jointly referenced to a single location using a single source signature estimate.

B. Predicted CIR arrival-time uncertainty
Matched filter CIR arrival-time uncertainty is predicted with respect to the direct-path or bottom-reflected-path.This uncertainty is inversely proportional to the root-mean-square (RMS) bandwidth of the source signal and SNR of the matched filter output: 4 The RMS frequency deviation or RMS signal bandwidth (Dx) rms is given by 26 ðDxÞ 2 rms ðx À where x and x 2 are called the mean and mean-square frequency deviations.
The SNR of a relative CIR peak of interest is defined as The amplitude of the peak of interest is denoted by ĝp ðt; s n ; fÞ and is determined by interpolation, along with its location.Assuming the peak of the envelope is best described by a parabola, a quadratic equation is used to fit the top three points.The location (time) and value (amplitude) of the maximum of this parabola then provide a better estimate of the true peak than the sampled (and quantized) data points alone.This method sometimes is referred to as subsample peak interpolation.The noise is estimated by the standard deviation of the matched filter output away from the peak and is denoted by r CIR .
We will show that these predicted arrival-time uncertainties [Eq.( 4)] are comparable with the observed shortterm variability of the time-series of relative CIR direct-path and bottom-reflected-path arrival-times (see Sec. IV).

A. Noise Correlation 2009 experiment
We use the Noise Correlation 2009 (NC09) experiment 27 data set for processing, recorded on Julian Day 31, 13:33-14:00 UTC.The experiment was carried out on the Coronado Bank 28 located southwest of Point Loma, San Diego, California, USA.Data were recorded at a f s ¼ 25 kHz sampling rate on three equally configured VLAs (see Fig. 1).
VLA 1 is the southeastern-most of the arrays and VLA 3 is the northwestern-most of the arrays.Assuming that VLA 1 is at the origin of coordinates with a bearing line of the arrays of 335 , the arrays are offset as noted in Table I, right-most column, along that bearing.
Representative hydrophone depths for VLAs 1-3 are shown in Fig. 2. Hydrophone 1 is 7 m above the seafloor and all N ¼ 16 elements on each VLA are used for processing.The element spacing of 1 m corresponds to a design frequency of 750 Hz at 1500 m/s.The mean sound speed over the arrays is c ¼ 1492 m/s and has been estimated using data obtained from nine CTD casts taken over nine experiment days and illustrate typically observed sea-surface warming effects.The surface mixed layer extends to a depth of approximately 25 m, followed by a downward refracting profile.
The 27 min track of the 52 m long R/V New Horizon is east of the VLA bearing line (see Fig. 1).The vessel travels northwest at a constant speed of 2.2 m/s and passes each VLA at approximately 100 m at the closest point of approach (CPA).

B. Data preprocessing
Two beamforming algorithms are used to preprocess the data in order to estimate relative CIRs.First, the direct raypath angle h D at each time step is estimated using sparse Bayesian learning (SBL).Appendix A 2 discusses the highresolution SBL method.The direct ray-path angle is used in Eq. ( 2) to estimate S(f) at the VLA reference element using conventional beamforming (CBF) over the bandwidth 0.07-2.5 kHz.See Appendixes A and A 1 for details.This source signal estimate subsequently is used to estimate the relative CIRs [see Eq. ( 3) below].
Relative CIRs are estimated approximately every 2.5 s for a total of i ¼ 1,…, I ¼ 270 realizations per element over the vessel transect of 1.4 km (11 min) from CPA at VLA 1 to CPA at VLA 3.An integration time of f ¼ 5 s is used for the results presented in Sec.IV.The source spectrum estimates overlap by 50% from one CIR realization to the next for a constant ray-path angle.Additional results presented in Appendix B investigate the effect of integration time  (1-44 s) on SNR [Eq.( 8)] and on uncertainty [Eq.( 4)] for the direct-path and the bottom-reflected-path arrival-times.

C. Inter-element time-delay estimation
The estimated relative CIRs ĝn ðt; s n Þj f¼5s can be used to observe the inter-element time-delays at each VLA.The direct-path relative arrival-times are assembled into a matrix ! 2 R NÂI .This matrix consists of N ¼ 16 time-series each with I ¼ 270 relative arrival-times.The reference element is hydrophone 16.These relative arrival-times are adjusted for the predicted moveout (the expected travel time differences for theoretical line array hydrophone positions without tilt).Let D 2 R NÂ1 be a column vector containing hydrophone distances to the reference element, and h D 2 R 1ÂI be a row vector containing the angles associated with the direct ray-path, where j is the array element index and i is the time index.Since the reference element is hydrophone 16 ðD 16 ¼ 0Þ; D 16i ¼ !16i ¼ 0 at every time step.The right-most term in Eq. ( 9) corresponds to plane wave beamforming time delays or phase shifts in the frequency domain, see Eq. (A1).Here, the direction of arrival angle h is negative when pointing upwards towards the sea surface.Thus, the computed matrix D 2 R NÂI captures relative arrival-times across the VLA adjusted for moveout for 270 ship positions.

D. Array element localization model
We hypothesize that actual VLA 2 hydrophone positions deviate from a simple line array.0][31][32] For example, an initial line array AEL model might assume that the elements are colinear with the array and their positions determined by a tilt angle and heading (direction of tilt) over time.
The NC09 data set provides an opportunity to analyze hydrophone position deviations from a simple line array configuration because of available array tilt and heading information recorded with engineering sensors mounted on each VLA.Heading information is included in the AEL model.Array tilt sensor data are not included explicitly in the model.During the period of data analysis, the mean tilt of VLA 2 was 0.5 with a standard deviation of 0.15 .
The AEL model investigates both horizontal and vertical displacement of the N ¼ 16 hydrophones from a tilted line array because the hydrophone umbilical cable may spiral around the tension-carrying strength member.Hence hydrophones can be displaced from their theoretical horizontal positions.The spiral also may cause a vertical displacement of hydrophones in addition to array tilt.
The unknown hydrophone positions may be determined successfully when interrogated with incoming plane waves arriving from different bearing angles (spanning a total of 180 is desirable).This is the case for a passing vessel when the approximate moveout associated with the CIR directpath relative arrival-times is removed [Eq.( 9)].
The AEL model attempts to explain the residual interelement time-delays in D. It creates a local, VLA-based reference frame and includes time-dependent parameters (such as VLA heading and ship position) significantly affecting the time-delays.Hydrophone positions are assumed timeinvariant with respect to this reference frame.
The right-most term in Eq. ( 9) attempts to remove significant dependence on h D of the inter-element time-delays.Effectively, all hydrophones are mapped onto the same horizontal plane and their positions are modeled using a polar coordinate system with radii r h 2 R NÂ1 and angles u 2 R NÂ1 , positive clockwise.Since the hydrophones are not exactly in that plane, the model includes parameters r v 2 R NÂ1 to account for a vertical offset for a total of three degrees of freedom per element.
The AEL-model of the inter-element time-delays where w 2 R NÂI captures the difference between two time-dependent angles u b and u h and the horizontal hydrophone angles u.The array element index is j and i is the time index [see Eq. ( 9)].The vessel-VLA 2 bearing angles u b 2 R 1ÂI are calculated using AIS (vessel GPS) and VLA 2 position information (see Fig. 1 and Table I).The VLA 2 heading angles u h 2 R 1ÂI vary 35 over the 11 min observation period.The sound speed c is 1492 m/s.The delay at hydrophone 16 is subtracted from T d , making it the reference element in T.
To estimate the unknown offset hydrophone parameters, the AEL-model inter-element time-delays T are compared to the relative arrival-times D in Eq. ( 9) using a vector valued, objective function ðT À DÞ 2 .Parameters are optimized using a non-linear, least-square algorithm (trust-region-reflective FIG. 2. Mean, first, and last sound speed profiles obtained from nine CTD casts during NC09.VLAs 1-3 are deployed at $150 m depth (see Table I).Instrumentation included a tilt, heading (compass), and depth sensor package at the top of each array.approach).Parameters r h ; r v , and u are constrained to 610 cm, 610 cm, and 0 -359 , respectively.The model output is only meaningful within the reference frame and does not provide absolute AEL.A detailed ray-arrival analysis (Fig. 4) is helpful to identify and extract direct ray-path angles (h D ) from 4-15 min for each VLA.The source (depth of 5 m) to hydrophone 16   34 The benefit of using sparse Bayesian learning (SBL) as a tool for a high-resolution angle estimate in the presence of coherent multipath appears evident when comparing results in Figs.4(d) and 4(h) to results using CBF in Figs.4(c) and 4(g).SBL processes approximately a 1 kHz bandwidth without suffering a loss in resolution or displaying aliased arrivals.This is useful if the VLA design frequency differs from the center frequency of the source (see also Appendix A 2). Angle identification is easier, especially when multiple ray-sets associated with multiple vessels are present or when identifying ray arrivals subject to an interference pattern [e.g., compare the BS interacting ray in panels (c)-(d)].

A. Data preprocessing and relative CIR estimation
The identified angles corresponding to the direct ray-path are used with conventional beamforming to estimate the source signature [Eq.( 2 An integration time of f ¼ 5 s is used to estimate these CIRs and a more detailed analysis of the effect of integration time on arrival-time accuracy is discussed in Appendix B. In particular, the standard deviation of the direct-path (bottomreflected-path) arrival-time uncertainty can be as small as 5 ( 8) ls for optimal integration times and by data addition over VLA sub-apertures.Integration time increases SNR up to a point (see Figs. 11-15).It is possible to obtain a ship signature estimate for a constant look-direction h D as long as the ray-path remains within the mainlobe of the beamformer.Thus, the integration time is limited by the rate of change of the direct ray-path [Figs. in Figs.7(a), 7(b1), and 7(c) for VLAs 1-3, respectively.It is apparent from these panels that the CIR relative arrivaltime structure varies approximately over 100 ls. 100 ls corresponds to a displacement of 15 cm at 1500 m/s.The shortterm standard deviations (2.5-10 ls) are calculated from 5 min, 10 s to 6 min, 10 s using the time-series of peaks corresponding to the direct-path D (discussed further below).One minute corresponds to 25 arrival-time peaks per hydrophone.Averaging standard deviations over all elements and three VLAs yields a root-mean-square (RMS) value of 7 ls.Uncertainties are comparable to predicted values using Eq. ( 4), illustrated in Fig. 11 and discussed in Appendix B.
Results in Fig. 7(b2) are further adjusted for VLA 2 AEL and also display standard deviations (4-8 ls) computed over all 270 (AEL-adjusted) direct-path arrival-times (4-15 min).The systemic uncertainty due to AEL is explained when comparing the data D [Eq. ( 9)] to the model T [Eq.(11)] timeseries.The relative CIR direct-path arrival-times (D) are shown in Fig. 8, black trace, for VLAs 1-3.The across-array mean is subtracted from the data, hence the 16th element is not a constant zero but carries the negative mean.In addition to the prior moveout correction, this adjustment helps to visualize the granular arrival-time structure across all the hydrophones.
The otherwise well-behaved relative arrival-time structure changes notably at CPA and is readily observed for VLA 2 (8 min), hydrophones 1-3.It appears a similar effect can be observed at CPA for VLA 1 (4 min) and VLA 3 (14 min, 30 s).This behavior points to the presence of AEL nuisance parameters and the data are compared to the output of the optimized AEL model T for VLA 2 (red trace).Indeed, this model correlates favorably with the time-series and mostly explains the trend.
Optimized AEL model parameters for VLA 2 are shown in Table II, along with standard deviations of the data and residual (AEL-adjusted data) time-series.Standard deviations are calculated over the entire time-series (4-15 min).The standard deviations of the residual (AEL-adjusted) time-series (bottom row of Table II) are on the order of centimeters.In comparison, Ref. 29 reports AEL uncertainties for horizontal (vertical) element displacement of approximately 20 cm (5-10 cm), which includes an inversion for ship position.

C. CIR bottom-reflected-path relative arrival-time structure
Panels (a)-(c) in Fig. 9 show a similar analysis as in Figs.7(a), 7(b1), and 7(c) but for the bottom-reflected-path arrivaltimes when the vessel position is close to VLA 3 (14 min, 30 s).These arrival-times are referenced to hydrophone 16 of the direct-path arrival and the moveout associated with the bottom-reflected-path is removed.Similar to the direct-path analysis, the CIR relative arrival-time structures vary over 100 ls.The short-term standard deviations are calculated from 14 to 15 min using the time-series of peaks corresponding to the bottom-reflected-path (similar to Fig. 8).
The increased standard deviations in Fig. 9(c) [37-122 ls] likely are due to the bottom-reflected angle being greater than the critical angle.Arrival-time uncertainty decreases when comparing results for VLA 3 to either VLA 2 in panel (b) [7-30 ls] or to VLA 1 in panel (a) [5-18 ls].Averaging standard deviations for all elements of VLAs 1 and 2 yields a RMS value of 17 ls.These short-term standard deviations are in line with predicted uncertainties illustrated in Fig. 15 and discussed in Appendix B. Results likely can be improved by removing systemic uncertainty due to AEL (not done here), considering the improvement from panel (b1) to panel (b2) in Fig. 7.

D. Application to making sound speed perturbation estimates
In this section we will demonstrate that arrival-time variabilities of the relative CIR bottom-reflected-path are the same size as those produced by sound speed perturbations.II).
The vessel is an uncontrolled source and its timing or position has a much larger uncertainty than the travel time perturbations of interest.Instead of estimating absolute travel times between the positions of the source and one receiver, relative arrival-times measured between multiple receivers can be used instead.However, the demands on arrival-time accuracy increase as the measurement distance decreases for estimating the same sound speed perturbations.
Accuracy requirements for absolute travel times (controlled source) or relative arrival-times can be evaluated analytically for the simple case of a uniform sound speed perturbation to a constant-velocity background sound speed.Let the reference sound speed be C 0 ¼ 1500 m/s with uniform perturbation Dc ( C 0 .The theoretical (or mean) travel time with respect to the reference sound speed from the source to the nth receiver is T n , and the travel time difference due to the sound speed perturbation is DT n .The slant range between the source and a receiver is R n .The measured travel time T 0 1 between the source and the first receiver is then where 1 is a realization of a random process with standard deviation r s [Eq.( 4)] and dT are correlated errors (e.g., timing errors for a single data acquisition system).We expanded T 0 1 around T 1 in Eq. ( 14) and neglect higher order terms [O(Dc 2 )] in the Taylor series expansion.Subtracting the mean travel time from Eq. ( 14) yields a first order estimate of the (absolute) travel time perturbation, If the source position and/or its timing is not known precisely, differential arrival-times taken between two receivers 4 can be used instead, which in this case cancels the common error dT.
The sensitivity with respect to the sound speed perturbation Dc (i.e., the quantity of interest) is a guide to what a given uncertainty corresponds to in terms of sound speed changes in this simple system with one unknown and one datum.The sensitivity for absolute travel times is R 1 =C 2 0 [Eq.( 16)] and is proportional to the full slant range.For ranges of several kilometers, uncertainties on order of milliseconds are acceptable for estimating sound speed perturbations within 1 m/s.
The sensitivity for relative travel times is proportional to the slant range difference [ðR 1 À R 2 Þ=C 2 0 in Eq. ( 18)].For the relative CIR direct-path, the distance represented by DR ¼ R 1 À R 2 is on the order of meters.Hence, if the distance decreases by a factor of 1000, travel-time accuracies on the order of microseconds instead of milliseconds are required because the sensitivity is of order ls per 1 m/s sound speed perturbation.
Specifically, one can observe in Fig. 5(d) that the travel time difference between the first and last hydrophone of the direct-path arrival is about 5 ms, which is consistent with DR ¼ 7:5 m.For this DR, a measurement having an arrivaltime uncertainty of 3 ls or less can produce an estimate with uncertainty of 1 m/s in a simple one-parameter inversion.In comparison, the range difference between the direct-path and bottom-reflected-path [Fig.5(d)] for hydrophone 16 is on the order of 20 ms (approximately 30 m).Relative CIR bottom-reflected-path arrival-time uncertainty must be 13 ls or less to produce a sound speed perturbation estimate with uncertainty of 1 m/s.While the actual environment is more complicated, this simple example indicates that it is useful to process multiple ray-paths to reduce the requirements on arrival-time accuracy.
The relationship between bottom-reflected-path relative CIR arrival-time variations and observed (extreme) sound speed perturbations is further investigated by comparing data to model results in Fig. 10.Travel times are calculated using Bellhop ray paths made from the minimum, maximum, and mean sound speed profiles observed during the NC09 experiment.The minimum and maximum (extreme) sound speed profiles correspond to the first and last CTD casts, respectively (see Fig. 2).Similar to the data, the arrivaltimes corresponding to the extreme and mean sound speed profiles are referenced to hydrophone 16 of the direct-path (i.e., the CIRs are relative CIRs).The travel times for the mean sound speed profile are removed from the extreme profiles.Hence the distance between the dashed line traces in Fig. 10 illustrates the travel time signature of observed minimum and maximum ocean sound speed perturbations and the possibility of estimating them from the observations.The data in Fig. 10 are also adjusted by the mean sound speed profile (this adjustment centers horizontal axes for relative CIR bottom-reflected-path arrival-times in between the minimum and maximum sound speed model results).In addition, all model arrival-times are adjusted by two ad hoc parameters, corresponding to a tilt angle and a travel time delay.The angles (0.9 and 0.6 , respectively, for VLAs 1 and 2) correspond to a tilted VLA, the delays (20 ls and À10 ls, respectively) correspond to a change in bottom depth.To allow for a comparison, the angle and delay for VLA 3 are selected such that the arrival-times approximately are centered at 0 ls.
For the direct-path arrival, analytical results indicate that an uncertainty of 3 ls would allow for distinguishing to within 1 m/s sound speed perturbations.Results in Figs.7(a), 7(b1), and 7(c) indicate an RMS error of 7 ls.While the AEL model [Fig.7(b2)] can correct for some of the systemic uncertainly, it does not fully capture the mean computed over the entire time-series (see Fig. 8, VLA 2, hydrophone 16).This may indicate that a more complicated model is required (e.g., with curved wavefronts around CPA).From an ocean state estimation viewpoint, the ocean travel time signature ("the signal") is small when compared to the uncertainty ("the noise") for the direct-path arrival.
The sensitivity to ocean perturbations increases when processing additional multipaths.The three panels in Fig. 10 show arrival-time uncertainties for VLAs 1-3.In panel (c), these variabilities are much larger than in the other panels and cannot be explained by ocean sound speed perturbation.The arrival-time variabilities shown in panels (a), (b) are the same size as those produced by sound speed perturbations.Specifically, model results in Fig. 10(a) indicate arrival-time perturbations of approximately 125 ls between minimum and maximum sound speed profiles observed during the experiment.Relative CIRs have a RMS error of 17 ls.The bottom-reflected-path yields information about the geoacoustic properties (e.g., bottom-depth) and the relative CIR has sufficient sensitivity to distinguish between the minimum and maximum observed sound speeds.

V. CONCLUSIONS
We demonstrated that relative channel impulse responses (CIRs) estimated from unknown but high SNR ship signature sources provide arrival-time accuracies on the order of 10 ls.The SNR and bandwidth of vessels allow for relative CIRs to be estimated with integration times of seconds, which is smaller than many processes contributing to the environmental variability.The radiated signatures of transiting ships provide a method for monitoring hydrophone positions (AEL) for vertical or horizontal line arrays.
While results demonstrate sensitivity on the order of 10 ls, the number of free parameters (e.g., hydrophone positions, bottom slope, depth, sound speed, etc.) significantly affect the arrival-time estimates.Systemic uncertainty due to AEL must be learned first and possibly tracked.The relative CIR direct-path provides marginal utility by itself for estimating sound speed perturbation.The relative CIR bottom-reflected-path (referenced to the direct-path) provides sufficient sensitivity to distinguish between the minimum and maximum sound speed arrival-time variabilities observed during the NC09 experiment.Using AIS as a constraint, geoacoustic parameters such as bottom depth over a region also can be learned using passing vessels with different transects.
coherent multipath arrivals at a VLA.When CS is implemented using multi-frequency sparse Bayesian learning (SBL) 35,36 as is the case here, sparsity (i.e., the number of unknown ray-paths) is determined automatically.SBL exhibits properties similar to the adaptive white noise constraint processor and is robust to mismatch. 37,38Next is a short overview of SBL, see also Ref. 38 for additional details on the identical code implementation.

a. Sparse Bayesian learning overview
The lth data snapshot y l can be expressed with an underdetermined system of linear equations, In Eq. (A5), A ¼ ½aðh 1 Þ; …; aðh M Þ is the dictionary containing M replica vectors and n l $ CN ðn l j0; r 2 I N Þ is complex Gaussian noise with noise variance r 2 .x l is the vector of complex source amplitudes with entries corresponding to the same angles as in h.SBL models the unknown source amplitudes as complex Gaussian random variables with prior density pðx l Þ ¼ CN ðx l j0; CÞ, where C is a diagonal covariance matrix, i.e., C ¼ diagðc 1 ; …; c M Þ ¼ diagðcÞ.The vector c is the source power for all h.Since the noise is Gaussian, the likelihood is expressed as pðy l jx l ; AÞ ¼ CN ðy l jAx l ; r 2 I N Þ.
Denote the collection of L snapshots at frequency f as Y f ¼ ½y f ;1 ; …; y f ;L .Let the corresponding collection of source amplitude vectors and dictionaries be denoted by X f and A f , respectively.Then where N f are additive noise contributions.For approximately stationary sources, x l,f are assumed independent over time.Hence, we have where C f ¼ diagðc f Þ is the covariance of the source amplitude.We assume that X f are independent for F processed frequencies.There is no assumption of sparsity in the frequency domain, which makes this formulation attractive for localizing a few broadband sources from many candidate replica vectors.A way to enhance sparsity of the solution is to set This prior model assumes that X f has the same statistical distribution with covariance C at each frequency.A sparse C would impose identical sparsity constraints on the vectors contained in X 1 ; …; X F .Because vectors contained in X 1 ,…,X F and N 1 ,…,N F are independent, the joint evidence p(Y 1 ,…,Y F ) over all frequencies is where the model covariance R We maximize the joint evidence by equating its derivative to zero, which gives the iterative update rule Since ĉnew m ðhÞ represents source power, its plot can be interpreted as a broadband beamformer output [see Figs.signal remains within the mainlobe of the beamformer.For a moving source, observed CIRs must be similar over this duration.The effect of integration time f on relative CIR direct-path arrival-time estimation accuracy is shown in Fig. 11.Results for each VLA are averaged using RMS (i.e., over 16 CIRs).The minima at VLAs 1-3 are 5, 6, and 5 ls, respectively.SNR is the main driver decreasing uncertainty of the direct arrival (Fig. 12), while RMS bandwidth of the source roughly remains a constant.
Note that integration time also can be interpreted as vessel movement (the velocity is 2.2 m/s), ranging from approximately 2.2 to 100 m.Example SNR and corresponding uncertainty using an integration time f ¼ 5 and f ¼ 20 s (corresponding to the white and black dashed lines in Fig. 11) are shown in Figs. 13 and 14.An integration time of f ¼ 5 s is used in Sec.IV, the predicted uncertainties are shown in Fig. 13(e), dotted trace.When time-series are grouped into sub-apertures and added, uncertainty further decrease from 10 to 5 ls in Fig. 14(d).Applying this processing before estimating relative CIRs [Eq.( 3)] primarily increases SNR and reduces the number of spatial observations.
The effect of integration time on relative CIR bottomreflected-path arrival-time estimation accuracy is shown in Fig. 15.The examined CIRs are the same as in Fig. 11, the dynamic range is extended to 50 ls.For identification, the bottom reflected peaks are assumed to be delayed by 1-35 ms when measured from the direct arrival in each channel.Similar to the direct-path arrival results in Fig. 11, the bottom-reflected-path relative arrival-time uncertainty decreases with increasing integration time.The minima at VLAs 1-3 are 8, 10, and 15 ls, respectively.Uncertainty increases most notably when the vessel is close to CPA, which corresponds to a grazing angle of approximately 28 and is comparable to results using surface noise to estimate the critical angle [see Fig.

FIG. 1 .
FIG. 1. (Color online) Bathymetry map of Coronado Bank off the coast of San Diego showing 3 VLA positions and the 27 min ship track marked every 4 min.The vessel travels northwest along the 150 m isobath and passes the most southeastern VLA 1 at 4 min.
Time-synchronized spectrograms of VLAs 1-3 [Figs.3(a)-3(c)] display the frequency vs time spectral structure of the recorded data.The received acoustic power at hydrophone 16 indicates a broadband signal observable at CPA for VLAs 1-3 at 4, 8, and 15 min, respectively.An interference pattern is visible before and after CPA in addition to a significant change in received power (about 20 dB over an interval of several minutes).These recorded data contain both the unknown CIR and the unknown source signature.

FIG. 3 .
FIG. 3. (Color online) Time-synchronized spectrograms showing sound pressure spectrum levels received at hydrophone 16 for (a) VLA 1, (b) VLA 2, and (c) VLA 3. The CPA of the vessel at each VLA is marked with a red arrow and occurred at 4, 8, and 15 min, respectively.The top horizontal axes display vessel-VLA distance, negative units indicate an approaching vessel.
eigenrays are computed with a range dependent Bellhop 33 model [Figs.4(b), 4(f)].We use a range dependent bottom obtained from the NOAA NGDC coastal digital elevation model.
)] and relative CIRs at each hydrophone [Eq.(3)].A comparison between the Bellhop eigenray arrival-times and the estimated relative CIRs ĝn ðt; s n Þ is shown in Fig. 5 at 5 min, 40 s (vessel-VLA 1-3 distances of 230, À280, and À1130 m, respectively, see Fig. 3).Both show a similar multipath arrival structure.A delayed arrival is visible for each arrival pair due to reflection off the surface, producing a Lloyd's mirror effect.VLA 3 is at longer range than the other VLAs and additional multipaths (BS, BSB) are present.The time-axes in panels (a)-(c) show absolute travel time from the position of the source to the VLAs.The relative arrival-times (s n ) for the estimated CIRs [panels (d)-(f)] are referenced to the direct-path arrival at hydrophone 16, so that arrival-time is identically 0 s.The observed ship signature [Eq.(2)] contains the dipole (Lloyd's mirror) effect, the anisotropic vessel radiation pattern and multipath-dependent Doppler shift due to motion.These transfer functions depend on time and on source-receiver geometry.In addition, the Lloyd's mirror is not present at close ranges for the bottom-reflected-path [Figs.5(d)and 5(e)], possibly because the vessel is not wellcharacterized as a point source.The relative CIR estimates at hydrophone 16 for the data in Fig.5are displayed in Figs.6(a)-6(c) with details of the direct-path arrivals in panels (d)-(f).The predicted timing accuracies [Eq.(4)] of the direct-path peaks are 16, 18, 32 ls, respectively.The bow (stern) of the vessel is approximately pointed towards VLA 2 (1).Panel (c) displays increased variability away from the peaks (r CIR ) when compared to panels (a), (b).The out-of-phase arrival at 5 ms in panel (c) is separated by 30 ls from the peak of the bottomreflected-path.
FIG. 5. (Color online) Bellhop eigenray arrival-times vs hydrophone number for the geometry at 5 min, 40 s (CIR set no. 42 of 270) showing multipath arrivals in panels (a)-(c) for VLAs 1-3, respectively.The delayed arrival pair is due to the Lloyd's mirror effect for the 5 m deep source.The vessel is positioned approximately equal distance between VLAs 1 and 2. Panels (d)-(f) show relative CIRs.The geometry is the same as in panels (a)-(c) but with the direct-path relative arrival-time s n at hydrophone 16 set to 0 s.The integration time f ¼ 5 s, the frequency band is 0.07-2.5 kHz.

FIG. 6 .
FIG.6.Panels (a)-(c) show normalized CIRs [hydrophone 16 in Figs.5(d)-5(f)] for VLAs 1-3, respectively.Right panels show details of the direct-path arrival with predicted peak accuracies of 16, 18, and 32 ls, respectively.Because the peaks have been normalized to 1.0, it is not possible to see the varying peak amplitude that determines the accuracy.

FIG. 8 .
FIG. 8. (Color online) Columns 1-3 show the direct-path relative arrival-times of 270 CIR sets (D) in black for VLAs 1-3, respectively.The red trace for VLA 2 is the output of the optimized array element localization model (T).Data are adjusted for the time delays associated with the direct-path arrival angle.The across-array mean is subtracted from the data, hence the 16th element is not a constant zero.The range of the left vertical axes are 690 ls, 100 ls corresponds to 15 cm displacement.CIR examples shown in Figs.5-7 correspond to CIR set No. 42 at 5 min, 40 s.The CPA of the vessel at each VLA is marked with a red arrow.
FIG. 11. (Color online) Panels (a)-(c) display RMS direct-path CIR relative arrival-time uncertainty r t for VLAs 1-3, respectively.The horizontal white (f ¼ 5 s) and black (f ¼ 20 s) lines mark the integration time used in Figs. 13 and 14, respectively.The CPA of the vessel is marked with a red arrow.

TABLE I .
VLA locations and separation from VLA 1.

TABLE II .
Optimized AEL model parameters and standard deviations (SDs) calculated over the entire 11 min time-series (4-15 min) for VLA 2 (see Fig.8).The three model parameters are discussed in Sec.III D. 7 ls correspond approximately to 1 cm displacement.