A general blind deconvolution algorithmic framework is developed for sources of opportunity (e.g., ships at known locations) in an ocean waveguide. Here, both channel impulse responses (CIRs) and unknown source signals need to be simultaneously estimated from only the recorded signals on a receiver array using blind deconvolution, which is generally an ill-posed problem without any a priori information or additional assumptions about the underlying structure of the CIRs. By exploiting the typical ray-like arrival-time structure of the CIRs between a surface source and the elements of a vertical line array (VLA) in ocean waveguides, a principle component analysis technique is applied to build a bilinear parametric model linking the amplitudes and arrival-times of the CIRs across all channels for a variety of admissible ocean environments. The bilinear channel representation further reduces the dimension of the channel parametric model compared to linear models. A truncated power interaction deconvolution algorithm is then developed by applying the bilinear channel model to the traditional subspace deconvolution method. Numerical and experimental results demonstrate the robustness of this blind deconvolution methodology.

Acoustic remote sensing methods typically use a known source signal to probe the environment, and then directly estimate the time-domain channel impulse responses (CIRs)—also referred to as the time-domain Green's function—from the measurements obtained on a receiver array, which are modeled as the convolution of the source signal and the CIRs when the channels are time invariant during the observation period. In this paper, for the sake of simplicity, the CIR refers to the band-limited version of the infinite bandwidth CIR within the bandwidth defined by the source and receiver frequency responses. However, in many situations, a known active source signal is not available because of either environmental concerns or hardware constraints, while only sources of opportunity (e.g., radiating ships at known locations) are available. Hence, sources of opportunity offer an attractive alternative for acoustic remote sensing applications (e.g., such as passive acoustic thermometry or ocean acoustic tomography, as well as geoacoustic inversions). The main challenge is that the precise radiated waveforms (i.e., acoustic signatures) of these opportunistic source signals are unknown; therefore, both the CIR and the unknown source signal need to be simultaneously estimated from only the recorded waveforms on a receiver array using blind deconvolution, which is generally an ill-posed problem if no a priori information or additional assumptions about the underlying structure of the CIRs are available. In ocean waveguide environments, the typical ray-like arrival structure of the broadband CIRs wavefronts between a surface source and the elements of a vertical line array (VLA;Godin, 2007; Roux et al., 2008) provides such a priori information about the CIRs, which can be used to turn this ill-posed problem into a feasible one. This paper presents a general algorithmic blind deconvolution framework for acoustic remote sensing using sources of opportunity and a priori information about the arrival-time structure of the CIRs in an ocean waveguide.

The general blind deconvolution problem under a multiple-channel receiver framework (hereafter referred to as the multichannel blind deconvolution problem) has been studied in various areas, and the selected solutions by previous studies are usually application dependent. In digital communication applications, multichannel blind deconvolution methods do not require the transmission of a training sequence, and they are adaptive to fast varying channels (Giannakis and Serpedin, 1997; Moulines et al., 1995; Xu et al., 1995). In image restoration, image denoising, image deblurring, and medical image reconstruction, multichannel deconvolution methods aim to recover the image (source) from convolutions with multiple unknown finite impulse response (FIR) filters (Katsaggelos, 1991; Souidene et al., 2009). Applications in seismology (Mendel, 2013) have also been reported. In underwater acoustics applications, the multichannel blind deconvolution problem has also been investigated using various approaches and assumptions on the CIR features (Abadi et al., 2012; Broadhead and Pflug, 2000; Byun et al., 2017; Cazzolato et al., 2001; Roan et al., 2003; Sabra and Dowling, 2004; Sabra et al., 2010; Smith, 2000). A subset of the authors of this study have previously developed a fast multichannel myopic deconvolution method via a low-rank recovery technique (Tian et al., 2017), in which the CIRs for channels were independently parametrized. However, this deconvolution method via low-rank recovery did not take into account the inherent relationships between the CIRs across all channels (e.g., such as the arrival-times of each CIR along a VLA being aligned along common ray-like wavefront arrivals), which can be exploited as a further constraint for jointly parametrizing the CIR between the source of opportunity and the receiver array elements for the selected channel model (here, an ocean waveguide).

