Estimation of Overspread Underwater Acoustic Channel Based on Low-Rank Matrix Recovery

Abstract

In this paper, the estimation of overspread, i.e., doubly spread underwater acoustic (UWA) channels of strong dispersion is considered. We show that although the UWA channel dispersion causes the degeneration of channel sparsity, it leads to a low-rank structure especially when the channel delay-Doppler-spread function is separable in delay and Doppler domain. Therefore, we introduce the low-rank criterion to estimate the UWA channels, which can help to improve the estimation performance in the case of strong dispersion. The estimator is based on the discrete delay-Doppler-spread function representation of channel, and is formulated as a low-rank matrix recovery problem which can be solved by the singular value projection technique. Simulation examples are carried out to demonstrate the effectiveness of the proposed low-rank-based channel estimator.

1. Introduction

Due to multi-path propagation and time-varying nature, the underwater acoustic (UWA) channel is known to be doubly spread in delay and Doppler domain [1,2,3,4], which is also referred to as doubly selective in the literature [5,6,7]. To combat the effects of delay-Doppler spread in time-varying UWA channels, accurate estimation of the multi-path delay, Doppler frequency, and the channel gain is needed, which is a challenging task for high- speed UWA communications [8,9].

Channel estimation is usually performed with the aid of training signals. One well-known scheme is the least-squares (LS) estimator. However, due to lack of exploiting any prior information of the channel, LS estimators usually require many training signal measurements to achieve good estimation, especially when the number of unknown parameters is large [10]. Therefore, various channel estimation schemes have been proposed by considering the sparsity of channels, in which the technique of compressed sensing (CS) plays an important role [11]. For example, the orthogonal matching pursuit (OMP) and basis pursuit (BP) algorithms are adopted in [12] to estimate sparse doubly spread channel, which are shown to outperform the LS approach and subspace methods. It was further shown that BP has better performance than OMP in the UWA environment at the expense of much higher computational complexity [13]. In [2], a sparse channel estimation technique is developed based on the delay-Doppler-spread function representation of the channel. In [3], a computationally efficient two-stage sparse channel estimation technique is developed by parameterizing the amplitude variation and delay variation of each path with polynomial approximation. In a recent work [14], a sparse channel estimation technique has been proposed based on the Alamouti’s spacetime block coding with transmit diversity scheme in the form of two transmit antennas and one receiver.

The performance of the CS-based channel estimation relies highly on the sparsity of the channel. However, in practice, the sparsity assumption does not always hold for UWA channels due to the leakage effect and rich scattering environment [11,15]. In fact, different from the land-based radio channels which are known to be underspread, the UWA channels can tend to be overspread [16], which means that the UWA channels exhibit a strong dispersion/diffusion along delay and Doppler domain. The performance CS-based estimator can be significantly degraded due to the poor sparsity of overspread UWA channels [11,17,18]. Several techniques have been proposed to deal with the less sparse case. For example, the basis expansion [6,11,17] has been considered to enhance the sparsity of channels. In [7,17,18], some block/group compressed sensing techniques are introduced to fit the block sparsity structure of channels.

In this paper, the estimation of a dispersive overspread UWA channel is considered. By formulating the channel in a discrete matrix with respect to the delay-Doppler-spread function (called the delay-Doppler-spread matrix (DDSM) in this paper), we show that the channel dispersion gives rise to a useful low-rank structure of DDSM, especially when the channel delay-Doppler-spread function is separable in delay and Doppler domain. This low-rank structure reveals that although the overspread channel is not sparse in delay-Doppler domain, the DDSM can be determined by only a small number of singular values. Therefore, we introduce the low-rank criterion to estimate the overspread channels, which can help to improve the estimation performance in case of strong dispersive UWA channels. The estimator is formulated in a low-rank matrix recovery problem and is solved by the singular value projection (SVP) technique. Numerical experiments are carried out to demonstrate the effectiveness of the proposed low-rank-based channel estimator.

