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.

## I. INTRODUCTION

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.

## II. THEORY

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.

### A. An example of the bilinear channel model

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

where *p*(*t*) denotes the baseline pulse waveform [Fig. 1(a)] that is determined by the frequency response of the receivers, and *t _{m}* [Fig. 1(c)] and

*a*[Fig. 1(b)] are, respectively, the values of the arrival-time and amplitude of the ray arrival for channel

_{m}*m*. While the true channel response is an impulse $am\delta (t\u2212tm)$ (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).

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 *t _{m}* 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

where *t*_{1} 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 *a _{m}* and one variable vector $t=[tm]$ to represent the linked ray arrival-times

*t*. Furthermore, we also discretize each time-shifted pulse profile $p(t\u2212tm)$ [see Eq. (1), $m=1,\u2026,M$] recorded on the

_{m}*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\xaf$,

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 *a _{m}* and the pulse shape $pm$. [e.g., given by Eq. (1) for the ray-model].

### B. General formulation for bilinear channel models

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\xaf$ [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\xafi$. 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\xafi$ 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 $\Phi \u2208\mathbb{R}MK\xd7W$ (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\xaf$, we can write it as

where $u\u2208\mathbb{R}W$ is the coefficient vector that fully describes the corresponding realization of arrival-times ** t**.

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.

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

where $(\Phi m)m=1M\u2208\mathbb{R}K\xd7W$ are the blocks of the subspace matrix $\Phi $ defined in Eq. (4), $u\u2208\mathbb{R}W$ 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,\u2026,aM]\u22a4\u2208\mathbb{R}M$ with $am>0$ is the arrival amplitude *a _{m}*. Furthermore, let $h\xaf\u2208\mathbb{R}MK$ denote the concatenation of the CIR coefficients for all channels, i.e.,

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

where

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

where

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.

## III. TRUNCATED POWER ITERATION ALGORITHM FOR MULTICHANNEL BLIND DECONVOLUTION

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=s\u2217hi$ and, therefore, for any $j\u2260i$,

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

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

where $My$ is the block-sparse matrix

In other words, we know that when there is no noise, $h\xaf$ 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\xaf$ 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

which, in turn, is equivalent to actually finding the eigenvector corresponding to the minimum eigenvalue of $My\u22a4My$. 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 $My\u22a4My$ 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,

where $\Phi \u0303$ 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 $[(a1\u2297u1)\u22a4,\u2026,(aN\u2297uN)\u22a4]\u22a4$ 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 $(\lambda I\u2212M)$ 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 $\Psi \u0303$, 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

**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**

*B**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

**, however, never has to be formed—nor stored—explicitly as it consists of an application of $\Psi \u0303$ (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 $\Psi \u0303$ (and $\Psi \u0303T$) is**

*B**O*(

*NMKW*), while the cost of applying $My$ (and $MyT$) is $M2/2$ length-

*L*convolutions, $O(M2L\u2009\u2009log\u2009\u2009L)$. As

*K*<

*L*, we have

For large arrays and observation times, the $M\u2009\u2009log\u2009\u2009L$ term above will dominate, so the cost of applying ** B** is roughly the same as

*M*

^{2}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

**is a length**

*X**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* input: $My\u22a4My,\u2009\Psi \u0303$, M, L, $u\u03021$,…, $u\u0302N,\u2009a\u03021$, …, $a\u0302N$ |

output: $h\u0302$ |

1 $X0$ ← $[(a\u03021\u2297u\u03021)\u22a4,\u2026,(a\u0302N\u2297u\u0302N)\u22a4]\u22a4$; |

2 ← $\Psi \u0303\u22a4(My\u22a4My)\Psi \u0303$; B |

3 λ ← $||B||*$; |

4 while stop condition not satisfied do |

5 ← $\lambda X0\u2212BX0$; X |

6 $u\u03021$,…, $u\u0302N,\u2009a\u03021$, …, $a\u0302N$ ← BlockRank1Approx(); X |

7 $X0$ ← $[(a\u03021\u2297u\u03021)\u22a4,\u2026,(a\u0302N\u2297u\u0302N)\u22a4]\u22a4$; |

8 $X0$ ← $X0/||X0||2$, |

9 end |

10 $h\u0302$ ← $\Psi \u0303[(a1\u2297u1)\u22a4,\u2026,(aN\u2297uN)\u22a4]\u22a4$; |

input: $My\u22a4My,\u2009\Psi \u0303$, M, L, $u\u03021$,…, $u\u0302N,\u2009a\u03021$, …, $a\u0302N$ |

output: $h\u0302$ |

1 $X0$ ← $[(a\u03021\u2297u\u03021)\u22a4,\u2026,(a\u0302N\u2297u\u0302N)\u22a4]\u22a4$; |

2 ← $\Psi \u0303\u22a4(My\u22a4My)\Psi \u0303$; B |

3 λ ← $||B||*$; |

4 while stop condition not satisfied do |

5 ← $\lambda X0\u2212BX0$; X |

6 $u\u03021$,…, $u\u0302N,\u2009a\u03021$, …, $a\u0302N$ ← BlockRank1Approx(); X |

7 $X0$ ← $[(a\u03021\u2297u\u03021)\u22a4,\u2026,(a\u0302N\u2297u\u0302N)\u22a4]\u22a4$; |

8 $X0$ ← $X0/||X0||2$, |

9 end |

10 $h\u0302$ ← $\Psi \u0303[(a1\u2297u1)\u22a4,\u2026,(aN\u2297uN)\u22a4]\u22a4$; |

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 $y\u03031,\u2026,y\u0303M$ can (in the noiseless case) be written as a linear operator applied to the tensor $u\u2297x\u2297a\u2208\mathbb{R}D\xd7L\xd7M$. With $y\u0303=[y\u03031\u22a4,\u2026,y\u0303M\u22a4]\u22a4$ as the set of concatenated channel outputs and

*P*as the length of $y\u0303$, there are tensors $(M1,M2,\u2026,MP)$ with the same dimensions as $u\u2297x\u2297a$ such that

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:\mathbb{R}D\xd7L\xd7M\u2192\mathbb{R}P$. The adjoint of $A$ applied to a measurement vector $y=[y1\u22a4,\u2026,yM\u22a4]\u22a4$ is the tensor given by

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

where $\xd73$ denotes the tensor multiplication along the third mode. In other words, we fold the tensor $A*y\u2208\mathbb{R}D\xd7L\xd7M$ 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 $a\u2208\mathbb{R}M$ of unknown channel gains has positive entries. In the noiseless case, we compute an estimate $u\u0302$ of ** u** as the first dominant eigenvector of $ZZ\u22a4$. By the construction of $A$, it follows from Eq. (14) that

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

**is proportional to $u\u2297x=ux\u22a4$. This is why one can estimate**

*Z***as the dominant eigenvector of $ZZ\u22a4$.**

*u*Once an estimate of $u\u0302$ 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.

input: $A*$, , $\Psi \u0303$, yM, L |

output $u\u0302,\u2009a\u0302$ : |

1 ← folding $A*y$ into a matrix by Eq. (14); Z |

2 ← $ZZ\u22a4$; Q |

3 $u\u03021,\u2009\u2026,\u2009u\u0302N$ ← FindMaxEigenvector(); Q |

4 $Vu\u0302$ ← $[IM\u2297u\u03021,\u2026,IM\u2297u\u0302N]$; |

5 $\Xi $ ← $Vu\u0302\u22a4\Psi \u0303\u22a4(My\u22a4My)\Psi \u0303Vu\u0302$; |

6 $a\u03021,\u2009\u2026,\u2009a\u0302N$ ← FindMinEigenvector($\Xi $); |

input: $A*$, , $\Psi \u0303$, yM, L |

output $u\u0302,\u2009a\u0302$ : |

1 ← folding $A*y$ into a matrix by Eq. (14); Z |

2 ← $ZZ\u22a4$; Q |

3 $u\u03021,\u2009\u2026,\u2009u\u0302N$ ← FindMaxEigenvector(); Q |

4 $Vu\u0302$ ← $[IM\u2297u\u03021,\u2026,IM\u2297u\u0302N]$; |

5 $\Xi $ ← $Vu\u0302\u22a4\Psi \u0303\u22a4(My\u22a4My)\Psi \u0303Vu\u0302$; |

6 $a\u03021,\u2009\u2026,\u2009a\u0302N$ ← FindMinEigenvector($\Xi $); |

## IV. RESULTS

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 $\u2248580$ 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 ($\u224856$ m) and one of small aperture ($\u224815$ 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.

Two different approaches are used to generate the *a priori* arrival-time information of the CIR used to construct the subspace basis $\Phi $ 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 ($\u2248\u200935\u2009\u2009dB$ 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).

First, using the method described in Sec. II, a subspace basis matrix $\Phi $—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 ($\u22480.3%$). This is to be expected as the subspace basis matrix $\Phi $—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 $\Phi $ 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 $\Phi $ 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 $\Phi $ 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 ($\u224812%$) 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).

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 $\Phi $ 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 $\Phi $, 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 $10\u2009\u2009log10(||hk||2/||hk\u2212hk,approx||2)$, where $hk,approx$ is the approximation of CIRs using the subspace basis $\Phi $. This first SNR quantifies how well the subspace basis $\Phi $, 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 $10\u2009\u2009log10(||hk||2/||hk\u2212hk,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).

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 $\Phi $ 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).

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 $\Phi $ 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.

## V. CONCLUSIONS

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.

## ACKNOWLEDGMENT

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.