In order to develop a more robust multichannel blind deconvolution method for ocean waveguides, this paper leverages the inherent linking of the amplitude and arrival-times of the CIRs across all channels of a VLA to develop a bilinear channel model for the CIRs. This approach further reduces the dimension of the channel parametric model compared to the previously proposed linear model (Tian et al., 2017) and, hence, ultimately improves upon the robustness of the proposed multichannel blind deconvolution method. To justify the selection of this bilinear channel model for the CIRs, it can be noted that the arrival-time structure of the broadband CIRs is composed of stable ray-like wavefront arrivals between a surface source and the elements of a VLA (Godin, 2007; Roux et al., 2008). The primary envisioned application for this method is the use of surface sources of opportunity [e.g., radiating ships with known locations inferred from the automatic identification system (AIS); Verlinden et al., 2015] to perform acoustic tomography (or thermometry) of the ocean, as well as geoacoustic inversions, passively. In this context, it is expected that prior information of the amplitude and arrival-time structure of the CIR is available and can be obtained either using a sufficiently accurate numerical model that describes the baseline acoustic propagation in the selected environment (as it is done for acoustic tomography applications) or from in situ measurements (collected in the vicinity of the actual location of source of opportunity of interest) using the ray-based blind deconvolution (RBD) method (Abadi et al., 2012; Byun et al., 2017; Sabra et al., 2010). In either case, it is expected that the proposed multichannel blind deconvolution method will provide a more accurate estimate of the actual CIR for the selected source of opportunity than the prior information available obtained from the CIR parametrization. The multichannel blind deconvolution problem solved in this paper uses a truncated power iteration deconvolution algorithm specifically developed for a constrained eigenvalue decomposition given by the cross-convolution formulation for solving the multichannel blind deconvolution problems, which has been researched extensively in the 1990s (Gürelli and Nikias, 1995; Moulines et al., 1995; Xu et al., 1995). However, both the empirical performance and theoretical guarantees of these classical methods were restricted to the noiseless case, and these methods are not robust in the presence of additive measurement noise. We overcome this limitation by imposing data-adaptive structural models in the estimation. The proposed numerical implementation is both fast and scalable. Similar approaches to model-based multichannel blind deconvolution have been developed in the context of various imaging modalities and image processing applications (Kazemi and Sacchi, 2014; Morrison et al., 2007; She et al., 2015; Sroubek et al., 2007; Sroubek and Milanfar, 2011). They proposed sophisticated iterative algorithms that provide estimates for both the source signal and the channels by alternating between optimization programs (She et al., 2015; Sroubek and Milanfar, 2011), taking advantage of expected sparse structure in one or both parts of the optimization. In this paper, we are interested only in recovering the channel responses and use a completely different model for their structure. We model the channels as jointly living in a low-dimensional subspace derived from a priori information about the CIRs (e.g., a window of arrival-times—the result is a “block rank-1” model for the channels. The algorithm we provide below is simple in that it estimates the channels only once, but, in principle, it could be used as a sub-step in an alternating optimization algorithm such as those in the references above. As we demonstrate in the numerical result section below, though, the algorithm already provides satisfactory reconstruction results. Thus, in order to keep the computational cost low, we did not pursue further iterative refinements. Finally, the proposed approach notably differs from previous sparse blind deconvolution techniques, which typically consider the single channel blind deconvolution problem where the CIR is sparse in the canonical basis (Choudhary and Mitra, 2012, 2014; Li et al., 2016). Instead, the proposed blind deconvolution technique in this paper relies on a different model than the sparsity model by using a low-dimensional subspace model, which is already leveraging the extra a priori information contained in the CIRs (e.g., specific arrival-time structure). As we explained above, this subspace model fits well to the data measured in an underwater acoustic channel and considered in this paper. Furthermore, even when using this extra a priori information, accurately estimating the arrival-times structure of the actual CIRs (e.g., for tomography purposes) under a low signal-to-noise ratio (SNR) is still quite challenging in practice.

The remainder of this paper is organized as follows. Section II presents the bilinear channel model—motivated by the acoustic remote sensing scenario using a VLA—and the formulation of the blind deconvolution problem using the cross convolution method under the bilinear channel representation. Section III introduces the truncated power iteration method algorithm and discusses the initialization step. Section IV then presents results for the case of bottom-mounted short VLAs in shallow water waveguides using either numerical simulations or in situ CIR estimates obtained experimentally with the RBD method to parametrize the CIRs arrival structure. Finally, Sec. V summarizes the findings of this paper.

In this section, we first introduce a bilinear channel model, which is motivated by an acoustic remote sensing scenario using a VLA listening to a source of opportunity in the far-field. A full theoretical analysis of this bilinear recovery problem can be found in Lee et al. (2018). We then mathematically recast the multichannel blind deconvolution problem as an optimization problem.

In order to illustrate the concept of the bilinear channel model, we first consider the simple case where the wavefronts composing the CIRs between the source of opportunity and the receiver array elements can be well approximated by propagating rays (i.e., for sufficiently high frequencies; Jensen et al., 2011). In the far-field of the source for any given ray, the individual arrivals along each wavefront can be considered as time-shifted and amplitude-scaled replicas of each other. The incident angle of this incoming ray is determined by the source-receiver configuration (i.e., VLA geometry and source and receiver range), as well as environmental parameters [e.g., local sound speed profile (SSP) at the VLA]. For instance, for the case of a single ray arrival, the CIR for channel m is written as

hm(t)=amp(ttm),
(1)

where p(t) denotes the baseline pulse waveform [Fig. 1(a)] that is determined by the frequency response of the receivers, and tm [Fig. 1(c)] and am [Fig. 1(b)] are, respectively, the values of the arrival-time and amplitude of the ray arrival for channel m. While the true channel response is an impulse amδ(ttm) (or, more generally, a linear combination impulse at different delays), we are combining it with the (known) receiver response p(t) to get our CIR model [Eq. (1)]. As we will see below, this allows us to build an efficient data-adaptive structural model on the channels. Furthermore, the arrival-time variations and amplitude variations of this ray arrival, shown in Fig. 1, were arbitrarily selected to represent an idealized ocean waveguide arrival for the sake of illustration).

FIG. 1.

(a) (Color online) Baseline pulse waveform used for the numerical simulations. (b) Variations of the received amplitude parameters am for a single ray assuming 16 receivers. (c) Variations of the ray arrival-times parameters tm. (d) Synthesized ray arrival using Eq. (1) with the parameters shown in (a)–(c).

FIG. 1.

(a) (Color online) Baseline pulse waveform used for the numerical simulations. (b) Variations of the received amplitude parameters am for a single ray assuming 16 receivers. (c) Variations of the ray arrival-times parameters tm. (d) Synthesized ray arrival using Eq. (1) with the parameters shown in (a)–(c).

Close modal

We assume that the arrival-times across channels are linked by a function tm=f(m), where the function f describes the relation of the arrival-times across channels. Then, we can build a parametric model for the CIRs using this relation. For instance, when the source of opportunity is in the far-field of the array and the local SSP along the array is quasi-constant, it can be assumed that the arrival-time tm falls along a line whose slope k is determined by the ray incident angle (Jensen et al., 2011). Furthermore, if the receivers are equally spaced in a vertical array, this relationship can be expressed as

tm=t1+(m1)k,
(2)

where t1 is the arrival-time for the first channel, and k is the slope of the line in Fig. 1(c).

Combining Eqs. (1) and (2) forms the basis for building a bilinear model for a given ray arrival of the CIR using one variable vector a=[am] to represent the discrete ray amplitudes am and one variable vector t=[tm] to represent the linked ray arrival-times tm. Furthermore, we also discretize each time-shifted pulse profile p(ttm) [see Eq. (1), m=1,,M] recorded on the M receivers into a vector pm of length K (the number of samples). Since the channels are treated jointly in the bilinear model, we concatenate the M discretized pulses pm into a single pulse basis vector p¯,

p¯=[p1,,pM]MK.
(3)

Similarly, we discretize each time-domain CIR hm(t) into a vector hm, where each arrival is determined by the product of the amplitude coefficient am and the pulse shape pm. [e.g., given by Eq. (1) for the ray-model].

The linear relation for arrival-times across channels was provided in Sec. II A as an illustration of a bilinear channel model, but it is not a necessary assumption in order for the method to work. Indeed, even if neither an explicit function f nor any other parametric function describing the relation of the arrival-times across channels is available a priori, a bilinear channel model can still be built by using a collection of realizations of CIRs (i.e., a library of a priori known CIRs obtained from experimental measurements or numerical simulations) to extract the actual pulse arrival-times, which are representative of the source-receiver configuration and variability of the environment of interest (e.g., potentially including the effects of the intermittent micro-multipath). Using this a priori information of arrival-times, the linear subspace containing the entire pulse vector p¯ [Eq. (3)] can still be learned, and a bilinear channel model can then be constructed. In ocean waveguide environments, the typical ray-like arrival structure of the broadband CIR wavefronts between a surface source and the elements of a VLA (Godin, 2007; Roux et al., 2008) provides such a priori information about the CIRs.