The rest of the paper is organized as follows. In Section 2, we give a discrete model of doubly spread channel-based on the delay-Doppler-spread function representation. In Section 3, we introduce the concept of underspread and overspread of wireless channels, and show the degeneration of sparsity in case of large dispersion spread. In Section 4, we show the low-rank structure of UWA channel and describe the proposed method. Finally, some computer simulation results are given to illustrate the behavior of the proposed method in Section 5, and the conclusions are drawn in Section 6.

Notations: Vectors are denoted by bold lower-case letters, while matrices use upper-case bold letters. $\mathbb{C}$, $\mathbb{R}$ denote the set of complex and real numbers, respectively. ${\mathit{A}}^T$, ${\mathit{A}}^H$ and ${\mathit{A}}^{†}$ are the transpose, conjugate transpose and Pseudo inverse of $\mathit{A}$, respectively. $\mathrm{vec}(\mathit{A})$ denotes the vectorization operator that stacks matrix $\mathit{A}$ column by column. $\mathrm{rank}(•)$ are the rank operators, respectively. $\mathbf{0}$ denotes the zero matrix. ${\parallel \mathit{A}\parallel }_2$, ${\parallel \mathit{A}\parallel }_F$ and ${\parallel \mathit{A}\parallel }_{*}$ denote the ${\ell }_2$-norm, Frobenius norm and nuclear norm of $\mathit{A}$, respectively.

2. System Model

Consider signal $x(t)$ propagating over a narrow-band time-varying UWA channel, the received signal $y(t)$ can be written as

$$y(t)=\int h(t,\tau )x(t-\tau )d\tau +\omega (t),$$

(1)

where $h(t,\tau )$ is the channel impulse response and $\omega (t)$ is additive noise. The delay-Doppler spreading function $S_h(\tau ,f)$, which represents the channel in terms of delay and Doppler frequency, is defined as the Fourier transform of $h(t,\tau )$

$$S_h(\tau ,v)=\int h(t,\tau )e^{-j2\pi vt}dt$$

(2)

where v is the variable in the Doppler domain. By combining (1) and (2), we have

$$y(t)=\int \int S_h(\tau ,v)x(t-\tau )e^{j2\pi vt}d\tau dv+w(t)$$

(3)

We assume that the bandwidth and the time duration of the transmitted signal are W and T, respectively, which corresponds to a signal length of $N=WT$. Moreover, we assume that the channel has limited spreads in delay-Doppler domain, i.e., the delay spreads are within $\tau \in [0,{\tau }_{\max }]$, and Doppler spreads are within $v\in [-v_{\max }/2,v_{\max }/2]$. Then, by uniformly sampling the delay-Doppler space with Nyquist sampling rate $(\Delta \tau ,\Delta v)=(1/W,1/T)$, we can obtain a discrete representation of (3) as [19]

$$\begin{array}{c}\hfill y\left[n\right]=\sum _{l=0}^{L-1}\sum _{m=-M}^M{S}_h[l,m]e^{j2\pi mn/N_0}x[n-l]+\omega \left[n\right],\\ \hfill n=0,\dots ,N_0-1.\end{array}$$

(4)

where $L=\lceil W{\tau }_{\max }\rceil +1$, $M=\lceil T{v}_{\max }/2\rceil$ and $N_0=N+L-1$.

Collecting $N_0$ points of the received signal, we can rewrite (4) in a matrix form as

$$\mathit{y}=\mathit{A}\mathrm{vec}\left[{\mathit{S}}_h\right]+\mathit{\omega }$$

(5)

where

$$\mathit{y}={[y\left[0\right],\dots ,y[N_0-1]]}^T$$

$$\mathit{\omega }={[\omega \left[0\right],\dots ,\omega [N_0-1]]}^T$$

and ${\mathit{S}}_h\in {\mathbb{C}}^{L\times 2M+1}$ is the DDSM with its $(l,m)$th element being $S_h[l,m]$.

$\mathit{A}\in {\mathbb{C}}^{N_0\times L(2M+1)}$ in (5) is a block matrix of the form

$$\mathit{A}=[{\mathsf{\Phi }}_{-M}\mathit{X},\dots ,{\mathsf{\Phi }}_M\mathit{X}].$$