Specifically, if i indexes the specific CIR realization for a given set of environmental parameters and source-receiver configuration, the associated pulse vector is denoted as p¯i. For instance, using the ray-model from Sec. II A [see Eq. (1)] and assuming a priori information is available to bound the fluctuations of ray arrival-times [and, hence, slope k—see Eq. (2)] due to environmental fluctuations (e.g., using historical values of SSPs at the test site or more sophisticated dynamic ocean simulations), we can determine the time-window containing a collection of possible arrival-times ti (Fig. 2) by varying the ray slope k accordingly. In general, although the total number of realizations can be large, the individual realizations p¯i are highly correlated; hence, they lie on a lower dimension manifold, which can be approximated by a linear subspace using the principal component analysis (PCA) technique (Jolliffe, 2011) to reduce the dimensionality of the search space for the blind deconvolution method. We denote the basis matrix of this subspace as ΦMK×W (where M is the number of channels and K is the sample length of each discretized CIR), which is constructed by using the first W orthonormalized vectors in the PCA expansion, and for a certain set p¯, we can write it as

p¯=Φu,
(4)

where uW is the coefficient vector that fully describes the corresponding realization of arrival-times t.

FIG. 2.

The black shaded region indicates the range of possible arrival-times variations for sound-speed fluctuations on the order of ±3 m/s for the experimental source-receiver configuration shown in Fig. 4.

FIG. 2.

The black shaded region indicates the range of possible arrival-times variations for sound-speed fluctuations on the order of ±3 m/s for the experimental source-receiver configuration shown in Fig. 4.

Close modal

For instance, assuming a total of 26 realizations [arbitrarily generated by varying the coefficients in Eq. (2) to account for sound-speed fluctuations on the order of ±3 m/s for the experimental source-receiver configuration shown in Fig. 4] and using the ray-model for the CIR, the first W = 6 basis vectors of the PCA expansion are shown in Fig. 3.

FIG. 3.

(Color online) First six basis vectors of the PCA expansion of the library of time-shifted ray arrivals shown in Fig. 2.

FIG. 3.

(Color online) First six basis vectors of the PCA expansion of the library of time-shifted ray arrivals shown in Fig. 2.

Close modal

If all discretized vectors of time-domain CIR hm s [as defined after Eq. (3)] are concatenated as follows, the bilinear channel model can be written as

[h1h2hM]=[a1Φ1a2Φ2aMΦM]u,
(5)

where (Φm)m=1MK×W are the blocks of the subspace matrix Φ defined in Eq. (4), uW is the coefficient vector that determines the pulse arrival-times for one ray path across all channels, and each element of the vector a=[a1,,aM]M with am>0 is the arrival amplitude am. Furthermore, let h¯MK denote the concatenation of the CIR coefficients for all channels, i.e.,

h¯=[h1,,hM].
(6)

Equation (5) can then be equivalently rewritten as a more concise form using the Kronecker product as

h¯=Ψ(au),
(7)

where

Ψ:=[Φ1000Φ2000ΦM].
(8)

The formulation above only represents a single wavefront of the CIR. For a more general multipath environment (i.e., where the number of identifiable paths is N > 1), the bilinear channel model still holds, and the concatenation of CIR vectors h can be written as

h¯=Ψ̃[(a1u1),,(aNuN)],
(9)

where

Ψ̃:=[Φ1100ΦN1000Φ1200Φ22000Φ1M00ΦNM].
(10)

Equation (9) describes the general bilinear channel model that exploits the inherent pulse arrival-times relation across channels (e.g., along a single ray-arrival or any more complex wavefront). This bilinear model differs from the linear channel model studied in previous research (Tian et al., 2017) and further reduces the dimension of the subspace for all CIRs by constraining the model parameter space compared to the linear model. However, under this bilinear channel model, a new multichannel blind deconvolution algorithm for the CIR recovery is needed and is presented in Sec. III.

Our blind deconvolution method combines the model in Eq. (7) with a classical cross-convolution technique (Gürelli and Nikias, 1995; Moulines et al., 1995; Xu et al., 1995). Cross-convolution takes advantage of an invariance that is a straightforward consequence of the commutivity of convolution operations. In the noiseless case, we have yi=shi and, therefore, for any ji,

yihj=shihj=shjhi=yjhi.

Thus, the true unknown channels hi and hj must obey the homogenous linear system of equations

TyihjTyjhi=0,

where Tyi is the (Toeplitz) matrix whose action corresponds to convolving with the channel output yi. This relation holds for each of the M(M1)/2 pairs of channels; we can express them jointly using

Myh¯=0M(M1)L/2×1,
(11)

where My is the block-sparse matrix

My=[Ty2Ty10000Ty30Ty10000Ty3Ty20000Ty40Ty2000000TyMTyM1].

In other words, we know that when there is no noise, h¯ is a null vector of My. In Xu et al. (1995), it is shown that under a mild algebraic condition on the channels (namely, that they have no common null in their frequency response), the true channel responses h¯ are the unique (up to a scaling coefficient) vectors that obey this condition. In the presence of noise, My, in general, will not have a null space. So instead of directly computing for the null space of My, we estimate the channels by solving an approximate constrained problem to approximate this null space

minimizeh¯||Myh¯||22subjectto||h¯||2=1,
(12)

which, in turn, is equivalent to actually finding the eigenvector corresponding to the minimum eigenvalue of MyMy. When there is zero-mean white noise added to the channel outputs, we know that the estimate above is consistent (as the corresponding perturbation to MyMy will average out as the number of measurements L goes to infinity), but without additional modeling assumptions in place, it will be very sensitive to noise when working with a finite number of samples.

We implicitly regularize this problem by applying the bilinear channel model in Eq. (9) to the optimization program in Eq. (12). Our proposed multichannel blind deconvolution technique is to solve the optimization problem,

minimizeu,aMyΨ̃[(a1u1),,(aNuN)]22,subjecttoΨ̃[(a1u1),,(aNuN)]2=1,am0,m=1,,M,
(13)

where Φ̃ is defined in Eq. (10). If, indeed, our ensemble of channel responses follows (or comes close to following) this model, restricting the set of channel responses that we search over can increase the robustness of the estimate.

The optimization problem stated in Eq. (13) can no longer be solved using a standard eigenvalue decomposition problem. It is also non-convex and so falls outside of the purview of most off-the-shelf solvers. But the problem structure that [(a1u1),,(aNuN)] is block rank-1 allows a natural approach that is a modification of the power method for finding the smallest eigenvector of a matrix. The power method for finding the smallest eigenvector of a (symmetric positive semi-definite) matrix M simply applies (λIM) to an initial vector repeatedly, renormalizing after each application. Our truncated power method is the same with an additional step of enforcing the block rank-1 structure. Our iterative truncated power method is summarized in Algorithm 1 and is not universally guaranteed to solve Eq. (13), but under certain conditions on the channel bases specified by Ψ̃, it is provably effective (Lee et al., 2018).

The matrices in the algorithm do get very large, but everything is structured in such a way that it scales well. Algorithm 1 is iterative with the cost of each iteration dominated by two specific computations: the application of the matrix B and the block rank-1 approximation (the norm of the matrix B in step 3 refers to the nuclear norm, that is, the L1-norm of the vector of its eigenvalues). More specifically, the matrix B is fixed and has dimensions NMW × NMW, where N is the number of paths, M is the number of channels, and W is the number of coefficients used in the joint channel representation. The number of paths N is, in general, pretty small for most ocean waveguides, especially in shallow water where bottom attenuation typically limits the number of observable multipath arrivals, especially for long range propagation. The matrix B, however, never has to be formed—nor stored—explicitly as it consists of an application of Ψ̃ (which generates the channels from their coefficients) and My (which convolves all of the ym with all of the different channels) and their adjoints. The cost of applying Ψ̃ (and Ψ̃T) is O(NMKW), while the cost of applying My (and MyT) is M2/2 length-L convolutions, O(M2LlogL). As K < L, we have

costofapplyingB=O(LM(NW+MlogL)).

For large arrays and observation times, the MlogL term above will dominate, so the cost of applying B is roughly the same as M2 convolutions of length L. This is similar to every blind deconvolution technique using spectral methods, e.g., as used in Xu et al. (1995). Furthermore, the block rank-1 approximation also scales. The optimization variable X is a length NMW vector; this step corresponds to taking this vector, reorganizing it as N different M × W matrices, and then finding the leading singular vectors in each of these matrices. Computing these exactly requires on the order of max(M3,W3) operations (and they can oftentimes be approximated to sufficient accuracy using the power method, which might require only a modest number of applications of the M × W matrix). For large arrays, we expect NW < M. It is also generally true that the number of samples taken per channel is more than the number of channels, L > M. Thus, the cost of the block rank-1 approximation is typically less than applying B.

Algorithm 1: Iterative truncated power method.

input:MyMy,Ψ̃, M, L, û1,…, ûN,â1, …, âN 
output:ĥ 
1 X0[(â1û1),,(âNûN)]
2 BΨ̃(MyMy)Ψ̃
3 λ||B||*
4 whilestop condition not satisfieddo 
5  XλX0BX0
6  û1,…, ûN,â1, …, âN ← BlockRank1Approx(X); 
7  X0[(â1û1),,(âNûN)]
8  X0X0/||X0||2
9 end 
10 ĥΨ̃[(a1u1),,(aNuN)]
input:MyMy,Ψ̃, M, L, û1,…, ûN,â1, …, âN 
output:ĥ 
1 X0[(â1û1),,(âNûN)]
2 BΨ̃(MyMy)Ψ̃
3 λ||B||*
4 whilestop condition not satisfieddo 
5  XλX0BX0
6  û1,…, ûN,â1, …, âN ← BlockRank1Approx(X); 
7  X0[(â1û1),,(âNûN)]
8  X0X0/||X0||2
9 end 
10 ĥΨ̃[(a1u1),,(aNuN)]

Furthermore, designing a good initialization scheme is critical so that the iterates will not get stuck in a local minimum. We propose a variant of the so-called “spectral” method. Let x denote the common source signal provided to the unknown multichannel system. Each output ỹ1,,ỹM can (in the noiseless case) be written as a linear operator applied to the tensor uxaD×L×M. With ỹ=[ỹ1,,ỹM] as the set of concatenated channel outputs and P as the length of ỹ, there are tensors (M1,M2,,MP) with the same dimensions as uxa such that

ỹ[l]=Ml,uxa=i1,i2,i3ml*[i1,i2,i3]*u[i1]x[i2]a[i3],

where ml[i1,i2,i3] are the entries of Ml. The action of these linear functionals can be collected into a single linear operator A:D×L×MP. The adjoint of A applied to a measurement vector y=[y1,,yM] is the tensor given by

A*y=ly[l]MlD×L×M.

Note that the applications of both A and its adjoint A* involve only a set of convolutions and can be implemented by using a fast algorithm. The standard spectral method computes the best rank-1 tensor approximation of A*y. However, there is no guaranteed algorithm that computes the best rank-1 approximation of a three-way tensor exactly. Thus, we instead propose to estimate just one factor u through the following procedure.

Let

Z=A*y×31M,
(14)

where ×3 denotes the tensor multiplication along the third mode. In other words, we fold the tensor A*yD×L×M by summing M slices of size D × L along the third mode with the same weight 1. The choice of the equal unit weight on each slice is justified by the fact that the vector aM of unknown channel gains has positive entries. In the noiseless case, we compute an estimate û of u as the first dominant eigenvector of ZZ. By the construction of A, it follows from Eq. (14) that

Z=A*A(uxa)×31M+A*w,

where w denotes the additive noise to the measurements. Intuitively, if A*A is approximated as the identity operator and there is no noise, then Z is proportional to ux=ux. This is why one can estimate u as the dominant eigenvector of ZZ.

Once an estimate of û is obtained, the estimation of a reduces to the conventional multichannel blind deconvolution with a linear model. Thus, it is computed by an ordinary eigenvalue decomposition. The resulting initialization algorithm is summarized in Algorithm 2.

Algorithm 2: Initialization.