where ${\mathsf{\Phi }}_m\in {\mathbb{C}}^{N_0\times N_0}$ is a diagonal matrix given by ${\mathsf{\Phi }}_m=\mathrm{diag}[e^{j2\pi m·0/N_0},\dots ,e^{j2\pi m(N_0-1)/N_0}]$, and $\mathit{X}\in {\mathbb{C}}^{N_0\times L}$ is a Toeplitz matrix whose first row and first column are given by $[x\left[0\right],{\mathbf{0}}_{L-1}^T]$ and ${[x\left[0\right],\dots ,x[N-1],{\mathbf{0}}_{L-1}^T]}^T$, respectively.

Therefore, given measurement vector $\mathit{y}$ and sensing matrix $\mathit{A}$, the doubly spread channel estimation problem has been transformed to the estimation of DDSM ${\mathit{S}}_h$ as shown in (5).

3. Overspread Channels

The doubly spread channel is commonly referred to be underspread if ${\tau }_{\max }{v}_{\max }\le 1$, and overspread otherwise [16]. It is widely known that wireless radio channels are always underspread. This is because the channel delay and Doppler are both inversely proportional to the speed of light [16]. However, situation can be very different for UWA channels, since the speed of sound in water is much slower than the light speed. Specifically, we have ${\tau }_{\max }=d_{\max }/c$ and $v_{\max }={\nu }_{\max }{f}_c/c$, where $d_{\max }$ is the maximum path length, ${\nu }_{\max }$ is the maximum relative velocity, $f_c$ is the carrier frequency and c is the speed of light or sound. If the channel is overspread, we need $d_{\max }{\nu }_{\max }{f}_c>c^2$. For radio communications, it is practically impossible since $c=3\times {10}^8$ m/s and $f_c$ is on the order of ${10}^6\sim {10}^7$ Hz. However, this can be true for UWA channel with $c=1.5\times {10}^3$ m/s and $f_c={10}^2\sim {10}^5$ Hz.

Please note that under the overspread assumption, i.e., ${\tau }_{\max }{v}_{\max }>1$, we have that $(2M+1)/T·L/W>1\Rightarrow L(2M+1)>WT=N$. This means that the estimation problem of ${\mathit{S}}_h$ in (5) is under-determined. Thus, the classical LS estimator, i.e.,

$$\begin{array}{c}\hfill {\widehat{S}}_h=\arg \underset{{\mathit{S}}_h}{min}{\parallel \mathit{y}-\mathit{A}\mathrm{vec}({\mathit{S}}_h)\parallel }_2^2={\mathit{A}}^{†}\mathit{y}\end{array}$$

(6)

may not give a satisfactory performance [20].

Numerous studies tackle the under-determined problem by considering the sparse property of the channel, where the technique of CS has been successfully applied. However, the CS-based estimator relies highly on the sparsity of the channel and the overspread channel may be of poor sparsity in delay-Doppler domain. Denoting ℑ as the support of $S_h(\tau ,v)$, i.e., $S_h(\tau ,v)=0$ for $(\tau ,v)\ne \Im$, then the area of its support $\left|\Im \right|$ can be used as a measurement of the dispersion spread of $S_h(\tau ,v)$. The overspread property of the underwater acoustic channel expresses the fact that the channel’s dispersion spread $\left|\Im \right|$ is relatively large. Hence the sparsity assumption may not directly hold for $S_h(\tau ,v)$ for overspread UWA channels, which motivates us to look for other criterion to estimate $S_h(\tau ,v)$ in case of large dispersion spread.

4. Proposed Method

4.1. Low-Rank-Based Estimator

Although the strong dispersive overspread UWA channels are not of good sparsity, they can enjoy a low-rank or approximately low-rank structure. Consider a doubly spread UWA channel under rich multi-path scattering environment, but tends to be dominated by a relatively small number of clusters of significant paths. Specifically, P dominant paths (clusters) are assumed with each dominant cluster formulated by many non-resolvable sub-paths which give rise to a relatively large dispersion spread. Please note that the well-known leakage effect [11] can also leads to a dispersion spread.

To show the low-rank structure of the UWA channel, we first consider the case that the delay and Doppler are separable in the spreading function, then $S_h(\tau ,v)$ can be expressed as

$$S_h(\tau ,v)=\sum _{p=1}^P{f}_p(\tau )g_p(v)$$

(7)

where $f_p(\tau )$ and $g_p(v)$ describe the delay and Doppler profile of the pth dominant path, respectively. For sparse channel models considered in the literature, $f_p(\tau )$ and $g_p(v)$ are usually assumed to be of a ’sharp’ shape function with a small spread, such as the Dirac delta function or sinc function [10,11,16]. In this paper, we did not have such restriction, and the model in (9) generally takes the case of large dispersion spreads into account.

Following the discrete sampling model in (4), we can rewrite the DDSM as

$${\mathit{S}}_h=\sum _{p=1}^P{\mathit{f}}_p{\mathit{g}}_p^T$$

(8)

where vector ${\mathit{f}}_p={[f_p(0),\dots ,f_p((L-1)\Delta \tau )]}^T\in {\mathbb{C}}^L$, and ${\mathit{g}}_p={[g_p(-M\Delta v),\dots ,g_p(0),\dots ,g_p(M\Delta v)]}^T\in {\mathbb{C}}^{2M+1}$. Then it is obvious that ${\mathit{f}}_p{\mathit{g}}_p^T$ is a rank-one matrix, and thus the spreading function matrix ${\mathit{S}}_h$ is a low-rank matrix with its rank no more than P. i.e.,

$$\mathrm{rank}({\mathit{S}}_h)\le P$$

Please note that the low-rank property of ${\mathit{S}}_h$ is independent of the shape and dispersion spread of the function $f_p(\tau )$ and $g_p(v)$ in this case.

In practice, when the spreading function may not perfectly separable in delay and Doppler domain, then the matrix ${\mathit{f}}_p{\mathit{g}}_p^T$ is generally not a strictly rank-one matrix. However, since the dispersion due to scattering or leakage effect tends to independent in delay and Doppler domain, the DDSM can still be well approximated by a low-rank matrix. We will further show the approximate low-rank property through numerical examples in the simulation section. Moreover, since

$$S_h(\tau ,v)=\sum _{({\tau }_i,v_j)\in \Im }{S}_h(\tau ,v)\delta (\tau -{\tau }_i)\delta (v-v_j)$$

(9)

then it is easy to verify

$$\mathrm{rank}({\mathit{S}}_h)\le \left|\Im \right|$$

(10)

which means that the rank of ${\mathit{S}}_h$ is no more than the sparsity of ${\mathit{S}}_h$. In fact, rank describes the 2-D structured sparsity of a matrix, which exploits the correlation between columns and rows within a matrix. While for the classic sparsity measurement, aka, $l_0$ or $l_1$ norm, is essentially a 1-D sparsity measurement which ignores the 2-D structure in the original matrix. Since the UWA channel is known to have a clustering structure and the DDSM is 2-D in nature, we propose to use low-rank criterion to estimate the DDSM. Please note that for certain extreme cases where the UWA is not low-rank, however, is not the scope of this paper.

To exploit the low-rank property, one can obtain an estimate of ${\mathit{S}}_h$, denoted as ${\widehat{\mathit{S}}}_h$. Mathematically, it can be expressed as

$$\begin{array}{cc}\hfill {\widehat{\mathit{S}}}_h& =\arg \underset{{\mathit{S}}_h}{min}\mathrm{rank}({\mathit{S}}_h)\hfill \\ \hfill s.t.& \parallel \mathit{y}-\mathit{A}\mathrm{vec}({\mathit{S}}_h){\parallel }_2^2\le {\epsilon }^2\hfill \end{array}$$

(11)

where $\epsilon$ measures signal reconstruction errors.

Although the optimization problem (11) is simple in form, it is difficult to be solved due to the non-convexity and discrete nature of the rank function. Fortunately, many techniques have been developed for the low-rank matrix recovery framework.Following the idea of SVP framework [21,22], we reformulate the low-rank overspread channel estimation problem as