input:A*, y, Ψ̃, M, L 
output:û,â 
1 Z ← folding A*y into a matrix by Eq. (14)
2 QZZ
3 û1,,ûN ← FindMaxEigenvector(Q); 
4 Vû[IMû1,,IMûN]
5 ΞVûΨ̃(MyMy)Ψ̃Vû
6 â1,,âN ← FindMinEigenvector(Ξ); 
input:A*, y, Ψ̃, M, L 
output:û,â 
1 Z ← folding A*y into a matrix by Eq. (14)
2 QZZ
3 û1,,ûN ← FindMaxEigenvector(Q); 
4 Vû[IMû1,,IMûN]
5 ΞVûΨ̃(MyMy)Ψ̃Vû
6 â1,,âN ← FindMinEigenvector(Ξ); 

In this section, we implement the truncated power iteration method presented in Sec. III to solve the multichannel blind deconvolution problem using the bilinear channel model introduced in Sec. II. The proposed approach is demonstrated using at-sea data from a shipping source of opportunity: the Anna Maersk—a passing container ship. Data were recorded on a 32-element VLA moored in between the north- and southbound shipping lanes of the Santa Barbara shipping channel (depth 580 m) during an experiment performed mid-September 2016—see Fig. 4(a). The selected VLA consisted of 2 16-element sub-arrays, one of large aperture (56 m) and one of small aperture (15 m) with uniform 3.75 m and 1 m element spacing, respectively. A representative SSP averaged from multiple CTD casts, which measure conductivity, temperature, and depth, is shown in Fig. 4(b). Data processing and signal analysis were performed in the frequency band from 150 Hz to 3 kHz.

FIG. 4.

(Color online) (a) Illustration of the Santa Barbara Channel Experiment. Four 32-elements are deployed in the Santa Barbara shipping channel. In our simulation, we used data collected from VLA3 from the source of opportunity, which is the passing container ship Anna Maersk represented by the blue line. (b) A representative SSP averaged from multiple measurements at the experiment site along with the bottom geoacoustic parameters used for numerical simulations. (c) Geometry of the bottom-mounted VLAs deployed in the SBC experiment.

FIG. 4.

(Color online) (a) Illustration of the Santa Barbara Channel Experiment. Four 32-elements are deployed in the Santa Barbara shipping channel. In our simulation, we used data collected from VLA3 from the source of opportunity, which is the passing container ship Anna Maersk represented by the blue line. (b) A representative SSP averaged from multiple measurements at the experiment site along with the bottom geoacoustic parameters used for numerical simulations. (c) Geometry of the bottom-mounted VLAs deployed in the SBC experiment.

Close modal