$$\begin{array}{cc}\hfill {\widehat{\mathit{S}}}_h& =\arg \underset{{\mathit{S}}_h}{min}\{J(\mathrm{vec}({\mathit{S}}_{\mathit{h}}))=\parallel \mathit{y}-\mathit{A}\mathrm{vec}({\mathit{S}}_h){\parallel }_2^2\},s.t.\mathrm{rank}({\mathit{S}}_{\mathit{h}})\le P\hfill \end{array}$$

(12)

The solution of (12) can be achieved by SVP at each iteration. In particular, in the ith iteration, the current result ${{\mathit{S}}_{\mathit{h}}}^i$ is projected onto a low-rank matrix ${\mathit{S}}_{\mathit{h}}^{\mathit{i}*}$, which is defined as ${\mathit{S}}_{\mathit{h}}^{\mathit{i}*}=\mathrm{SVP}({\mathit{S}}_{\mathit{h}}^{\mathit{i}})={\sum }_{p=1}^P{\mathit{u}}_p{\mathit{s}}_p{\mathit{v}}_p$, where ${\left\{{\mathit{s}}_p\right\}}_{p=1}^P$ is the P most significant singular values of ${\mathit{S}}_{\mathit{h}}^{\mathit{i}}$. The SVP algorithm is shown in Algorithm 1 [21,22].

Algorithm 1: The SVP algorithm.

Input: $\mathit{y}$, $\mathit{A}$ Initialization: P, ${\mathit{S}}_{\mathit{h}}^0=\mathbf{0}$, i = 1 Repeat: 1.$\mu \leftarrow -\frac{{\mathit{d}}_{i-1}^H{\mathit{d}}_i}{{\mathit{d}}_i^H(2{\mathit{A}}^H\mathit{A}){\mathit{d}}_i},\mathbf{where}{\mathit{d}}_i\mathbf{is}\mathbf{the}\mathbf{gradient}\mathbf{of}J(\mathbf{vec}({\mathit{S}}_{\mathit{h}}^{\mathit{i}}))$ 2.$\mathbf{vec}({{\mathit{S}}_{\mathit{h}}}^i)\leftarrow \mathbf{vec}({{\mathit{S}}_{\mathit{h}}}^i)+\mu {\mathit{d}}_i$; 3.${{\mathit{S}}_{\mathit{h}}}^i\leftarrow \mathbf{SVP}({{\mathit{S}}_{\mathit{h}}}^i)$; 4.$i\leftarrow i+1$; Until   convergence Output:${\mathit{S}}_{\mathit{h}}$

It should be noted that we only use the SVP method to demonstrate the validity of the proposed low-rank estimator, and it could be replaced by other suitable low-rank matrix recovery techniques.

4.2. Complexity

The computational complexity of the proposed estimator depends on the techniques used for solving the low-rank matrix recovery problem. For example, the SVP solver approximately has a complexity order of $O((\max (L,2M+1)P^2+N_0^{2}L(2M+1))N_{iter})$, where $N_{iter}$ is the number of iterations [22]. For sparsity-based methods, for instance the OMP method, the complexity is $O((N_{0}L(2M+1)+N_0)N_{iter})$ [15]. We can see that the low-rank-based method is generally has more computational burden than the sparsity-based methods. However, as we will see in the simulation section, the low-rank-based method may outperform the sparsity-based method in case of overspread UWA channels.

5. Simulation Results and Performance Comparison

In this section, we illustrate the performance of the proposed method through numerical examples. The time-varying UWA channel is generated with statistical underwater acoustic channel model proposed in [1]. Specifically, the time-domain channel impulse response $h(t,\tau )$ is generated with an observation period of $T=2$ s and sampling rate of $W=1.6$ KHz. Then we can obtain the spreading function $S_h(v,\tau )$ by applying the Fourier transform of $h(t,\tau )$ with respect to t. We consider $S_h(v,\tau )$ is limited in delay domain as $\tau \in [0,50]$ ms and Doppler domain as $v\in [-20,21]$ Hz, which corresponds to L = 50 ms × 1.6 KHz = 80 and 2M + 1 = 2 s × 41 Hz = 81.

Both the underspread (sparse) and overspread (strong dispersion) channel models are considered in the first example. The DDSMs of the UWA channel are shown in Figure 1a,b (with Doppler and time delay normalized with respect to the frequency sampling rate $\Delta v=1/T$ and $\Delta \tau =1/W$, respectively). To show the low-rank property of the channel, we plot the singular values of DDSM (normalized with the maximum value) as shown in Figure 2. For comparison, we also show the sparsity of DDSM by plotting the absolute values of each element in DDSM in descending order (normalized with the maximum value). We can see that the number of dominant singular values of ${\mathit{S}}_h$ is smaller than that of its elements, especially in case of the strong disperse channel.

In the second example, we estimate both the sparse and dispersion UWA channel through a training signal length of $N=WT=3200$. The signal-to-noise ratio (SNR) is 20 dB. Besides the proposed low-rank estimator, we also consider the classic LS estimator and the sparsity-based methods for performance comparison. The orthogonal matching pursuit (OMP) algorithm [12,23,24] and the alternating direction method of multipliers (ADMM) [25] is adopted to solve the sparsity-based estimator. Figure 3 and Figure 4 show the 2-D plots of the estimated DDSM of by different methods under sparse and dispersion channel cases, respectively.

We can see from Figure 3 and Figure 4 that the LS-based method leads to many pseudo paths as well as limited accuracy under both sparse and dispersion channel cases. That is because LS method does not exploit either sparsity or low-rank property of the DDSM. The OMP and ADMM methods work well for the sparse channel case due to its sparsity exploitation. For the dispersion case, however, they suffer from performance degradation due to the poor sparsity of the overspread channel. The proposed low-rank estimator can achieve a comparable performance for the sparse channel case, and better performance than the compared methods for the dispersion channel case.

In the third example, we compare the normalized mean squared error (NMSE) performance of the considered estimators. The NMSE is defined as

$$\mathrm{NMSE}=E\left\{\frac{\parallel {\mathit{S}}_h-{\widehat{\mathit{S}}}_h{\parallel }_F^2}{\parallel {\mathit{S}}_h{\parallel }_F^2}\right\}.$$

The NMSE performance versus different SNR is shown in Figure 5. We can see that the proposed low-rank-based estimator has similar performance to the ADMM method in case of sparse channel, and outperforms the two sparsity-based estimators and the LS estimator in case of dispersive channel under relatively higher SNRs.

In the fourth example, we investigate the estimation performance versus signal length N. The SNR is fixed to 20 dB, and N varies from 500 to 3000. As shown in Figure 6, the proposed Low-rank-based method again outperforms the other compared methods in case of dispersive channel under moderate or large number of training signal.

Finally, we investigate the impact of user parameters in the proposed method under. The key parameter in SVP algorithm is the rank projection parameter P. The NMSE performance of the proposed method versus parameter P is shown in Figure 7. We can see that the proposed method achieves the best performance when P is near 5 to 9, which is approximately the number of the significant singular values of the DDSM. Moreover, we can also indicate from Figure 7 that it is better to overestimate rather than underestimate P when the exact rank of the DDSM is not known in practice.

6. Conclusions

In this paper, the estimation of a dispersive overspread UWA channel is considered. We formulate the channel input–output relationship in a discrete delay-Doppler-spread function representation. Then we show that although the channel dispersion causes the degeneration of channel sparsity, it gives rise to a useful low-rank structure, when the channel delay-Doppler-spread function is separable in delay and Doppler domain. This low-rank structure reveals that the matrix representation of delay-Doppler-spread function can be determined by only a small number of singular values. Therefore, we introduce the low-rank criterion to estimate the overspread channels, which can help to improve the estimation performance in case of strong dispersive UWA channels. The estimator is formulated in low-rank matrix recovery problem and is solved by the SVP technique. Simulation examples are carried out to demonstrate the effectiveness of the proposed low-rank-based channel estimator. In future work, it is interesting to further consider a joint sparse and low-rank criterion for UWA channel estimation.