Two different approaches are used to generate the a priori arrival-time information of the CIR used to construct the subspace basis Φ in Eq. (4). For the first data-derived approach, the a priori information for the arrival-times and amplitudes values for the direct path of the CIR were directly obtained using a library of 500 estimated CIRs independently obtained from the RBD method (Byun et al., 2017) since the true CIRs between the Anna Maersk location and the VLA elements were not available otherwise. An example of such estimated CIR from RBD is displayed in Fig. 5, which shows that the first direct arrival is the only dominant wavefront while the subsequent arrival associated with the bottom bounce is much weaker (35dB below the direct arrival's amplitude). Hence, only the direct-arrival is used hereafter to construct and parametrize the CIR library. Specifically, the library of 500 estimated CIRs corresponds to ship locations linearly spaced throughout the course of one minute (12:15:18–12:16:18 on September 16) and covering ranges from 1.71 to 1.94 km as the vessel Anna Maersk cruises away from the VLA at a nearly constant speed of 10 m/s. Speed and position information of the Anna Maersk were obtained using data from the AIS. The superimposed arrival-times profiles for the direct-arrival all 500 estimated CIRs are shown in Fig. 6(a), where the delay-time for the first receiver was arbitrarily set to 20 ms. The relative nature of the CIRs is a consequence of the RBD method because the absolute arrival time of the ray cannot be recovered (Byun et al., 2017).

FIG. 5.

(Color online) Envelope of a relative CIR estimated via the RBD method (Sabra et al., 2010). At the time of the snapshot, the source of opportunity (Anna Maersk) was nearly 1.8 km from the VLA. Note that the first arrival on the top-most hydrophone element has been arbitrarily aligned with t = 20 ms.

FIG. 5.

(Color online) Envelope of a relative CIR estimated via the RBD method (Sabra et al., 2010). At the time of the snapshot, the source of opportunity (Anna Maersk) was nearly 1.8 km from the VLA. Note that the first arrival on the top-most hydrophone element has been arbitrarily aligned with t = 20 ms.

Close modal
FIG. 6.

(Color online) Recovery of M = 31 CIRs using the bilinear channel model (one channel discarded from the VLA because of a hardware problem). (a) 500 snapshots of arrival-times for CIRs when the ship is cruising away from VLA3 at range 1.71 km–1.94 km. (b) The estimated CIRs (red) plot on top of the original CIRs, and the recovery SNR is 25 dB.

FIG. 6.

(Color online) Recovery of M = 31 CIRs using the bilinear channel model (one channel discarded from the VLA because of a hardware problem). (a) 500 snapshots of arrival-times for CIRs when the ship is cruising away from VLA3 at range 1.71 km–1.94 km. (b) The estimated CIRs (red) plot on top of the original CIRs, and the recovery SNR is 25 dB.

Close modal

First, using the method described in Sec. II, a subspace basis matrix Φ—reduced to a 42-dimensional space here as the first 42 singular values of the PCA represented over 99% of the total signal energy— was constructed using the PCA of the data library containing the 500 arrival-times profiles [shown superimposed in Fig. 6(a)], which were obtained from the a priori CIRs computed with the RBD method (Byun et al., 2017). For one arbitrarily chosen ship location (located near the middle of the section of the Anna Maersk vessel's track shown in Fig. 4(a) but not part of the 500 locations used to construct the aforementioned CIR data library), Fig. 6(b) compares the a priori CIR hk (obtained from the ray-based deconvolution) to the estimated CIR hkEstimation obtained using the multichannel blind deconvolution method (introduced in Secs. II and III) applied independently to the raw shipping noise recordings on the VLA. The percentage error between these two CIR (RBD estimate and multichannel-blind method) is very low (0.3%). This is to be expected as the subspace basis matrix Φ—directly obtained from in situ data—provides a very accurate description of the CIR arrival-time structure over the selected Anna Maersk track section, which, in turn, enables a near-perfect recovery of the CIR using the proposed multichannel blind deconvolution method. It should be noted that since no ground truth was available for the CIR (as no active source was used here), the very good agreement between the two independently obtained CIRs shown in Fig. 6(b) can only be interpreted as a simple experimental validation of the proposed multichannel blind deconvolution method but does not serve as an exact error analysis.

On the other hand, building the subspace basis matrix Φ directly from data-derived estimates of the CIRs may not always be available. Instead, now we use a model-based approach as the second approach to generate the a priori arrival-time information of the CIR. Since the VLA is in the far-field of the Anna Maersk for the selected track section, we can well model the direct-path arrival as a straight line using the ray-approximation of the CIR presented in Sec. II A [see Eq. (2)]. Under this assumption, the estimated arrival-time uncertainty for each ray-path of the CIRs is obtained by simulating the direct arrival-time and amplitude using the software BELLHOP (Porter, 2011) for the same 500 source locations using the environmental parameters shown in Fig. 4(b) and assuming a source depth of 5 m. The simulated direct arrival profile computed from Bellhop are shown superimposed in Fig. 7(a). For this simple linear arrival model, all 500 simulated direct-path arrivals are highly correlated; hence, the associated subspace matrix Φ for this model-based approach has a much smaller dimension of 12 than does the subspace matrix obtained with the data-derived approach (whose dimension was 42). However, this model-based subspace matrix Φ is less accurate than the data-derived one as it does not actually account for all experimental variations and actual deviations of the direct path wavefront from a simple straight line, which explains the relatively larger percentage error (12%) shown Fig. 7(b) for this model-based approach between the CIR obtained from the ray-based deconvolution to the estimated CIR obtained using the multichannel blind deconvolution method for one arbitrarily selected ship location (located near the middle of the section of the Anna Maersk vessel's track shown in Fig. 4(a) but not part of the 500 locations used to construct the aforementioned CIR model-based library). Hence, these results demonstrate that the model-based approach can provide acceptable deconvolution estimates of the CIR for surface sources of opportunity—such as a shipping vessel—provided that the environmental and geoacoustic parameters used as input for the numerical simulations are sufficiently accurate (which is a classical and known limitation for all types of model-based approaches).

FIG. 7.

(Color online) Recovery of M = 31 CIRs using the bilinear channel model (one channel discarded from the VLA because of a hardware problem). (a) Linear approximations of arrival-times for CIRs when the ship is 1.8 km away from VLA3, this approximation accounts for the predicted environmental variability affecting the CIR at the test site for the track section of the Anna Maersk vessel shown in Fig. 4(a). (b) The estimated CIRs (red) plot on top of the original CIRs, and the recovery SNR is 12 dB.

FIG. 7.

(Color online) Recovery of M = 31 CIRs using the bilinear channel model (one channel discarded from the VLA because of a hardware problem). (a) Linear approximations of arrival-times for CIRs when the ship is 1.8 km away from VLA3, this approximation accounts for the predicted environmental variability affecting the CIR at the test site for the track section of the Anna Maersk vessel shown in Fig. 4(a). (b) The estimated CIRs (red) plot on top of the original CIRs, and the recovery SNR is 12 dB.

Close modal

The results from Figs. 6(b) and 7(b) indicate that the primary source of error for the multichannel blind deconvolution method is the inaccuracy of the subspace basis matrix Φ such that it cannot accurately represent the actual arrival-times of the in situ CIR for the location of interest. In order to further quantify this effect, numerical simulations were conducted using the same model-based subspace basis matrix Φ, which was obtained from the arrival-time profiles shown in Fig. 7(a), and applying the multichannel blind deconvolution method to deconvolve a simulated CIR hk whose arrival-times deviate from the straight line approximation by adding an increasing amount of random fluctuations [expressed as an off-line coefficient, indicating the percentage of fluctuations of the linear slope angle compared to the reference CIR obtained using the environmental parameters shown in Fig. 4(b)]. Such arrival-time fluctuations can be representative of, for instance, array element position uncertainties or localized sound-speed fluctuations unaccounted for by the reference SSP shown in Fig. 4(b). The simulated received signals on the VLA were synthesized by convolving the simulated CIR hk with a Gaussian white noise signal filtered the same bandwidth (150–3000 Hz) representative of shipping noise spectra. Despite the noise being predominately low frequency (<1 kHz), using a large frequency band improves the SNR of the RBD estimated CIR (Durofchalk et al., 2019). A Monte Carlo simulation with 100 trials for each noise level was performed.

Two types of relative errors can be defined and expressed in terms of the SNR (i.e., by taking the decimal logarithm of the inverse of the relative error). First, the model approximation SNR is defined as 10log10(||hk||2/||hkhk,approx||2), where hk,approx is the approximation of CIRs using the subspace basis Φ. This first SNR quantifies how well the subspace basis Φ, constructed from linear arrival-time profiles, describes the actual CIR hk, which has the random variations in arrival-times. Second, the CIR estimation SNR is defined as 10log10(||hk||2/||hkhk,est||2), where hk,est is the estimated CIR using the multichannel blind deconvolution method. This second SNR quantifies how accurate the CIR recovery is. As expected, both SNR values decrease (i.e., the corresponding error increases) for increasing percentage of the arrival-times off-line error (see Fig. 8).

FIG. 8.

(Color online) Illustration of the blind deconvolution algorithm robustness performance when arrival-times are not exactly on a line, whereas the bilinear model is formed with the assumption that all arrival-times fall along a line. The red dots indicate the approximation errors when using such a channel model, and the blue dots are the estimation errors using such a channel model to perform blind deconvolution. The model errors are measured by an off-line coefficient, which indicates the ratio between the off-line time and the arrival-time differences among channels.

FIG. 8.

(Color online) Illustration of the blind deconvolution algorithm robustness performance when arrival-times are not exactly on a line, whereas the bilinear model is formed with the assumption that all arrival-times fall along a line. The red dots indicate the approximation errors when using such a channel model, and the blue dots are the estimation errors using such a channel model to perform blind deconvolution. The model errors are measured by an off-line coefficient, which indicates the ratio between the off-line time and the arrival-time differences among channels.

Close modal

Here, for the selected simulation parameters and source-receiver configuration, as long as the arrival-time fluctuations of the actual CIR hk remains below 7%, the model approximation SNR remains approximately over 9 dB (i.e., the actual CIR hk fits the subspace approximation with less than 12% relative error) and the CIR estimation SNR remains approximately over 7 dB (i.e., the actual CIR hk fits the subspace approximation with less than 20% relative error). Furthermore, the multichannel blind deconvolution method is also able to cope with larger arrival-times variations in the library CIRs (e.g., for the case of a longer section of the shipping track), but this increases the dimension of the subspace matrix Φ basis and, thus, makes the recovery less computationally tractable. For instance, when doubling the spread of the arrival-angle of the direct path [see Fig. 9(a)], Fig. 9(b) shows that both the model approximation SNR and CIR estimation SNR values get worse than the results in Fig. 8. In this case, as long as the arrival-time fluctuations of the actual CIR hk remains below 7%, the model approximation SNR remains approximately over 8.5 dB (i.e., the actual CIR hk fits the subspace approximation with less than 14% relative error) and the CIR estimation SNR remains approximately over 4 dB (i.e., the actual CIR hk fits the subspace approximation with less than 40% relative error).

FIG. 9.

(Color online) (a) A subspace with a larger dimension is used in this simulation, which is generated from a wider range of arrival time fluctuations when compared to Fig. 7(a). (b) The red dots indicate the approximation errors when using such a channel model, and the blue dots are the estimation errors using such a channel model to perform blind deconvolution. The model errors are measured by an off-line coefficient, which indicates the ratio between the off-line time and the arrival-time differences among channels.

FIG. 9.

(Color online) (a) A subspace with a larger dimension is used in this simulation, which is generated from a wider range of arrival time fluctuations when compared to Fig. 7(a). (b) The red dots indicate the approximation errors when using such a channel model, and the blue dots are the estimation errors using such a channel model to perform blind deconvolution. The model errors are measured by an off-line coefficient, which indicates the ratio between the off-line time and the arrival-time differences among channels.

Close modal

Finally, as discussed at the beginning of Sec. III, classical cross-convolution methods for multichannel blind deconvolution (which are based on the commutativity of the convolution) are very sensitive to the noise level present in the measurement (referred to hereafter as measurement noise). To address this limitation, Sec. III introduced the truncated power iteration method under the bilinear channel model to solve the multichannel blind deconvolution problem. Using the same parameters and subspace basis Φ constructed to generate the results from Fig. 6(b), Fig. 10 demonstrates the robustness of the proposed approach. The measurement SNR (which is used to quantify the noise level along the horizontal axis) is defined as the decimal logarithm of the ratio of the amplitude of the boat signal to the variable amplitude of the additive white Gaussian noise used to synthesize the noisy measurements. For each measurement SNR value, 100 Monte Carlo realizations of noisy measurements were synthesized—assuming that the noise waveforms are uncorrelated among each receiver—and the corresponding averaged value of the CIR estimation SNR, obtained using the truncated power iteration method, is displayed along the vertical axis. Figure 10 demonstrates that the CIR estimation SNR smoothly decays as the measurement SNR increases, which confirms the stability of the truncated power iteration method for solving the multichannel blind deconvolution problem in the presence of an increasing level of additive noise in the original measurements.

FIG. 10.

(Color online) Evolution of the SNR of the estimated CIR from the multichannel blind deconvolution vs the measurement SNR, which quantifies the amount of uncorrelated white noise added to the raw receiver data used for deconvolution.

FIG. 10.

(Color online) Evolution of the SNR of the estimated CIR from the multichannel blind deconvolution vs the measurement SNR, which quantifies the amount of uncorrelated white noise added to the raw receiver data used for deconvolution.

Close modal

This study introduced a multichannel blind deconvolution algorithm that imposes a bilinear model on the arrival structure of the CIRs and demonstrated its robustness using experimental recordings of a shipping source of opportunity on a bottom-mounted short VLAs in an ocean waveguide. Such a bilinear model was obtained, for example, by embedding a parametric model for the arrival-time structure of the CIR wavefronts—i.e., jointly across receivers—into a low-dimensional subspace using PCA while the receiver-dependent arrival amplitudes are treated as independent variables. Furthermore, under the bilinear channel model, the truncated power iteration method was developed as a modification to the cross-convolution method for blind deconvolution algorithms introduced by previous studies in order to mitigate their high sensitivity to measurement noise. The bilinear system model imposes a strong prior on the arrival structure of the CIR, which enables us to efficiently narrow the search space using the PCA and ultimately estimate the unknown CIR with short observation times (i.e., using short snapshot duration), which can be advantageous for the case of a moving source. Hence, it is expected that this multichannel blind deconvolution method will perform well when sufficient a priori information of the CIR arrival structure is available; and, thus, this method may potentially be relevant for passive ocean acoustic tomography using sources of opportunity as tomography applications typically rely on accurately knowing the CIR (i.e., a known reference set of environmental parameters) to invert for first-order fluctuations of the environmental parameters from in situ observations.

This work was sponsored by the Office of Naval Research (Code 321) under Contract No. N00014-17-1-2189. The authors also wish to acknowledge the Marine Physical Laboratory, Scripps Institute of Oceanography at the University of California (UC) San Diego for making available the SBCEx-16 experimental data set used in the study. The work by K.L. and J.R. has been supported in part by the National Science Foundation (NSF CCF) 17-18771.

1.
Abadi
,
S. H.
,
Rouseff
,
D.
, and
Dowling
,
D. R.
(
2012
). “
Blind deconvolution for robust signal estimation and approximate source localization
,”
J. Acoust. Soc. Am.
131
(
4
),
2599
2610
.
2.
Broadhead
,
M. K.
, and
Pflug
,
L. A.
(
2000
). “
Performance of some sparseness criterion blind deconvolution methods in the presence of noise
,”
J. Acoust. Soc. Am.
107
(
2
),
885
893
.
3.
Byun
,
S.-H.
,
Verlinden
,
C. M.
, and
Sabra
,
K. G.
(
2017
). “
Blind deconvolution of shipping sources in an ocean waveguide
,”
J. Acoust. Soc. Am.
141
(
2
),
797
807
.
4.
Cazzolato
,
B. S.
,
Nelson
,
P.
,
Joseph
,
P.
, and
Brind
,
R. J.
(
2001
). “
Numerical simulation of optimal deconvolution in a shallow-water environment
,”
J. Acoust. Soc. Am.
110
(
1
),
170
185
.
5.
Choudhary
,
S.
, and
Mitra
,
U.
(
2012
). “
Sparse recovery from convolved output in underwater acoustic relay networks
,” in
Proceedings of the 2012 Asia Pacific Signal and Information Processing Association Annual Summit and Conference
, pp.
1
8
.
6.
Choudhary
,
S.
, and
Mitra
,
U.
(
2014
). “
Sparse blind deconvolution: What cannot be done
,” in
2014 IEEE International Symposium on Information Theory
, pp.
3002
3006
.
7.
Durofchalk
,
N. C.
,
Zhang
,
X.
, and
Sabra
,
K. G.
(
2019
). “
Analysis of a ray-based blind deconvolution algorithm on ships of opportunity in the Santa Barbara channel
,”
J. Acoust. Soc. Am.
145
,
1935
1935
.
8.
Giannakis
,
G. B.
, and
Serpedin
,
E.
(
1997
). “
Linear multichannel blind equalizers of nonlinear fir volterra channels
,”
IEEE Trans. Signal Process.
45
(
1
),
67
81
.
9.
Godin
,
O. A.
(
2007
). “
Restless rays, steady wave fronts
,”
J. Acoust. Soc. Am.
122
(
6
),
3353
3363
.
10.
Gürelli
,
M.
, and
Nikias
,
C. L.
(
1995
). “
Evam: An eigenvector-based algorithm for multichannel blind deconvolution of input colored signals
,”
IEEE Trans. Signal Process.
43
(
1
),
134
149
.
11.
Jensen
,
F. B.
,
Kuperman
,
W. A.
,
Porter
,
M. B.
, and
Schmidt
,
H.
(
2011
).
Computational Ocean Acoustics
(
Springer Science and Business Media
,
New York
).
12.
Jolliffe
,
I.
(
2011
).
Principal Component Analysis
(
Springer
,
Berlin
), pp.
1094
1096
.
13.
Katsaggelos
,
A. K.
(
1991
).
Digital Image Restoration
(
Springer
,
Berlin
).
14.
Kazemi
,
N.
, and
Sacchi
,
M. D.
(
2014
). “
Sparse multichannel blind deconvolution
,”
Geophysics
79
(
5
),
V143
V152
.
15.
Lee
,
K.
,
Tian
,
N.
, and
Romberg
,
J.
(
2018
). “
Fast and guaranteed blind multichannel deconvolution under a bilinear system model
,”
IEEE Trans. Inform. Theory
64
(
7
),
4792
4818
.
16.
Li
,
Y.
,
Lee
,
K.
, and
Bresler
,
Y.
(
2016
). “
Identifiability in blind deconvolution with subspace or sparsity constraints
,”
IEEE Trans. Inform. Theory
62
(
7
),
4266
4275
.
17.
Mendel
,
J. M.
(
2013
).
Optimal Seismic Deconvolution: An Estimation-Based Approach
(
Elsevier
,
New York
).
18.
Morrison
,
R. L.
,
Jacob
,
M.
, and
Do
,
M. N.
(
2007
). “
Multichannel estimation of coil sensitivities in parallel MRI
,” in
2007 4th IEEE International Symposium on Biomedical Imaging: From Nano to Macro
, pp.
117
120
.
19.
Moulines
,
E.
,
Duhamel
,
P.
,
Cardoso
,
J.-F.
, and
Mayrargue
,
S.
(
1995
). “
Subspace methods for the blind identification of multichannel fir filters
,”
IEEE Trans. Signal Process.
43
(
2
),
516
525
.
20.
Porter
,
Michael
. (
2011
). The bellhop manual and user's guide: Preliminary draft; see http://oalib.hlsresearch.com/Rays/.
21.
Roan
,
M. J.
,
Gramann
,
M. R.
,
Erling
,
J. G.
, and
Sibul
,
L. H.
(
2003
). “
Blind deconvolution applied to acoustical systems identification with supporting experimental results
,”
J. Acoust. Soc. Am.
114
(
4
),
1988
1996
.
22.
Roux
,
P.
,
Cornuelle
,
B. D.
,
Kuperman
,
W. A.
, and
Hodgkiss
,
W. S.
(
2008
). “
The structure of raylike arrivals in a shallow-water waveguide
,”
J. Acoust. Soc. Am.
124
(
6
),
3430
3439
.
23.
Sabra
,
K.
, and
Dowling
,
D.
(
2004
). “
Blind deconvolution in oceanic waveguides using artificial time reversal
,”
J. Acoust. Soc. Am.
116
,
262
271
.
24.
Sabra
,
K. G.
,
Song
,
H.-C.
, and
Dowling
,
D. R.
(
2010
). “
Ray-based blind deconvolution in ocean sound channels
,”
J. Acoust. Soc. Am.
127
(
2
),
EL42
EL47
.
25.
She
,
H.
,
Chen
,
R.-R.
,
Liang
,
D.
,
Chang
,
Y.
, and
Ying
,
L.
(
2015
). “
Image reconstruction from phased-array data based on multichannel blind deconvolution
,”
Magn. Reson. Imaging
33
(
9
),
1106
1113
.
26.
Smith
,
G. B.
(
2000
). “
Blind deconvolution for multipath mitigation in shallow water acoustics
,”
J. Acoust. Soc. Am.
107
(
5
),
2868
.
27.
Souidene
,
W.
,
Abed-Meraim
,
K.
, and
Beghdadi
,
A.
(
2009
). “
A new look to multichannel blind image deconvolution
,”
IEEE Trans. Image Process.
18
(
7
),
1487
1500
.
28.
Sroubek
,
F.
,
Cristóbal
,
G.
, and
Flusser
,
J.
(
2007
). “
A unified approach to superresolution and multichannel blind deconvolution
,”
IEEE Trans. Image Process.
16
(
9
),
2322
2332
.
29.
Sroubek
,
F.
, and
Milanfar
,
P.
(
2011
). “
Robust multichannel blind deconvolution via fast alternating minimization
,”
IEEE Trans. Image Process.
21
(
4
),
1687
1700
.
30.
Tian
,
N.
,
Byun
,
S.-H.
,
Sabra
,
K.
, and
Romberg
,
J.
(
2017
). “
Multichannel myopic deconvolution in underwater acoustic channels via low-rank recovery
,”
J. Acoust. Soc. Am.
141
(
5
),
3337
3348
.
31.
Verlinden
,
C. M.
,
Sarkar
,
J.
,
Hodgkiss
,
W.
,
Kuperman
,
W.
, and
Sabra
,
K.
(
2015
). “
Passive acoustic source localization using sources of opportunity
,”
J. Acoust. Soc. Am.
138
(
1
),
EL54
EL59
.
32.
Xu
,
G.
,
Liu
,
H.
,
Tong
,
L.
, and
Kailath
,
T.
(
1995
). “
A least-squares approach to blind channel identification
,”
IEEE Trans. Signal Process.
43
(
12
),
2982
2993
.